33 #define MU 398600.8 // Earth gravitational constant (km3/s2)
34 #define RADIUSEARTHKM 6378.135 // Earth radius (km)
35 #define XKE 0.07436691613317 // 60.0 / sqrt(RADIUSEARTHKM^3/MU)
36 #define TUMIN 13.44683969695931 // sqrt(RADIUSEARTHKM^3/MU) / 60.0
37 #define J2 0.001082616 // The second gravitational zonal harmonic of the Earth
38 #define J3 -0.00000253881 // The third gravitational zonal harmonic of the Earth
39 #define J4 -0.00000165597 // The fourth gravitational zonal harmonic of the Earth
40 #define J3OJ2 -2.34506972e-3 // J3 / J2
43 #define PI 3.14159265358979323846
44 #define TWOPI 6.2831853071795864769 // 2*PI
45 #define PIO2 1.5707963267948966192 // PI/2
46 #define X3PIO2 4.7123889803846898577 // 3*PI/2
47 #define X2O3 .66666666666666666667 // 2/3
48 #define DEG2RAD 1.745329251994330e-2 // Deg -> Rad
51 #define MINPD 1440 // Minutes per day
52 #define MEANALT 0.84 // Mean altitude (km)
53 #define SR 6.96000e5 // Solar radius - km (IAU 76)
54 #define AU 1.49597870691e8 // Astronomical unit - km (IAU 76)
55 #define XPDOTP 229.1831180523293 // 1440.0 / (2.0 * pi)
56 #define SS 1.0122292801892716288 // Parameter for the SGP4 density function
57 #define QZMS2T 1.8802791590152706439e-9 // (( 120.0 - 78.0) / RADIUSEARTHKM )^4
58 #define F 3.35281066474748e-3 // Flattening factor
59 #define MFACTOR 7.292115e-5
65 m_number = line1.mid( 2, 5 ).toInt();
66 m_class = line1.at( 7 );
67 m_id = line1.mid( 9, 8 );
68 m_epoch = line1.mid( 18, 14 ).toDouble();
69 m_first_deriv = line1.mid( 33, 10 ).toDouble() / (
XPDOTP *
MINPD );
70 m_second_deriv = line1.mid( 44, 6 ).toDouble() * ( 1.0e-5 / pow( 10.0, line1.mid( 51, 1 ).toDouble() ) ) / (
XPDOTP *
MINPD*
MINPD );
71 m_bstar = line1.mid( 53, 6 ).toDouble() * 1.0e-5 / pow( 10.0, line1.mid( 60, 1 ).toDouble() );
72 m_ephem_type = line1.mid( 62, 1 ).toInt();
73 m_elem_number = line1.mid( 64, 4 ).toInt();
74 m_inclination = line2.mid( 8, 8 ).toDouble() *
DEG2RAD;
75 m_ra = line2.mid( 17, 8 ).toDouble() *
DEG2RAD;
76 m_eccentricity = line2.mid( 26, 7 ).toDouble() * 1.0e-7;
77 m_arg_perigee = line2.mid( 34, 8 ).toDouble() *
DEG2RAD;
78 m_mean_anomaly = line2.mid( 43, 8 ).toDouble() *
DEG2RAD;
79 m_mean_motion = line2.mid( 52, 11 ).toDouble() *
TWOPI /
MINPD;
80 m_nb_revolution = line2.mid( 63, 5 ).toInt();
91 double day = modf( m_epoch * 1.e-3, &m_epoch_year) * 1.e3;
92 if ( m_epoch_year < 57. )
93 m_epoch_year += 2000.;
95 m_epoch_year += 1900.;
96 double year= m_epoch_year - 1.;
103 m_tle_jd = i + 1720994.5 + B + day;
113 void Satellite::init()
115 double ao , cosio , sinio , cosio2,
116 omeosq, posq , rp , rteosq, eccsq , con42 ,
117 cnodm , snodm , cosim , sinim , cosomm, sinomm, cc1sq ,
118 cc2 , cc3 , coef , coef1 , cosio4, day , dndt ,
119 em , emsq , eeta , etasq , gam ,
120 inclm , nm , perige, pinvsq, psisq , qzms24,
121 rtemsq, s1 , s2 , s3 , s4 , s5 , s6 ,
122 s7 , sfour , ss1(0), ss2(0), ss3(0), ss4(0), ss5(0),
123 ss6(0), ss7(0), sz1(0), sz2(0), sz3(0), sz11(0), sz12(0),
124 sz13(0), sz21(0), sz22(0), sz23(0), sz31(0), sz32(0), sz33(0),
125 tc , temp , temp1 , temp2 , temp3 , tsi , xpidot,
126 xhdot1, z1 , z2 , z3 , z11 , z12 , z13 ,
127 z21 , z22 , z23 , z31 , z32 , z33 , ak ,
128 d1 , del , adel , po , ds70 , ts70 , tfrac ,
129 c1 , thgr70, fk5r , c1p2p;
133 aycof =0. ; con41 =0. ; cc1 =0. ; cc4 =0. ; cc5 =0. ; d2 =0. ; d3 =0. ; d4 =0. ;
134 delmo =0. ; eta =0. ; argpdot =0. ; omgcof =0. ; sinmao =0. ; t =0. ; t2cof =0. ; t3cof =0. ;
135 t4cof =0. ; t5cof =0. ; x1mth2 =0. ; x7thm1 =0. ; mdot =0. ; nodedot =0. ; xlcof =0. ; xmcof =0. ;
140 d2201 =0. ; d2211 =0. ; d3210 =0. ; d3222 =0. ; d4410 =0. ; d4422 =0. ; d5220 =0. ; d5232 =0. ;
141 d5421 =0. ; d5433 =0. ; dedt =0. ; del1 =0. ; del2 =0. ; del3 =0. ; didt =0. ; dmdt =0. ;
142 dnodt =0. ; domdt =0. ; e3 =0. ; ee2 =0. ; peo =0. ; pgho =0. ; pho =0. ; pinco =0. ;
143 plo =0. ; se2 =0. ; se3 =0. ; sgh2 =0. ; sgh3 =0. ; sgh4 =0. ; sh2 =0. ; sh3 =0. ;
144 si2 =0. ; si3 =0. ; sl2 =0. ; sl3 =0. ; sl4 =0. ; gsto =0. ; xfact =0. ; xgh2 =0. ;
145 xgh3 =0. ; xgh4 =0. ; xh2 =0. ; xh3 =0. ; xi2 =0. ; xi3 =0. ; xl2 =0. ; xl3 =0. ;
146 xl4 =0. ; xlamo =0. ; zmol =0. ; zmos =0. ; atime =0. ; xli =0. ; xni =0.;
150 m_is_visible =
false;
153 const double temp4 = 1.5e-12;
157 eccsq = m_eccentricity * m_eccentricity;
158 omeosq = 1.0 - eccsq;
159 rteosq = sqrt( omeosq );
160 cosio = cos( m_inclination );
161 cosio2 = cosio * cosio;
164 ak = pow(
XKE / m_mean_motion,
X2O3 );
165 d1 = 0.75 *
J2 * ( 3.0 * cosio2 - 1.0 ) / ( rteosq * omeosq );
166 del = d1 / ( ak * ak );
167 adel = ak * ( 1.0 - del*del - del * ( 1.0 / 3.0 + 134.0 * del*del / 81.0 ) );
168 del = d1 / ( adel * adel );
169 m_mean_motion = m_mean_motion / (1.0 + del);
171 ao = pow(
XKE / m_mean_motion,
X2O3 );
172 sinio = sin( m_inclination );
174 con42 = 1.0 - 5.0 * cosio2;
175 con41 = -con42 - ( 2.0 * cosio2 );
177 rp = ao * ( 1.0 - m_eccentricity );
181 ts70 = m_tle_jd - 2433281.5 - 7305.0;
182 ds70 = floor( ts70 + 1.0e-8 );
185 c1 = 1.72027916940703639e-2;
186 thgr70= 1.7321343856509374;
187 fk5r = 5.07551419432269442e-15;
189 gsto = fmod( thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r,
TWOPI );
193 if ( ( omeosq >= 0.0 ) || ( m_mean_motion >= 0.0 ) ) {
201 if ( perige < 156.0 ) {
202 sfour = perige - 78.0;
210 tsi = 1.0 / ( ao - sfour );
211 eta = ao * m_eccentricity * tsi;
213 eeta = m_eccentricity * eta;
214 psisq = fabs( 1.0 - etasq );
215 coef = qzms24 * pow( tsi, 4.0 );
216 coef1 = coef / pow( psisq, 3.5 );
217 cc2 = coef1 * m_mean_motion * ( ao * ( 1.0 + 1.5 * etasq + eeta *
218 ( 4.0 + etasq ) ) + 0.375 *
J2 * tsi / psisq * con41 *
219 ( 8.0 + 3.0 * etasq * ( 8.0 + etasq ) ) );
222 if ( m_eccentricity > 1.0e-4 )
223 cc3 = -2.0 * coef * tsi *
J3OJ2 * m_mean_motion * sinio / m_eccentricity;
224 x1mth2 = 1.0 - cosio2;
225 cc4 = 2.0 * m_mean_motion * coef1 * ao * omeosq *
226 ( eta * ( 2.0 + 0.5 * etasq ) + m_eccentricity *
227 ( 0.5 + 2.0 * etasq ) -
J2 * tsi / ( ao * psisq ) *
228 ( -3.0 * con41 * ( 1.0 - 2.0 * eeta + etasq *
229 ( 1.5 - 0.5 * eeta ) ) + 0.75 * x1mth2 *
230 ( 2.0 * etasq - eeta * ( 1.0 + etasq ) ) * cos( 2.0 * m_arg_perigee ) ) );
231 cc5 = 2.0 * coef1 * ao * omeosq * ( 1.0 + 2.75 * ( etasq + eeta ) + eeta * etasq );
233 cosio4 = cosio2 * cosio2;
234 temp1 = 1.5 *
J2 * pinvsq * m_mean_motion;
235 temp2 = 0.5 * temp1 *
J2 * pinvsq;
236 temp3 = -0.46875 *
J4 * pinvsq * pinvsq * m_mean_motion;
237 mdot = m_mean_motion + 0.5 * temp1 * rteosq * con41 + 0.0625 *
238 temp2 * rteosq * ( 13.0 - 78.0 * cosio2 + 137.0 * cosio4 );
239 argpdot= -0.5 * temp1 * con42 + 0.0625 * temp2 *
240 ( 7.0 - 114.0 * cosio2 + 395.0 * cosio4 ) +
241 temp3 * ( 3.0 - 36.0 * cosio2 + 49.0 * cosio4 );
242 xhdot1 = -temp1 * cosio;
243 nodedot= xhdot1 + ( 0.5 * temp2 * ( 4.0 - 19.0 * cosio2 ) +
244 2.0 * temp3 * ( 3.0 - 7.0 * cosio2 ) ) * cosio;
245 xpidot = argpdot + nodedot;
246 omgcof = m_bstar * cc3 * cos( m_arg_perigee );
248 if ( m_eccentricity > 1.0e-4)
249 xmcof = -
X2O3 * coef * m_bstar / eeta;
250 nodecf = 3.5 * omeosq * xhdot1 * cc1;
253 if ( fabs( 1.0 + cosio ) > 1.5e-12 )
254 xlcof = -0.25 *
J3OJ2 * sinio * ( 3.0 + 5.0 * cosio ) / ( 1.0 + cosio );
256 xlcof = -0.25 *
J3OJ2 * sinio * ( 3.0 + 5.0 * cosio ) / temp4;
257 aycof = -0.5 *
J3OJ2 * sinio;
258 delmo = pow( ( 1.0 + eta * cos( m_mean_anomaly ) ), 3 );
259 sinmao = sin( m_mean_anomaly );
260 x7thm1 = 7.0 * cosio2 - 1.0;
263 if ( ( TWOPI / m_mean_motion ) >= 225.0 ) {
267 inclm = m_inclination;
271 const double zes = 0.01675;
272 const double zel = 0.05490;
273 const double c1ss = 2.9864797e-6;
274 const double c1l = 4.7968065e-7;
275 const double zsinis = 0.39785416;
276 const double zcosis = 0.91744867;
277 const double zcosgs = 0.1945905;
278 const double zsings = -0.98088458;
281 double a1 , a2 , a3 , a4 , a5 , a6 , a7 ,
282 a8 , a9 , a10 , betasq, cc , ctem , stem ,
283 x1 , x2 , x3 , x4 , x5 , x6 , x7 ,
284 x8 , xnodce, xnoi , zcosg , zcosgl, zcosh , zcoshl,
285 zcosi , zcosil, zsing , zsingl, zsinh , zsinhl, zsini ,
292 sinomm = sin( m_arg_perigee );
293 cosomm = cos( m_arg_perigee );
294 sinim = sin( m_inclination );
295 cosim = cos( m_inclination );
298 rtemsq = sqrt( betasq );
306 day = m_tle_jd - 2433281.5 + 18261.5 + tc / 1440.0;
307 xnodce = fmod( 4.5236020 - 9.2422029e-4 * day, TWOPI );
308 stem = sin( xnodce );
309 ctem = cos( xnodce );
310 zcosil = 0.91375164 - 0.03568096 * ctem;
311 zsinil = sqrt( 1.0 - zcosil * zcosil );
312 zsinhl = 0.089683511 * stem / zsinil;
313 zcoshl = sqrt( 1.0 - zsinhl * zsinhl );
314 gam = 5.8351514 + 0.0019443680 * day;
315 zx = 0.39785416 * stem / zsinil;
316 zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
317 zx = atan2( zx, zy );
318 zx = gam + zx - xnodce;
332 for ( lsflg = 1; lsflg <= 2; ++lsflg ) {
333 a1 = zcosg * zcosh + zsing * zcosi * zsinh;
334 a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
335 a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
337 a9 = zsing * zsinh + zcosg * zcosi * zcosh;
339 a2 = cosim * a7 + sinim * a8;
340 a4 = cosim * a9 + sinim * a10;
341 a5 = -sinim * a7 + cosim * a8;
342 a6 = -sinim * a9 + cosim * a10;
344 x1 = a1 * cosomm + a2 * sinomm;
345 x2 = a3 * cosomm + a4 * sinomm;
346 x3 = -a1 * sinomm + a2 * cosomm;
347 x4 = -a3 * sinomm + a4 * cosomm;
353 z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3;
354 z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4;
355 z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4;
356 z1 = 3.0 * ( a1 * a1 + a2 * a2 ) + z31 * emsq;
357 z2 = 6.0 * ( a1 * a3 + a2 * a4 ) + z32 * emsq;
358 z3 = 3.0 * ( a3 * a3 + a4 * a4 ) + z33 * emsq;
359 z11 = -6.0 * a1 * a5 + emsq * ( -24.0 * x1 * x7-6.0 * x3 * x5 );
360 z12 = -6.0 * ( a1 * a6 + a3 * a5 ) + emsq *
361 ( -24.0 * ( x2 * x7 + x1 * x8 ) - 6.0 * ( x3 * x6 + x4 * x5 ) );
362 z13 = -6.0 * a3 * a6 + emsq * ( -24.0 * x2 * x8 - 6.0 * x4 * x6 );
363 z21 = 6.0 * a2 * a5 + emsq * ( 24.0 * x1 * x5 - 6.0 * x3 * x7 );
364 z22 = 6.0 * ( a4 * a5 + a2 * a6 ) + emsq *
365 ( 24.0 * ( x2 * x5 + x1 * x6 ) - 6.0 * ( x4 * x7 + x3 * x8 ) );
366 z23 = 6.0 * a4 * a6 + emsq * ( 24.0 * x2 * x6 - 6.0 * x4 * x8 );
367 z1 = z1 + z1 + betasq * z31;
368 z2 = z2 + z2 + betasq * z32;
369 z3 = z3 + z3 + betasq * z33;
371 s2 = -0.5 * s3 / rtemsq;
373 s1 = -15.0 * em * s4;
374 s5 = x1 * x3 + x2 * x4;
375 s6 = x2 * x3 + x1 * x4;
376 s7 = x2 * x4 - x1 * x3;
403 zcosh = zcoshl * cnodm + zsinhl * snodm;
404 zsinh = snodm * zcoshl - cnodm * zsinhl;
409 zmol = fmod( 4.7199672 + 0.22997150 * day - gam, TWOPI );
410 zmos = fmod( 6.2565837 + 0.017201977 * day, TWOPI );
413 se2 = 2.0 * ss1 * ss6;
414 se3 = 2.0 * ss1 * ss7;
415 si2 = 2.0 * ss2 * sz12;
416 si3 = 2.0 * ss2 * ( sz13 - sz11 );
417 sl2 = -2.0 * ss3 * sz2;
418 sl3 = -2.0 * ss3 * ( sz3 - sz1 );
419 sl4 = -2.0 * ss3 * ( -21.0 - 9.0 * emsq ) * zes;
420 sgh2 = 2.0 * ss4 * sz32;
421 sgh3 = 2.0 * ss4 * ( sz33 - sz31 );
422 sgh4 = -18.0 * ss4 * zes;
423 sh2 = -2.0 * ss2 * sz22;
424 sh3 = -2.0 * ss2 * ( sz23 - sz21 );
429 xi2 = 2.0 * s2 * z12;
430 xi3 = 2.0 * s2 * ( z13 - z11 );
431 xl2 = -2.0 * s3 * z2;
432 xl3 = -2.0 * s3 * ( z3 - z1 );
433 xl4 = -2.0 * s3 * ( -21.0 - 9.0 * emsq ) * zel;
434 xgh2 = 2.0 * s4 * z32;
435 xgh3 = 2.0 * s4 * ( z33 - z31 );
436 xgh4 = -18.0 * s4 * zel;
437 xh2 = -2.0 * s2 * z22;
438 xh3 = -2.0 * s2 * ( z23 - z21 );
442 ses , sghl , sghs , shll, shs,
443 sinzf, sis , sls , zf , zm;
446 const double zns = 1.19459e-5;
447 const double znl = 1.5835218e-4;
451 zf = zm + 2.0 * zes * sin( zm );
453 f2 = 0.5 * sinzf * sinzf - 0.25;
454 f3 = -0.5 * sinzf * cos( zf );
455 ses = se2* f2 + se3 * f3;
456 sis = si2 * f2 + si3 * f3;
457 sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
458 sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
459 shs = sh2 * f2 + sh3 * f3;
462 zf = zm + 2.0 * zel * sin( zm );
464 f2 = 0.5 * sinzf * sinzf - 0.25;
465 f3 = -0.5 * sinzf * cos( zf );
466 sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
467 shll = xh2 * f2 + xh3 * f3;
470 double ainv2 , aonv=0.0, cosisq, eoc , f220 , f221 , f311 ,
471 f321 , f322 , f330 , f441 , f442 , f522 , f523 ,
472 f542 , f543 , g200 , g201 , g211 , g300 , g310 ,
473 g322 , g410 , g422 , g520 , g521 , g532 , g533 ,
474 sgs , sini2 , temp , temp1 , theta , xno2 , emo ,
478 const double q22 = 1.7891679e-6;
479 const double q31 = 2.1460748e-6;
480 const double q33 = 2.2123015e-7;
481 const double root22 = 1.7891679e-6;
482 const double root44 = 7.3636953e-9;
483 const double root54 = 2.1765803e-9;
484 const double rptim = 4.37526908801129966e-3;
485 const double root32 = 3.7393792e-7;
486 const double root52 = 1.1428639e-7;
490 if ( ( nm < 0.0052359877 ) && ( nm > 0.0034906585 ) )
492 if ( ( nm >= 8.26e-3 ) && ( nm <= 9.24e-3 ) && ( em >= 0.5 ) )
496 ses = ss1 * zns * ss5;
497 sis = ss2 * zns * ( sz11 + sz13 );
498 sls = -zns * ss3 * ( sz1 + sz3 - 14.0 - 6.0 * emsq );
499 sghs = ss4 * zns * ( sz31 + sz33 - 6.0 );
500 shs = -zns * ss2 * ( sz21 + sz23 );
501 if ( ( inclm < 5.2359877e-2 ) || ( inclm >
PI - 5.2359877e-2 ) )
505 sgs = sghs - cosim * shs;
508 dedt = ses + s1 * znl * s5;
509 didt = sis + s2 * znl * ( z11 + z13 );
510 dmdt = sls - znl * s3 * ( z1 + z3 - 14.0 - 6.0 * emsq );
511 sghl = s4 * znl * ( z31 + z33 - 6.0 );
512 shll = -znl * s2 * ( z21 + z23 );
513 if ( ( inclm < 5.2359877e-2 ) || ( inclm >
PI - 5.2359877e-2 ) )
517 if ( sinim != 0.0 ) {
518 domdt = domdt - cosim / sinim * shll;
519 dnodt = dnodt + shll / sinim;
524 theta = fmod( gsto + tc * rptim, TWOPI );
532 cosisq = cosim * cosim;
538 g201 = -0.306 - ( em - 0.64 ) * 0.440;
541 g211 = 3.616 - 13.2470 * em + 16.2900 * emsq;
542 g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc;
543 g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc;
544 g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc;
545 g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc;
546 g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc;
548 g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc;
549 g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc;
550 g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc;
551 g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc;
552 g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc;
554 g520 =-5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc;
556 g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq;
559 g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc;
560 g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc;
561 g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc;
563 g533 =-37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc;
564 g521 =-51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc;
565 g532 =-40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc;
568 sini2= sinim * sinim;
569 f220 = 0.75 * ( 1.0 + 2.0 * cosim + cosisq );
571 f321 = 1.875 * sinim * ( 1.0 - 2.0 * cosim - 3.0 * cosisq );
572 f322 = -1.875 * sinim * ( 1.0 + 2.0 * cosim - 3.0 * cosisq );
573 f441 = 35.0 * sini2 * f220;
574 f442 = 39.3750 * sini2 * sini2;
575 f522 = 9.84375 * sinim * (sini2 * ( 1.0 - 2.0 * cosim- 5.0 * cosisq ) +
576 0.33333333 * ( -2.0 + 4.0 * cosim + 6.0 * cosisq ) );
577 f523 = sinim * ( 4.92187512 * sini2 * ( -2.0 - 4.0 * cosim +
578 10.0 * cosisq ) + 6.56250012 * ( 1.0+2.0 * cosim - 3.0 * cosisq ) );
579 f542 = 29.53125 * sinim * ( 2.0 - 8.0 * cosim + cosisq *
580 ( -12.0 + 8.0 * cosim + 10.0 * cosisq ) );
581 f543 = 29.53125 * sinim * ( -2.0 - 8.0 * cosim + cosisq *
582 ( 12.0 + 8.0 * cosim - 10.0 * cosisq ) );
585 temp1 = 3.0 * xno2 * ainv2;
586 temp = temp1 * root22;
587 d2201 = temp * f220 * g201;
588 d2211 = temp * f221 * g211;
589 temp1 = temp1 * aonv;
590 temp = temp1 * root32;
591 d3210 = temp * f321 * g310;
592 d3222 = temp * f322 * g322;
593 temp1 = temp1 * aonv;
594 temp = 2.0 * temp1 * root44;
595 d4410 = temp * f441 * g410;
596 d4422 = temp * f442 * g422;
597 temp1 = temp1 * aonv;
598 temp = temp1 * root52;
599 d5220 = temp * f522 * g520;
600 d5232 = temp * f523 * g532;
601 temp = 2.0 * temp1 * root54;
602 d5421 = temp * f542 * g521;
603 d5433 = temp * f543 * g533;
604 xlamo = fmod( m_mean_anomaly + m_ra + m_ra - theta - theta, TWOPI );
605 xfact = mdot + dmdt + 2.0 * ( nodedot + dnodt - rptim ) - m_mean_motion;
611 g200 = 1.0 + emsq * ( -2.5 + 0.8125 * emsq );
612 g310 = 1.0 + 2.0 * emsq;
613 g300 = 1.0 + emsq * ( -6.0 + 6.60937 * emsq );
614 f220 = 0.75 * ( 1.0 + cosim ) * ( 1.0 + cosim );
615 f311 = 0.9375 * sinim * sinim * ( 1.0 + 3.0 * cosim ) - 0.75 * ( 1.0 + cosim );
617 f330 = 1.875 * f330 * f330 * f330;
618 del1 = 3.0 * nm * nm * aonv * aonv;
619 del2 = 2.0 * del1 * f220 * g200 * q22;
620 del3 = 3.0 * del1 * f330 * g300 * q33 * aonv;
621 del1 = del1 * f311 * g310 * q31 * aonv;
622 xlamo = fmod( m_mean_anomaly + m_ra + m_arg_perigee - theta, TWOPI );
623 xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - m_mean_motion;
629 nm = m_mean_motion + dndt;
636 d2 = 4.0 * ao * tsi * cc1sq;
637 temp = d2 * tsi * cc1 / 3.0;
638 d3 = ( 17.0 * ao + sfour ) * temp;
639 d4 = 0.5 * temp * ao * tsi * ( 221.0 * ao + 31.0 * sfour ) * cc1;
640 t3cof = d2 + 2.0 * cc1sq;
641 t4cof = 0.25 * ( 3.0 * d3 + cc1 * ( 12.0 * d2 + 10.0 * cc1sq ) );
642 t5cof = 0.2 * ( 3.0 * d4 + 12.0 * cc1 * d3 + 6.0 * d2 * d2 + 15.0 * cc1sq * ( 2.0 * d2 + cc1sq ) );
653 int Satellite::sgp4(
double tsince )
658 double am , axnl , aynl , betal , cosim , cnod ,
659 cos2u, coseo1, cosi , cosip , cosisq, cossu , cosu,
660 delm , delomg, em , emsq , ecose , el2 , eo1 ,
661 ep , esine , argpm, argpp , argpdf, pl , mrt = 0.0,
662 mvt , rdotl , rl , rvdot , rvdotl, sinim , dndt,
663 sin2u, sineo1, sini , sinip , sinsu , sinu ,
664 snod , su , t2 , t3 , t4 , tem5 , temp,
665 temp1, temp2 , tempa, tempe , templ , u , ux ,
666 uy , uz , vx , vy , vz , inclm , mm ,
667 nm , nodem , xinc , xincp , xl , xlm , mp ,
668 xmdf , xmx , xmy , nodedf, xnode , nodep , tc ,
669 sat_posx, sat_posy , sat_posz, sat_posw, sat_velx ,
670 sat_vely , sat_velz , sinlat, obs_posx, obs_posy,
671 obs_posz, obs_posw, obs_velx, obs_vely, obs_velz,
672 coslat, thetageo, sintheta, costheta, c, sq, achcp,
675 const double temp4 = 1.5e-12;
682 xmdf = m_mean_anomaly + mdot * tsince;
683 argpdf = m_arg_perigee + argpdot * tsince;
684 nodedf = m_ra + nodedot * tsince;
687 t2 = tsince * tsince;
688 nodem = nodedf + nodecf * t2;
689 tempa = 1.0 - cc1 * tsince;
690 tempe = m_bstar * cc4 * tsince;
694 delomg = omgcof * tsince;
695 delm = xmcof * ( pow( ( 1.0 + eta * cos( xmdf ) ), 3 ) - delmo);
696 temp = delomg + delm;
698 argpm = argpdf - temp;
701 tempa = tempa - d2 * t2 - d3 * t3 - d4 * t4;
702 tempe = tempe + m_bstar * cc5 * ( sin( mm ) - sinmao );
703 templ = templ + t3cof * t3 + t4 * ( t4cof + tsince * t5cof );
708 inclm = m_inclination;
710 if ( method ==
'd' ) {
714 double delt, ft, theta, x2li, x2omi, xl, xldot, xnddt, xndt, xomi;
717 const double fasx2 = 0.13130908;
718 const double fasx4 = 2.8843198;
719 const double fasx6 = 0.37448087;
720 const double g22 = 5.7686396;
721 const double g32 = 0.95240898;
722 const double g44 = 1.8014998;
723 const double g52 = 1.0508330;
724 const double g54 = 4.4108898;
725 const double rptim = 4.37526908801129966e-3;
726 const double step = 720.0;
727 const double step2 = step * step / 2;
731 theta = fmod( gsto + tc * rptim, TWOPI );
732 em = em + dedt * tsince;
734 inclm = inclm + didt * tsince;
735 argpm = argpm + domdt * tsince;
736 nodem = nodem + dnodt * tsince;
737 mm = mm + dmdt * tsince;
742 if ( ( atime == 0.0 ) || ( tsince * atime <= 0.0 ) || ( fabs( tsince ) < fabs( atime ) ) ) {
755 while ( iretn == 381 ) {
758 xndt = del1 * sin( xli - fasx2 ) + del2 * sin( 2.0 * ( xli - fasx4 ) ) +
759 del3 * sin( 3.0 * ( xli - fasx6 ) );
761 xnddt = del1 * cos( xli - fasx2 ) +
762 2.0 * del2 * cos( 2.0 * ( xli - fasx4 ) ) +
763 3.0 * del3 * cos( 3.0 * ( xli - fasx6 ) );
764 xnddt = xnddt * xldot;
767 xomi = m_arg_perigee + argpdot * atime;
770 xndt = d2201 * sin( x2omi + xli - g22 ) + d2211 * sin( xli - g22 ) +
771 d3210 * sin( xomi + xli - g32 ) + d3222 * sin( -xomi + xli - g32 )+
772 d4410 * sin( x2omi + x2li - g44 )+ d4422 * sin( x2li - g44 ) +
773 d5220 * sin( xomi + xli - g52 ) + d5232 * sin( -xomi + xli - g52 )+
774 d5421 * sin( xomi + x2li - g54 ) + d5433 * sin( -xomi + x2li - g54 );
776 xnddt = d2201 * cos( x2omi + xli - g22 ) + d2211 * cos( xli - g22 ) +
777 d3210 * cos( xomi + xli - g32 ) + d3222 * cos( -xomi + xli - g32 ) +
778 d5220 * cos( xomi + xli - g52 ) + d5232 * cos( -xomi + xli - g52 ) +
779 2.0 * ( d4410 * cos( x2omi + x2li - g44 ) +
780 d4422 * cos( x2li - g44 ) + d5421 * cos( xomi + x2li - g54 ) +
781 d5433 * cos( -xomi + x2li - g54 ) );
782 xnddt = xnddt * xldot;
785 if ( fabs( tsince - atime ) >= step ) {
792 if ( iretn == 381 ) {
793 xli = xli + xldot * delt + xndt * step2;
794 xni = xni + xndt * delt + xnddt * step2;
795 atime = atime + delt;
799 nm = xni + xndt * ft + xnddt * ft * ft * 0.5;
800 xl = xli + xldot * ft + xndt * ft * ft * 0.5;
803 mm = xl - 2.0 * nodem + 2.0 * theta;
804 dndt = nm - m_mean_motion;
806 mm = xl - nodem - argpm + theta;
807 dndt = nm - m_mean_motion;
810 nm = m_mean_motion + dndt;
815 kDebug() <<
"Mean motion less than 0.0";
819 am = pow( (
XKE / nm ),
X2O3 ) * tempa * tempa;
820 nm =
XKE / pow( am, 1.5 );
823 if ( ( em >= 1.0 ) || ( em < -0.001 ) ) {
824 kDebug() <<
"Eccentricity >= 1.0 or < -0.001";
831 mm = mm + m_mean_motion * templ;
832 xlm = mm + argpm + nodem;
836 nodem = fmod( nodem, TWOPI );
837 argpm = fmod( argpm, TWOPI );
838 xlm = fmod( xlm, TWOPI );
839 mm = fmod( xlm - argpm - nodem, TWOPI );
842 sinim = sin( inclm );
843 cosim = cos( inclm );
853 if ( method ==
'd' ) {
854 double alfdp, betdp, cosip, cosop, dalf, dbet, dls,
855 f2 , f3 , pe , pgh , ph , pinc, pl ,
856 sel , ses , sghl , sghs , shll, shs , sil,
857 sinip, sinop, sinzf, sis , sll , sls , xls,
861 const double zns = 1.19459e-5;
862 const double zes = 0.01675;
863 const double znl = 1.5835218e-4;
864 const double zel = 0.05490;
867 zm = zmos + zns * tsince;
868 zf = zm + 2.0 * zes * sin( zm );
870 f2 = 0.5 * sinzf * sinzf - 0.25;
871 f3 = -0.5 * sinzf * cos( zf );
872 ses = se2* f2 + se3 * f3;
873 sis = si2 * f2 + si3 * f3;
874 sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
875 sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
876 shs = sh2 * f2 + sh3 * f3;
877 zm = zmol + znl * tsince;
879 zf = zm + 2.0 * zel * sin( zm );
881 f2 = 0.5 * sinzf * sinzf - 0.25;
882 f3 = -0.5 * sinzf * cos( zf );
883 sel = ee2 * f2 + e3 * f3;
884 sil = xi2 * f2 + xi3 * f3;
885 sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf;
886 sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
887 shll = xh2 * f2 + xh3 * f3;
899 xincp = xincp + pinc;
901 sinip = sin( xincp );
902 cosip = cos( xincp );
905 if ( xincp >= 0.2 ) {
907 pgh = pgh - cosip * ph;
913 sinop = sin( nodep );
914 cosop = cos( nodep );
915 alfdp = sinip * sinop;
916 betdp = sinip * cosop;
917 dalf = ph * cosop + pinc * cosip * sinop;
918 dbet = -ph * sinop + pinc * cosip * cosop;
919 alfdp = alfdp + dalf;
920 betdp = betdp + dbet;
921 nodep = fmod( nodep, TWOPI );
924 xls = mp + argpp + cosip * nodep;
925 dls = pl + pgh - pinc * nodep * sinip;
928 nodep = atan2( alfdp, betdp );
929 if ( ( nodep < 0.0) )
931 if ( fabs( xnoh - nodep ) >
PI ) {
938 argpp = xls - mp - cosip * nodep;
947 if ( ( ep < 0.0 ) || ( ep > 1.0 ) ) {
948 kDebug() <<
"Eccentricity < 0.0 or > 1.0";
954 if ( method ==
'd' ) {
955 sinip = sin( xincp );
956 cosip = cos( xincp );
957 aycof = -0.5 *
J3OJ2 * sinip;
958 if ( fabs( cosip + 1.0 ) > 1.5e-12 )
959 xlcof = -0.25 *
J3OJ2 * sinip * ( 3.0 + 5.0 * cosip ) / ( 1.0 + cosip );
961 xlcof = -0.25 *
J3OJ2 * sinip * ( 3.0 + 5.0 * cosip ) / temp4;
963 axnl = ep * cos( argpp );
964 temp = 1.0 / ( am * ( 1.0 - ep * ep ) );
965 aynl = ep* sin( argpp ) + temp * aycof;
966 xl = mp + argpp + nodep + temp * xlcof * axnl;
969 u = fmod( xl - nodep, TWOPI );
973 while ( ( fabs( tem5 ) >= 1.0e-12 ) && ( ktr <= 10 ) ) {
976 tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl;
977 tem5 = ( u - aynl * coseo1 + axnl * sineo1 - eo1 ) / tem5;
978 if ( fabs( tem5 ) >= 0.95 )
979 tem5 = tem5 > 0.0 ? 0.95 : -0.95;
985 ecose = axnl * coseo1 + aynl * sineo1;
986 esine = axnl * sineo1 - aynl * coseo1;
987 el2 = axnl * axnl + aynl * aynl;
988 pl = am * ( 1.0 - el2 );
991 kDebug() <<
"Semi-latus rectum < 0.0";
995 rl = am * ( 1.0 - ecose );
996 rdotl = sqrt( am ) * esine / rl;
997 rvdotl = sqrt( pl ) / rl;
998 betal = sqrt( 1.0 - el2 );
999 temp = esine / ( 1.0 + betal );
1000 sinu = am / rl * ( sineo1 - aynl - axnl * temp );
1001 cosu = am / rl * ( coseo1 - axnl + aynl * temp );
1002 su = atan2( sinu, cosu );
1003 sin2u = ( cosu + cosu ) * sinu;
1004 cos2u = 1.0 - 2.0 * sinu * sinu;
1006 temp1 = 0.5 *
J2 * temp;
1007 temp2 = temp1 * temp;
1010 if ( method ==
'd' ) {
1011 cosisq = cosip * cosip;
1012 con41 = 3.0 * cosisq - 1.0;
1013 x1mth2 = 1.0 - cosisq;
1014 x7thm1 = 7.0 * cosisq - 1.0;
1016 mrt = rl * ( 1.0 - 1.5 * temp2 * betal * con41 ) +
1017 0.5 * temp1 * x1mth2 * cos2u;
1018 su = su - 0.25 * temp2 * x7thm1 * sin2u;
1019 xnode = nodep + 1.5 * temp2 * cosip * sin2u;
1020 xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u;
1021 mvt = rdotl - nm * temp1 * x1mth2 * sin2u /
XKE;
1022 rvdot = rvdotl + nm * temp1 * ( x1mth2 * cos2u +
1023 1.5 * con41 ) /
XKE;
1028 snod = sin( xnode );
1029 cnod = cos( xnode );
1034 ux = xmx * sinsu + cnod * cossu;
1035 uy = xmy * sinsu + snod * cossu;
1037 vx = xmx * cossu - cnod * sinsu;
1038 vy = xmy * cossu - snod * sinsu;
1045 sat_posw = sqrt( sat_posx*sat_posx + sat_posy*sat_posy + sat_posz*sat_posz );
1046 sat_velx = ( mvt * ux + rvdot * vx ) * vkmpersec;
1047 sat_vely = ( mvt * uy + rvdot * vy ) * vkmpersec;
1048 sat_velz = ( mvt * uz + rvdot * vz ) * vkmpersec;
1049 m_velocity = sqrt( sat_velx*sat_velx + sat_vely*sat_vely + sat_velz*sat_velz );
1060 kDebug() <<
"Satellite has decayed";
1065 sinlat = sin( data->
geo()->
lat()->radians() );
1066 coslat = cos( data->
geo()->
lat()->radians() );
1067 thetageo = data->
geo()->
LMST( jul_utc );
1068 sintheta = sin( thetageo );
1069 costheta = cos( thetageo );
1070 c = 1.0 / sqrt( 1.0 +
F * (
F - 2.0 ) * sinlat * sinlat );
1071 sq = ( 1.0 -
F ) * ( 1.0 -
F ) * c;
1073 obs_posx = achcp * costheta;
1074 obs_posy = achcp * sintheta;
1076 obs_posw = sqrt( obs_posx*obs_posx + obs_posy*sat_posy + obs_posz*obs_posz );
1077 obs_velx = -
MFACTOR * obs_posy;
1078 obs_vely =
MFACTOR * obs_posx;
1081 m_altitude = sat_posw - obs_posw +
MEANALT;
1084 double range_posx = sat_posx - obs_posx;
1085 double range_posy = sat_posy - obs_posy;
1086 double range_posz = sat_posz - obs_posz;
1087 m_range = sqrt( range_posx*range_posx + range_posy*range_posy + range_posz*range_posz );
1092 double top_s = sinlat*costheta*range_posx + sinlat*sintheta*range_posy - coslat*range_posz;
1093 double top_e = -sintheta*range_posx + costheta*range_posy;
1094 double top_z = coslat*costheta*range_posx + coslat*sintheta*range_posy + sinlat*range_posz;
1096 double azimut = atan( -top_e / top_s );
1101 double elevation = arcSin( top_z / m_range );
1112 double mjd, year, T, M, L, e, C, O, Lsa, nu, R, eps;
1114 mjd = jul_utc - 2415020.0;
1115 year = 1900.0 + mjd / 365.25;
1116 T = ( mjd + deltaET( year ) / (
MINPD * 60.0 ) ) / 36525.0;
1117 M =
DEG2RAD * ( Modulus( 358.47583 + Modulus( 35999.04975 * T, 360.0 ) - ( 0.000150 + 0.0000033 * T ) * T*T, 360.0 ) );
1118 L =
DEG2RAD * ( Modulus( 279.69668 + Modulus( 36000.76892 * T, 360.0 ) + 0.0003025 * T*T, 360.0 ) );
1119 e = 0.01675104 - ( 0.0000418 + 0.000000126 * T ) * T;
1120 C =
DEG2RAD * ( ( 1.919460 - ( 0.004789 + 0.000014 * T ) * T ) *
1121 sin( M ) + ( 0.020094 - 0.000100 * T) *
1122 sin( 2 * M ) + 0.000293 * sin( 3 * M ) );
1123 O =
DEG2RAD * ( Modulus( 259.18 - 1934.142 * T, 360.0 ) );
1124 Lsa = Modulus( L + C -
DEG2RAD * ( 0.00569 -0.00479 * sin( O ) ), TWOPI );
1125 nu = Modulus( M + C, TWOPI);
1126 R = 1.0000002 * ( 1.0 - e*e ) / ( 1.0 + e * cos( nu ) );
1127 eps =
DEG2RAD * ( 23.452294 - ( 0.0130125 + ( 0.00000164 - 0.000000503 * T ) * T ) * T + 0.00256 * cos( O ) );
1130 double sun_posx = R * cos( Lsa );
1131 double sun_posy = R * sin( Lsa ) * cos( eps );
1132 double sun_posz = R * sin( Lsa ) * sin( eps );
1133 double sun_posw = R;
1136 double sd_sun, sd_earth, delta, depth;
1140 double rho_x = sun_posx - sat_posx;
1141 double rho_y = sun_posy - sat_posy;
1142 double rho_z = sun_posz - sat_posz;
1143 double rho_w = sqrt( rho_x*rho_x + rho_y*rho_y + rho_z*rho_z );
1144 sd_sun = arcSin(
SR / rho_w );
1145 double earth_x = -1.0 * sat_posx;
1146 double earth_y = -1.0 * sat_posy;
1147 double earth_z = -1.0 * sat_posz;
1148 double earth_w = sat_posw;
1149 delta =
PIO2 - arcSin( ( sun_posx*earth_x + sun_posy*earth_y + sun_posz*earth_z ) / ( sun_posw*earth_w ) );
1150 depth = sd_earth - sd_sun - delta;
1153 m_is_eclipsed = sd_earth >= sd_sun && depth >= 0;
1154 m_is_visible = !m_is_eclipsed && sun->
alt().
Degrees() <= -12.0 && elevation >= 0.0;
1159 double Satellite::arcSin(
double arg )
1161 if ( fabs( arg ) >= 1. )
1164 else if ( arg < 0. )
1169 return( atan( arg / sqrt( 1. - arg*arg ) ) );
1172 double Satellite::deltaET(
double year )
1176 delta_et=26.465+0.747622*(year-1950)+1.886913*sin(TWOPI*(year-1975)/33);
1181 double Satellite::Modulus(
double arg1,
double arg2)
1198 return m_is_visible;
1203 return m_is_selected;
1211 void Satellite::initPopupMenu(
KSPopupMenu *pmenu ) {
static QStringList selectedSatellites()
Get Selected satellites.
KStarsData is the backbone of KStars.
Child class of KSPlanetBase; encapsulates information about the Sun.
void setAz(dms az)
Sets Az, the Azimuth.
static KStarsData * Instance()
const double & Degrees() const
void setLongName(const QString &longname=QString())
Set the object's long name.
Satellite(QString name, QString line1, QString line2)
Constructor.
double LMST(double jd)
Local Mean Sidereal Time.
const KStarsDateTime & utc() const
void HorizontalToEquatorial(const dms *LST, const dms *lat)
Determine the (RA, Dec) coordinates of the SkyPoint from its (Altitude, Azimuth) coordinates, given the local sidereal time and the observer's latitude.
void setName2(const QString &name2=QString())
Set the object's secondary name.
SkyMapComposite * skyComposite()
void setType(int t)
Set the object's type identifier to the argument.
void setMag(float m)
Set the object's sorting magnitude.
virtual SkyObject * findByName(const QString &name)
Search the children of this SkyMapComposite for a SkyObject whose name matches the argument...
void setAlt(dms alt)
Sets Alt, the Altitude.
void setName(const QString &name)
Set the object's primary name.
void updatePos()
Update satellite position.
void setSelected(bool selected)
Select or not the satellite.