satellite-sgp4unit.cc
Go to the documentation of this file.
1 /* ----------------------------------------------------------------
2  *
3  * sgp4unit.h
4  *
5  * this file contains the sgp4 procedures for analytical propagation
6  * of a satellite. the code was originally released in the 1980 and 1986
7  * spacetrack papers. a detailed discussion of the theory and history
8  * may be found in the 2006 aiaa paper by vallado, crawford, hujsak,
9  * and kelso.
10  *
11  * companion code for
12  * fundamentals of astrodynamics and applications
13  * 2007
14  * by david vallado
15  *
16  * (w) 719-573-2600, email dvallado@agi.com
17  *
18  * current :
19  * 30 Dec 11 david vallado
20  * consolidate updated code version
21  * 30 Aug 10 david vallado
22  * delete unused variables in initl
23  * replace pow inetger 2, 3 with multiplies for speed
24  * changes :
25  * 3 Nov 08 david vallado
26  * put returns in for error codes
27  * 29 sep 08 david vallado
28  * fix atime for faster operation in dspace
29  * add operationmode for afspc (a) or improved (i)
30  * performance mode
31  * 20 apr 07 david vallado
32  * misc fixes for constants
33  * 11 aug 06 david vallado
34  * chg lyddane choice back to strn3, constants, misc doc
35  * 15 dec 05 david vallado
36  * misc fixes
37  * 26 jul 05 david vallado
38  * fixes for paper
39  * note that each fix is preceded by a
40  * comment with "sgp4fix" and an explanation of
41  * what was changed
42  * 10 aug 04 david vallado
43  * 2nd printing baseline working
44  * 14 may 01 david vallado
45  * 2nd edition baseline
46  * 80 norad
47  * original baseline
48  *
49  * code from https://gitlab.inesctec.pt/pmms/ns3-satellite
50  * ---------------------------------------------------------------- */
51 
52 #include "satellite-sgp4unit.h"
53 
54 /* ----------- local functions - only ever used internally by sgp4 ---------- */
55 static void dpper(double e3,
56  double ee2,
57  double peo,
58  double pgho,
59  double pho,
60  double pinco,
61  double plo,
62  double se2,
63  double se3,
64  double sgh2,
65  double sgh3,
66  double sgh4,
67  double sh2,
68  double sh3,
69  double si2,
70  double si3,
71  double sl2,
72  double sl3,
73  double sl4,
74  double t,
75  double xgh2,
76  double xgh3,
77  double xgh4,
78  double xh2,
79  double xh3,
80  double xi2,
81  double xi3,
82  double xl2,
83  double xl3,
84  double xl4,
85  double zmol,
86  double zmos,
87  double inclo,
88  char init,
89  double& ep,
90  double& inclp,
91  double& nodep,
92  double& argpp,
93  double& mp,
94  char opsmode);
95 
96 static void dscom(double epoch,
97  double ep,
98  double argpp,
99  double tc,
100  double inclp,
101  double nodep,
102  double np,
103  double& snodm,
104  double& cnodm,
105  double& sinim,
106  double& cosim,
107  double& sinomm,
108  double& cosomm,
109  double& day,
110  double& e3,
111  double& ee2,
112  double& em,
113  double& emsq,
114  double& gam,
115  double& peo,
116  double& pgho,
117  double& pho,
118  double& pinco,
119  double& plo,
120  double& rtemsq,
121  double& se2,
122  double& se3,
123  double& sgh2,
124  double& sgh3,
125  double& sgh4,
126  double& sh2,
127  double& sh3,
128  double& si2,
129  double& si3,
130  double& sl2,
131  double& sl3,
132  double& sl4,
133  double& s1,
134  double& s2,
135  double& s3,
136  double& s4,
137  double& s5,
138  double& s6,
139  double& s7,
140  double& ss1,
141  double& ss2,
142  double& ss3,
143  double& ss4,
144  double& ss5,
145  double& ss6,
146  double& ss7,
147  double& sz1,
148  double& sz2,
149  double& sz3,
150  double& sz11,
151  double& sz12,
152  double& sz13,
153  double& sz21,
154  double& sz22,
155  double& sz23,
156  double& sz31,
157  double& sz32,
158  double& sz33,
159  double& xgh2,
160  double& xgh3,
161  double& xgh4,
162  double& xh2,
163  double& xh3,
164  double& xi2,
165  double& xi3,
166  double& xl2,
167  double& xl3,
168  double& xl4,
169  double& nm,
170  double& z1,
171  double& z2,
172  double& z3,
173  double& z11,
174  double& z12,
175  double& z13,
176  double& z21,
177  double& z22,
178  double& z23,
179  double& z31,
180  double& z32,
181  double& z33,
182  double& zmol,
183  double& zmos);
184 
185 static void dsinit(gravconsttype whichconst,
186  double cosim,
187  double emsq,
188  double argpo,
189  double s1,
190  double s2,
191  double s3,
192  double s4,
193  double s5,
194  double sinim,
195  double ss1,
196  double ss2,
197  double ss3,
198  double ss4,
199  double ss5,
200  double sz1,
201  double sz3,
202  double sz11,
203  double sz13,
204  double sz21,
205  double sz23,
206  double sz31,
207  double sz33,
208  double t,
209  double tc,
210  double gsto,
211  double mo,
212  double mdot,
213  double no,
214  double nodeo,
215  double nodedot,
216  double xpidot,
217  double z1,
218  double z3,
219  double z11,
220  double z13,
221  double z21,
222  double z23,
223  double z31,
224  double z33,
225  double ecco,
226  double eccsq,
227  double& em,
228  double& argpm,
229  double& inclm,
230  double& mm,
231  double& nm,
232  double& nodem,
233  int& irez,
234  double& atime,
235  double& d2201,
236  double& d2211,
237  double& d3210,
238  double& d3222,
239  double& d4410,
240  double& d4422,
241  double& d5220,
242  double& d5232,
243  double& d5421,
244  double& d5433,
245  double& dedt,
246  double& didt,
247  double& dmdt,
248  double& dndt,
249  double& dnodt,
250  double& domdt,
251  double& del1,
252  double& del2,
253  double& del3,
254  double& xfact,
255  double& xlamo,
256  double& xli,
257  double& xni);
258 
259 static void dspace(int irez,
260  double d2201,
261  double d2211,
262  double d3210,
263  double d3222,
264  double d4410,
265  double d4422,
266  double d5220,
267  double d5232,
268  double d5421,
269  double d5433,
270  double dedt,
271  double del1,
272  double del2,
273  double del3,
274  double didt,
275  double dmdt,
276  double dnodt,
277  double domdt,
278  double argpo,
279  double argpdot,
280  double t,
281  double tc,
282  double gsto,
283  double xfact,
284  double xlamo,
285  double no,
286  double& atime,
287  double& em,
288  double& argpm,
289  double& inclm,
290  double& xli,
291  double& mm,
292  double& xni,
293  double& nodem,
294  double& dndt,
295  double& nm);
296 
297 static void initl(long int satn,
298  gravconsttype whichconst,
299  double ecco,
300  double epoch,
301  double inclo,
302  double& no,
303  char& method,
304  double& ainv,
305  double& ao,
306  double& con41,
307  double& con42,
308  double& cosio,
309  double& cosio2,
310  double& eccsq,
311  double& omeosq,
312  double& posq,
313  double& rp,
314  double& rteosq,
315  double& sinio,
316  double& gsto,
317  char opsmode);
318 
319 /* -----------------------------------------------------------------------------
320  *
321  * procedure dpper
322  *
323  * this procedure provides deep space long period periodic contributions
324  * to the mean elements. by design, these periodics are zero at epoch.
325  * this used to be dscom which included initialization, but it's really a
326  * recurring function.
327  *
328  * author : david vallado 719-573-2600 28 jun 2005
329  *
330  * inputs :
331  * e3 -
332  * ee2 -
333  * peo -
334  * pgho -
335  * pho -
336  * pinco -
337  * plo -
338  * se2 , se3 , sgh2, sgh3, sgh4, sh2, sh3, si2, si3, sl2, sl3, sl4 -
339  * t -
340  * xh2, xh3, xi2, xi3, xl2, xl3, xl4 -
341  * zmol -
342  * zmos -
343  * ep - eccentricity 0.0 - 1.0
344  * inclo - inclination - needed for lyddane modification
345  * nodep - right ascension of ascending node
346  * argpp - argument of perigee
347  * mp - mean anomaly
348  *
349  * outputs :
350  * ep - eccentricity 0.0 - 1.0
351  * inclp - inclination
352  * nodep - right ascension of ascending node
353  * argpp - argument of perigee
354  * mp - mean anomaly
355  *
356  * locals :
357  * alfdp -
358  * betdp -
359  * cosip , sinip , cosop , sinop ,
360  * dalf -
361  * dbet -
362  * dls -
363  * f2, f3 -
364  * pe -
365  * pgh -
366  * ph -
367  * pinc -
368  * pl -
369  * sel , ses , sghl , sghs , shl , shs , sil , sinzf , sis ,
370  * sll , sls
371  * xls -
372  * xnoh -
373  * zf -
374  * zm -
375  *
376  * coupling :
377  * none.
378  *
379  * references :
380  * hoots, roehrich, norad spacetrack report #3 1980
381  * hoots, norad spacetrack report #6 1986
382  * hoots, schumacher and glover 2004
383  * vallado, crawford, hujsak, kelso 2006
384  * ----------------------------------------------------------------------------*/
385 
386 static void
387 dpper(double e3,
388  double ee2,
389  double peo,
390  double pgho,
391  double pho,
392  double pinco,
393  double plo,
394  double se2,
395  double se3,
396  double sgh2,
397  double sgh3,
398  double sgh4,
399  double sh2,
400  double sh3,
401  double si2,
402  double si3,
403  double sl2,
404  double sl3,
405  double sl4,
406  double t,
407  double xgh2,
408  double xgh3,
409  double xgh4,
410  double xh2,
411  double xh3,
412  double xi2,
413  double xi3,
414  double xl2,
415  double xl3,
416  double xl4,
417  double zmol,
418  double zmos,
419  double inclo,
420  char init,
421  double& ep,
422  double& inclp,
423  double& nodep,
424  double& argpp,
425  double& mp,
426  char opsmode)
427 {
428  /* --------------------- local variables ------------------------ */
429  const double twopi = 2.0 * pi;
430  double alfdp, betdp, cosip, cosop, dalf, dbet, dls, f2, f3, pe, pgh, ph, pinc, pl, sel, ses,
431  sghl, sghs, shll, shs, sil, sinip, sinop, sinzf, sis, sll, sls, xls, xnoh, zf, zm, zel, zes,
432  znl, zns;
433 
434  /* ---------------------- constants ----------------------------- */
435  zns = 1.19459e-5;
436  zes = 0.01675;
437  znl = 1.5835218e-4;
438  zel = 0.05490;
439 
440  /* --------------- calculate time varying periodics ----------- */
441  zm = zmos + zns * t;
442  // be sure that the initial call has time set to zero
443  if (init == 'y')
444  zm = zmos;
445  zf = zm + 2.0 * zes * sin(zm);
446  sinzf = sin(zf);
447  f2 = 0.5 * sinzf * sinzf - 0.25;
448  f3 = -0.5 * sinzf * cos(zf);
449  ses = se2 * f2 + se3 * f3;
450  sis = si2 * f2 + si3 * f3;
451  sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
452  sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
453  shs = sh2 * f2 + sh3 * f3;
454  zm = zmol + znl * t;
455  if (init == 'y')
456  zm = zmol;
457  zf = zm + 2.0 * zel * sin(zm);
458  sinzf = sin(zf);
459  f2 = 0.5 * sinzf * sinzf - 0.25;
460  f3 = -0.5 * sinzf * cos(zf);
461  sel = ee2 * f2 + e3 * f3;
462  sil = xi2 * f2 + xi3 * f3;
463  sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf;
464  sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
465  shll = xh2 * f2 + xh3 * f3;
466  pe = ses + sel;
467  pinc = sis + sil;
468  pl = sls + sll;
469  pgh = sghs + sghl;
470  ph = shs + shll;
471 
472  if (init == 'n')
473  {
474  pe = pe - peo;
475  pinc = pinc - pinco;
476  pl = pl - plo;
477  pgh = pgh - pgho;
478  ph = ph - pho;
479  inclp = inclp + pinc;
480  ep = ep + pe;
481  sinip = sin(inclp);
482  cosip = cos(inclp);
483 
484  /* ----------------- apply periodics directly ------------ */
485  // sgp4fix for lyddane choice
486  // strn3 used original inclination - this is technically feasible
487  // gsfc used perturbed inclination - also technically feasible
488  // probably best to readjust the 0.2 limit value and limit discontinuity
489  // 0.2 rad = 11.45916 deg
490  // use next line for original strn3 approach and original inclination
491  // if (inclo >= 0.2)
492  // use next line for gsfc version and perturbed inclination
493  if (inclp >= 0.2)
494  {
495  ph = ph / sinip;
496  pgh = pgh - cosip * ph;
497  argpp = argpp + pgh;
498  nodep = nodep + ph;
499  mp = mp + pl;
500  }
501  else
502  {
503  /* ---- apply periodics with lyddane modification ---- */
504  sinop = sin(nodep);
505  cosop = cos(nodep);
506  alfdp = sinip * sinop;
507  betdp = sinip * cosop;
508  dalf = ph * cosop + pinc * cosip * sinop;
509  dbet = -ph * sinop + pinc * cosip * cosop;
510  alfdp = alfdp + dalf;
511  betdp = betdp + dbet;
512  nodep = fmod(nodep, twopi);
513  // sgp4fix for afspc written intrinsic functions
514  // nodep used without a trigonometric function ahead
515  if ((nodep < 0.0) && (opsmode == 'a'))
516  nodep = nodep + twopi;
517  xls = mp + argpp + cosip * nodep;
518  dls = pl + pgh - pinc * nodep * sinip;
519  xls = xls + dls;
520  xnoh = nodep;
521  nodep = atan2(alfdp, betdp);
522  // sgp4fix for afspc written intrinsic functions
523  // nodep used without a trigonometric function ahead
524  if ((nodep < 0.0) && (opsmode == 'a'))
525  nodep = nodep + twopi;
526  if (fabs(xnoh - nodep) > pi)
527  {
528  if (nodep < xnoh)
529  nodep = nodep + twopi;
530  else
531  nodep = nodep - twopi;
532  }
533  mp = mp + pl;
534  argpp = xls - mp - cosip * nodep;
535  }
536  } // if init == 'n'
537 } // end dpper
538 
539 /*-----------------------------------------------------------------------------
540  *
541  * procedure dscom
542  *
543  * this procedure provides deep space common items used by both the secular
544  * and periodics subroutines. input is provided as shown. this routine
545  * used to be called dpper, but the functions inside weren't well organized.
546  *
547  * author : david vallado 719-573-2600 28 jun 2005
548  *
549  * inputs :
550  * epoch -
551  * ep - eccentricity
552  * argpp - argument of perigee
553  * tc -
554  * inclp - inclination
555  * nodep - right ascension of ascending node
556  * np - mean motion
557  *
558  * outputs :
559  * sinim , cosim , sinomm , cosomm , snodm , cnodm
560  * day -
561  * e3 -
562  * ee2 -
563  * em - eccentricity
564  * emsq - eccentricity squared
565  * gam -
566  * peo -
567  * pgho -
568  * pho -
569  * pinco -
570  * plo -
571  * rtemsq -
572  * se2, se3 -
573  * sgh2, sgh3, sgh4 -
574  * sh2, sh3, si2, si3, sl2, sl3, sl4 -
575  * s1, s2, s3, s4, s5, s6, s7 -
576  * ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3 -
577  * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 -
578  * xgh2, xgh3, xgh4, xh2, xh3, xi2, xi3, xl2, xl3, xl4 -
579  * nm - mean motion
580  * z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33 -
581  * zmol -
582  * zmos -
583  *
584  * locals :
585  * a1, a2, a3, a4, a5, a6, a7, a8, a9, a10 -
586  * betasq -
587  * cc -
588  * ctem, stem -
589  * x1, x2, x3, x4, x5, x6, x7, x8 -
590  * xnodce -
591  * xnoi -
592  * zcosg , zsing , zcosgl , zsingl , zcosh , zsinh , zcoshl , zsinhl ,
593  * zcosi , zsini , zcosil , zsinil ,
594  * zx -
595  * zy -
596  *
597  * coupling :
598  * none.
599  *
600  * references :
601  * hoots, roehrich, norad spacetrack report #3 1980
602  * hoots, norad spacetrack report #6 1986
603  * hoots, schumacher and glover 2004
604  * vallado, crawford, hujsak, kelso 2006
605  * ----------------------------------------------------------------------------*/
606 
607 static void
608 dscom(double epoch,
609  double ep,
610  double argpp,
611  double tc,
612  double inclp,
613  double nodep,
614  double np,
615  double& snodm,
616  double& cnodm,
617  double& sinim,
618  double& cosim,
619  double& sinomm,
620  double& cosomm,
621  double& day,
622  double& e3,
623  double& ee2,
624  double& em,
625  double& emsq,
626  double& gam,
627  double& peo,
628  double& pgho,
629  double& pho,
630  double& pinco,
631  double& plo,
632  double& rtemsq,
633  double& se2,
634  double& se3,
635  double& sgh2,
636  double& sgh3,
637  double& sgh4,
638  double& sh2,
639  double& sh3,
640  double& si2,
641  double& si3,
642  double& sl2,
643  double& sl3,
644  double& sl4,
645  double& s1,
646  double& s2,
647  double& s3,
648  double& s4,
649  double& s5,
650  double& s6,
651  double& s7,
652  double& ss1,
653  double& ss2,
654  double& ss3,
655  double& ss4,
656  double& ss5,
657  double& ss6,
658  double& ss7,
659  double& sz1,
660  double& sz2,
661  double& sz3,
662  double& sz11,
663  double& sz12,
664  double& sz13,
665  double& sz21,
666  double& sz22,
667  double& sz23,
668  double& sz31,
669  double& sz32,
670  double& sz33,
671  double& xgh2,
672  double& xgh3,
673  double& xgh4,
674  double& xh2,
675  double& xh3,
676  double& xi2,
677  double& xi3,
678  double& xl2,
679  double& xl3,
680  double& xl4,
681  double& nm,
682  double& z1,
683  double& z2,
684  double& z3,
685  double& z11,
686  double& z12,
687  double& z13,
688  double& z21,
689  double& z22,
690  double& z23,
691  double& z31,
692  double& z32,
693  double& z33,
694  double& zmol,
695  double& zmos)
696 {
697  /* -------------------------- constants ------------------------- */
698  const double zes = 0.01675;
699  const double zel = 0.05490;
700  const double c1ss = 2.9864797e-6;
701  const double c1l = 4.7968065e-7;
702  const double zsinis = 0.39785416;
703  const double zcosis = 0.91744867;
704  const double zcosgs = 0.1945905;
705  const double zsings = -0.98088458;
706  const double twopi = 2.0 * pi;
707 
708  /* --------------------- local variables ------------------------ */
709  int lsflg;
710  double a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, betasq, cc, ctem, stem, x1, x2, x3, x4, x5, x6,
711  x7, x8, xnodce, xnoi, zcosg, zcosgl, zcosh, zcoshl, zcosi, zcosil, zsing, zsingl, zsinh,
712  zsinhl, zsini, zsinil, zx, zy;
713 
714  nm = np;
715  em = ep;
716  snodm = sin(nodep);
717  cnodm = cos(nodep);
718  sinomm = sin(argpp);
719  cosomm = cos(argpp);
720  sinim = sin(inclp);
721  cosim = cos(inclp);
722  emsq = em * em;
723  betasq = 1.0 - emsq;
724  rtemsq = sqrt(betasq);
725 
726  /* ----------------- initialize lunar solar terms --------------- */
727  peo = 0.0;
728  pinco = 0.0;
729  plo = 0.0;
730  pgho = 0.0;
731  pho = 0.0;
732  day = epoch + 18261.5 + tc / 1440.0;
733  xnodce = fmod(4.5236020 - 9.2422029e-4 * day, twopi);
734  stem = sin(xnodce);
735  ctem = cos(xnodce);
736  zcosil = 0.91375164 - 0.03568096 * ctem;
737  zsinil = sqrt(1.0 - zcosil * zcosil);
738  zsinhl = 0.089683511 * stem / zsinil;
739  zcoshl = sqrt(1.0 - zsinhl * zsinhl);
740  gam = 5.8351514 + 0.0019443680 * day;
741  zx = 0.39785416 * stem / zsinil;
742  zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
743  zx = atan2(zx, zy);
744  zx = gam + zx - xnodce;
745  zcosgl = cos(zx);
746  zsingl = sin(zx);
747 
748  /* ------------------------- do solar terms --------------------- */
749  zcosg = zcosgs;
750  zsing = zsings;
751  zcosi = zcosis;
752  zsini = zsinis;
753  zcosh = cnodm;
754  zsinh = snodm;
755  cc = c1ss;
756  xnoi = 1.0 / nm;
757 
758  for (lsflg = 1; lsflg <= 2; lsflg++)
759  {
760  a1 = zcosg * zcosh + zsing * zcosi * zsinh;
761  a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
762  a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
763  a8 = zsing * zsini;
764  a9 = zsing * zsinh + zcosg * zcosi * zcosh;
765  a10 = zcosg * zsini;
766  a2 = cosim * a7 + sinim * a8;
767  a4 = cosim * a9 + sinim * a10;
768  a5 = -sinim * a7 + cosim * a8;
769  a6 = -sinim * a9 + cosim * a10;
770 
771  x1 = a1 * cosomm + a2 * sinomm;
772  x2 = a3 * cosomm + a4 * sinomm;
773  x3 = -a1 * sinomm + a2 * cosomm;
774  x4 = -a3 * sinomm + a4 * cosomm;
775  x5 = a5 * sinomm;
776  x6 = a6 * sinomm;
777  x7 = a5 * cosomm;
778  x8 = a6 * cosomm;
779 
780  z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3;
781  z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4;
782  z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4;
783  z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq;
784  z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq;
785  z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq;
786  z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5);
787  z12 = -6.0 * (a1 * a6 + a3 * a5) +
788  emsq * (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5));
789  z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6);
790  z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7);
791  z22 = 6.0 * (a4 * a5 + a2 * a6) +
792  emsq * (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8));
793  z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8);
794  z1 = z1 + z1 + betasq * z31;
795  z2 = z2 + z2 + betasq * z32;
796  z3 = z3 + z3 + betasq * z33;
797  s3 = cc * xnoi;
798  s2 = -0.5 * s3 / rtemsq;
799  s4 = s3 * rtemsq;
800  s1 = -15.0 * em * s4;
801  s5 = x1 * x3 + x2 * x4;
802  s6 = x2 * x3 + x1 * x4;
803  s7 = x2 * x4 - x1 * x3;
804 
805  /* ----------------------- do lunar terms ------------------- */
806  if (lsflg == 1)
807  {
808  ss1 = s1;
809  ss2 = s2;
810  ss3 = s3;
811  ss4 = s4;
812  ss5 = s5;
813  ss6 = s6;
814  ss7 = s7;
815  sz1 = z1;
816  sz2 = z2;
817  sz3 = z3;
818  sz11 = z11;
819  sz12 = z12;
820  sz13 = z13;
821  sz21 = z21;
822  sz22 = z22;
823  sz23 = z23;
824  sz31 = z31;
825  sz32 = z32;
826  sz33 = z33;
827  zcosg = zcosgl;
828  zsing = zsingl;
829  zcosi = zcosil;
830  zsini = zsinil;
831  zcosh = zcoshl * cnodm + zsinhl * snodm;
832  zsinh = snodm * zcoshl - cnodm * zsinhl;
833  cc = c1l;
834  }
835  }
836 
837  zmol = fmod(4.7199672 + 0.22997150 * day - gam, twopi);
838  zmos = fmod(6.2565837 + 0.017201977 * day, twopi);
839 
840  /* ------------------------ do solar terms ---------------------- */
841  se2 = 2.0 * ss1 * ss6;
842  se3 = 2.0 * ss1 * ss7;
843  si2 = 2.0 * ss2 * sz12;
844  si3 = 2.0 * ss2 * (sz13 - sz11);
845  sl2 = -2.0 * ss3 * sz2;
846  sl3 = -2.0 * ss3 * (sz3 - sz1);
847  sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes;
848  sgh2 = 2.0 * ss4 * sz32;
849  sgh3 = 2.0 * ss4 * (sz33 - sz31);
850  sgh4 = -18.0 * ss4 * zes;
851  sh2 = -2.0 * ss2 * sz22;
852  sh3 = -2.0 * ss2 * (sz23 - sz21);
853 
854  /* ------------------------ do lunar terms ---------------------- */
855  ee2 = 2.0 * s1 * s6;
856  e3 = 2.0 * s1 * s7;
857  xi2 = 2.0 * s2 * z12;
858  xi3 = 2.0 * s2 * (z13 - z11);
859  xl2 = -2.0 * s3 * z2;
860  xl3 = -2.0 * s3 * (z3 - z1);
861  xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel;
862  xgh2 = 2.0 * s4 * z32;
863  xgh3 = 2.0 * s4 * (z33 - z31);
864  xgh4 = -18.0 * s4 * zel;
865  xh2 = -2.0 * s2 * z22;
866  xh3 = -2.0 * s2 * (z23 - z21);
867 } // end dscom
868 
869 /*-----------------------------------------------------------------------------
870  *
871  * procedure dsinit
872  *
873  * this procedure provides deep space contributions to mean motion dot due
874  * to geopotential resonance with half day and one day orbits.
875  *
876  * author : david vallado 719-573-2600 28 jun 2005
877  *
878  * inputs :
879  * cosim, sinim-
880  * emsq - eccentricity squared
881  * argpo - argument of perigee
882  * s1, s2, s3, s4, s5 -
883  * ss1, ss2, ss3, ss4, ss5 -
884  * sz1, sz3, sz11, sz13, sz21, sz23, sz31, sz33 -
885  * t - time
886  * tc -
887  * gsto - greenwich sidereal time rad
888  * mo - mean anomaly
889  * mdot - mean anomaly dot (rate)
890  * no - mean motion
891  * nodeo - right ascension of ascending node
892  * nodedot - right ascension of ascending node dot (rate)
893  * xpidot -
894  * z1, z3, z11, z13, z21, z23, z31, z33 -
895  * eccm - eccentricity
896  * argpm - argument of perigee
897  * inclm - inclination
898  * mm - mean anomaly
899  * xn - mean motion
900  * nodem - right ascension of ascending node
901  *
902  * outputs :
903  * em - eccentricity
904  * argpm - argument of perigee
905  * inclm - inclination
906  * mm - mean anomaly
907  * nm - mean motion
908  * nodem - right ascension of ascending node
909  * irez - flag for resonance 0-none, 1-one day, 2-half day
910  * atime -
911  * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 -
912  * dedt -
913  * didt -
914  * dmdt -
915  * dndt -
916  * dnodt -
917  * domdt -
918  * del1, del2, del3 -
919  * ses , sghl , sghs , sgs , shl , shs , sis , sls
920  * theta -
921  * xfact -
922  * xlamo -
923  * xli -
924  * xni
925  *
926  * locals :
927  * ainv2 -
928  * aonv -
929  * cosisq -
930  * eoc -
931  * f220, f221, f311, f321, f322, f330, f441, f442, f522, f523, f542, f543 -
932  * g200, g201, g211, g300, g310, g322, g410, g422, g520, g521, g532, g533 -
933  * sini2 -
934  * temp -
935  * temp1 -
936  * theta -
937  * xno2 -
938  *
939  * coupling :
940  * getgravconst
941  *
942  * references :
943  * hoots, roehrich, norad spacetrack report #3 1980
944  * hoots, norad spacetrack report #6 1986
945  * hoots, schumacher and glover 2004
946  * vallado, crawford, hujsak, kelso 2006
947  * ----------------------------------------------------------------------------*/
948 
949 static void
950 dsinit(gravconsttype whichconst,
951  double cosim,
952  double emsq,
953  double argpo,
954  double s1,
955  double s2,
956  double s3,
957  double s4,
958  double s5,
959  double sinim,
960  double ss1,
961  double ss2,
962  double ss3,
963  double ss4,
964  double ss5,
965  double sz1,
966  double sz3,
967  double sz11,
968  double sz13,
969  double sz21,
970  double sz23,
971  double sz31,
972  double sz33,
973  double t,
974  double tc,
975  double gsto,
976  double mo,
977  double mdot,
978  double no,
979  double nodeo,
980  double nodedot,
981  double xpidot,
982  double z1,
983  double z3,
984  double z11,
985  double z13,
986  double z21,
987  double z23,
988  double z31,
989  double z33,
990  double ecco,
991  double eccsq,
992  double& em,
993  double& argpm,
994  double& inclm,
995  double& mm,
996  double& nm,
997  double& nodem,
998  int& irez,
999  double& atime,
1000  double& d2201,
1001  double& d2211,
1002  double& d3210,
1003  double& d3222,
1004  double& d4410,
1005  double& d4422,
1006  double& d5220,
1007  double& d5232,
1008  double& d5421,
1009  double& d5433,
1010  double& dedt,
1011  double& didt,
1012  double& dmdt,
1013  double& dndt,
1014  double& dnodt,
1015  double& domdt,
1016  double& del1,
1017  double& del2,
1018  double& del3,
1019  double& xfact,
1020  double& xlamo,
1021  double& xli,
1022  double& xni)
1023 {
1024  /* --------------------- local variables ------------------------ */
1025  const double twopi = 2.0 * pi;
1026 
1027  double ainv2, aonv = 0.0, cosisq, eoc, f220, f221, f311, f321, f322, f330, f441, f442, f522,
1028  f523, f542, f543, g200, g201, g211, g300, g310, g322, g410, g422, g520, g521,
1029  g532, g533, ses, sgs, sghl, sghs, shs, shll, sis, sini2, sls, temp, temp1, theta,
1030  xno2, q22, q31, q33, root22, root44, root54, rptim, root32, root52, x2o3, xke,
1031  znl, emo, zns, emsqo, tumin, mu, radiusearthkm, j2, j3, j4, j3oj2;
1032 
1033  // ensure that all variables are initialized to make the compiler happy
1034  ainv2 = cosisq = eoc = f220 = f221 = f311 = f321 = f322 = f330 = f441 = f442 = f522 = f523 =
1035  f542 = f543 = g200 = g201 = g211 = g300 = g310 = g322 = g410 = g422 = g520 = g521 = g532 =
1036  g533 = ses = sgs = sghl = sghs = shs = shll = sis = sini2 = sls = temp = temp1 = theta =
1037  xno2 = q22 = q31 = q33 = root22 = root44 = root54 = rptim = root32 = root52 = x2o3 =
1038  xke = znl = emo = zns = emsqo = tumin = mu = radiusearthkm = j2 = j3 = j4 =
1039  j3oj2 = 0.0;
1040 
1041  q22 = 1.7891679e-6;
1042  q31 = 2.1460748e-6;
1043  q33 = 2.2123015e-7;
1044  root22 = 1.7891679e-6;
1045  root44 = 7.3636953e-9;
1046  root54 = 2.1765803e-9;
1047  rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec
1048  root32 = 3.7393792e-7;
1049  root52 = 1.1428639e-7;
1050  x2o3 = 2.0 / 3.0;
1051  znl = 1.5835218e-4;
1052  zns = 1.19459e-5;
1053 
1054  // sgp4fix identify constants and allow alternate values
1055  getgravconst(whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2);
1056 
1057  /* -------------------- deep space initialization ------------ */
1058  irez = 0;
1059  if ((nm < 0.0052359877) && (nm > 0.0034906585))
1060  irez = 1;
1061  if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5))
1062  irez = 2;
1063 
1064  /* ------------------------ do solar terms ------------------- */
1065  ses = ss1 * zns * ss5;
1066  sis = ss2 * zns * (sz11 + sz13);
1067  sls = -zns * ss3 * (sz1 + sz3 - 14.0 - 6.0 * emsq);
1068  sghs = ss4 * zns * (sz31 + sz33 - 6.0);
1069  shs = -zns * ss2 * (sz21 + sz23);
1070  // sgp4fix for 180 deg incl
1071  if ((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2))
1072  shs = 0.0;
1073  if (sinim != 0.0)
1074  shs = shs / sinim;
1075  sgs = sghs - cosim * shs;
1076 
1077  /* ------------------------- do lunar terms ------------------ */
1078  dedt = ses + s1 * znl * s5;
1079  didt = sis + s2 * znl * (z11 + z13);
1080  dmdt = sls - znl * s3 * (z1 + z3 - 14.0 - 6.0 * emsq);
1081  sghl = s4 * znl * (z31 + z33 - 6.0);
1082  shll = -znl * s2 * (z21 + z23);
1083  // sgp4fix for 180 deg incl
1084  if ((inclm < 5.2359877e-2) || (inclm > pi - 5.2359877e-2))
1085  shll = 0.0;
1086  domdt = sgs + sghl;
1087  dnodt = shs;
1088  if (sinim != 0.0)
1089  {
1090  domdt = domdt - cosim / sinim * shll;
1091  dnodt = dnodt + shll / sinim;
1092  }
1093 
1094  /* ----------- calculate deep space resonance effects -------- */
1095  dndt = 0.0;
1096  theta = fmod(gsto + tc * rptim, twopi);
1097  em = em + dedt * t;
1098  inclm = inclm + didt * t;
1099  argpm = argpm + domdt * t;
1100  nodem = nodem + dnodt * t;
1101  mm = mm + dmdt * t;
1102  // sgp4fix for negative inclinations
1103  // the following if statement should be commented out
1104  // if (inclm < 0.0)
1105  // {
1106  // inclm = -inclm;
1107  // argpm = argpm - pi;
1108  // nodem = nodem + pi;
1109  // }
1110 
1111  /* -------------- initialize the resonance terms ------------- */
1112  if (irez != 0)
1113  {
1114  aonv = pow(nm / xke, x2o3);
1115 
1116  /* ---------- geopotential resonance for 12 hour orbits ------ */
1117  if (irez == 2)
1118  {
1119  cosisq = cosim * cosim;
1120  emo = em;
1121  em = ecco;
1122  emsqo = emsq;
1123  emsq = eccsq;
1124  eoc = em * emsq;
1125  g201 = -0.306 - (em - 0.64) * 0.440;
1126 
1127  if (em <= 0.65)
1128  {
1129  g211 = 3.616 - 13.2470 * em + 16.2900 * emsq;
1130  g310 = -19.302 + 117.3900 * em - 228.4190 * emsq + 156.5910 * eoc;
1131  g322 = -18.9068 + 109.7927 * em - 214.6334 * emsq + 146.5816 * eoc;
1132  g410 = -41.122 + 242.6940 * em - 471.0940 * emsq + 313.9530 * eoc;
1133  g422 = -146.407 + 841.8800 * em - 1629.014 * emsq + 1083.4350 * eoc;
1134  g520 = -532.114 + 3017.977 * em - 5740.032 * emsq + 3708.2760 * eoc;
1135  }
1136  else
1137  {
1138  g211 = -72.099 + 331.819 * em - 508.738 * emsq + 266.724 * eoc;
1139  g310 = -346.844 + 1582.851 * em - 2415.925 * emsq + 1246.113 * eoc;
1140  g322 = -342.585 + 1554.908 * em - 2366.899 * emsq + 1215.972 * eoc;
1141  g410 = -1052.797 + 4758.686 * em - 7193.992 * emsq + 3651.957 * eoc;
1142  g422 = -3581.690 + 16178.110 * em - 24462.770 * emsq + 12422.520 * eoc;
1143  if (em > 0.715)
1144  g520 = -5149.66 + 29936.92 * em - 54087.36 * emsq + 31324.56 * eoc;
1145  else
1146  g520 = 1464.74 - 4664.75 * em + 3763.64 * emsq;
1147  }
1148  if (em < 0.7)
1149  {
1150  g533 = -919.22770 + 4988.6100 * em - 9064.7700 * emsq + 5542.21 * eoc;
1151  g521 = -822.71072 + 4568.6173 * em - 8491.4146 * emsq + 5337.524 * eoc;
1152  g532 = -853.66600 + 4690.2500 * em - 8624.7700 * emsq + 5341.4 * eoc;
1153  }
1154  else
1155  {
1156  g533 = -37995.780 + 161616.52 * em - 229838.20 * emsq + 109377.94 * eoc;
1157  g521 = -51752.104 + 218913.95 * em - 309468.16 * emsq + 146349.42 * eoc;
1158  g532 = -40023.880 + 170470.89 * em - 242699.48 * emsq + 115605.82 * eoc;
1159  }
1160 
1161  sini2 = sinim * sinim;
1162  f220 = 0.75 * (1.0 + 2.0 * cosim + cosisq);
1163  f221 = 1.5 * sini2;
1164  f321 = 1.875 * sinim * (1.0 - 2.0 * cosim - 3.0 * cosisq);
1165  f322 = -1.875 * sinim * (1.0 + 2.0 * cosim - 3.0 * cosisq);
1166  f441 = 35.0 * sini2 * f220;
1167  f442 = 39.3750 * sini2 * sini2;
1168  f522 = 9.84375 * sinim *
1169  (sini2 * (1.0 - 2.0 * cosim - 5.0 * cosisq) +
1170  0.33333333 * (-2.0 + 4.0 * cosim + 6.0 * cosisq));
1171  f523 = sinim * (4.92187512 * sini2 * (-2.0 - 4.0 * cosim + 10.0 * cosisq) +
1172  6.56250012 * (1.0 + 2.0 * cosim - 3.0 * cosisq));
1173  f542 = 29.53125 * sinim *
1174  (2.0 - 8.0 * cosim + cosisq * (-12.0 + 8.0 * cosim + 10.0 * cosisq));
1175  f543 = 29.53125 * sinim *
1176  (-2.0 - 8.0 * cosim + cosisq * (12.0 + 8.0 * cosim - 10.0 * cosisq));
1177  xno2 = nm * nm;
1178  ainv2 = aonv * aonv;
1179  temp1 = 3.0 * xno2 * ainv2;
1180  temp = temp1 * root22;
1181  d2201 = temp * f220 * g201;
1182  d2211 = temp * f221 * g211;
1183  temp1 = temp1 * aonv;
1184  temp = temp1 * root32;
1185  d3210 = temp * f321 * g310;
1186  d3222 = temp * f322 * g322;
1187  temp1 = temp1 * aonv;
1188  temp = 2.0 * temp1 * root44;
1189  d4410 = temp * f441 * g410;
1190  d4422 = temp * f442 * g422;
1191  temp1 = temp1 * aonv;
1192  temp = temp1 * root52;
1193  d5220 = temp * f522 * g520;
1194  d5232 = temp * f523 * g532;
1195  temp = 2.0 * temp1 * root54;
1196  d5421 = temp * f542 * g521;
1197  d5433 = temp * f543 * g533;
1198  xlamo = fmod(mo + nodeo + nodeo - theta - theta, twopi);
1199  xfact = mdot + dmdt + 2.0 * (nodedot + dnodt - rptim) - no;
1200  em = emo;
1201  emsq = emsqo;
1202  }
1203 
1204  /* ---------------- synchronous resonance terms -------------- */
1205  if (irez == 1)
1206  {
1207  g200 = 1.0 + emsq * (-2.5 + 0.8125 * emsq);
1208  g310 = 1.0 + 2.0 * emsq;
1209  g300 = 1.0 + emsq * (-6.0 + 6.60937 * emsq);
1210  f220 = 0.75 * (1.0 + cosim) * (1.0 + cosim);
1211  f311 = 0.9375 * sinim * sinim * (1.0 + 3.0 * cosim) - 0.75 * (1.0 + cosim);
1212  f330 = 1.0 + cosim;
1213  f330 = 1.875 * f330 * f330 * f330;
1214  del1 = 3.0 * nm * nm * aonv * aonv;
1215  del2 = 2.0 * del1 * f220 * g200 * q22;
1216  del3 = 3.0 * del1 * f330 * g300 * q33 * aonv;
1217  del1 = del1 * f311 * g310 * q31 * aonv;
1218  xlamo = fmod(mo + nodeo + argpo - theta, twopi);
1219  xfact = mdot + xpidot - rptim + dmdt + domdt + dnodt - no;
1220  }
1221 
1222  /* ------------ for sgp4, initialize the integrator ---------- */
1223  xli = xlamo;
1224  xni = no;
1225  atime = 0.0;
1226  nm = no + dndt;
1227  }
1228 } // end dsinit
1229 
1230 /*-----------------------------------------------------------------------------
1231  *
1232  * procedure dspace
1233  *
1234  * this procedure provides deep space contributions to mean elements for
1235  * perturbing third body. these effects have been averaged over one
1236  * revolution of the sun and moon. for earth resonance effects, the
1237  * effects have been averaged over no revolutions of the satellite.
1238  * (mean motion)
1239  *
1240  * author : david vallado 719-573-2600 28 jun 2005
1241  *
1242  * inputs :
1243  * d2201, d2211, d3210, d3222, d4410, d4422, d5220, d5232, d5421, d5433 -
1244  * dedt -
1245  * del1, del2, del3 -
1246  * didt -
1247  * dmdt -
1248  * dnodt -
1249  * domdt -
1250  * irez - flag for resonance 0-none, 1-one day, 2-half day
1251  * argpo - argument of perigee
1252  * argpdot - argument of perigee dot (rate)
1253  * t - time
1254  * tc -
1255  * gsto - gst
1256  * xfact -
1257  * xlamo -
1258  * no - mean motion
1259  * atime -
1260  * em - eccentricity
1261  * ft -
1262  * argpm - argument of perigee
1263  * inclm - inclination
1264  * xli -
1265  * mm - mean anomaly
1266  * xni - mean motion
1267  * nodem - right ascension of ascending node
1268  *
1269  * outputs :
1270  * atime -
1271  * em - eccentricity
1272  * argpm - argument of perigee
1273  * inclm - inclination
1274  * xli -
1275  * mm - mean anomaly
1276  * xni -
1277  * nodem - right ascension of ascending node
1278  * dndt -
1279  * nm - mean motion
1280  *
1281  * locals :
1282  * delt -
1283  * ft -
1284  * theta -
1285  * x2li -
1286  * x2omi -
1287  * xl -
1288  * xldot -
1289  * xnddt -
1290  * xndt -
1291  * xomi -
1292  *
1293  * coupling :
1294  * none -
1295  *
1296  * references :
1297  * hoots, roehrich, norad spacetrack report #3 1980
1298  * hoots, norad spacetrack report #6 1986
1299  * hoots, schumacher and glover 2004
1300  * vallado, crawford, hujsak, kelso 2006
1301  * ----------------------------------------------------------------------------*/
1302 
1303 static void
1304 dspace(int irez,
1305  double d2201,
1306  double d2211,
1307  double d3210,
1308  double d3222,
1309  double d4410,
1310  double d4422,
1311  double d5220,
1312  double d5232,
1313  double d5421,
1314  double d5433,
1315  double dedt,
1316  double del1,
1317  double del2,
1318  double del3,
1319  double didt,
1320  double dmdt,
1321  double dnodt,
1322  double domdt,
1323  double argpo,
1324  double argpdot,
1325  double t,
1326  double tc,
1327  double gsto,
1328  double xfact,
1329  double xlamo,
1330  double no,
1331  double& atime,
1332  double& em,
1333  double& argpm,
1334  double& inclm,
1335  double& xli,
1336  double& mm,
1337  double& xni,
1338  double& nodem,
1339  double& dndt,
1340  double& nm)
1341 {
1342  const double twopi = 2.0 * pi;
1343  int iretn;
1344  double delt, ft, theta, x2li, x2omi, xl, xldot, xnddt, xndt, xomi, g22, g32, g44, g52, g54,
1345  fasx2, fasx4, fasx6, rptim, step2, stepn, stepp;
1346 
1347  fasx2 = 0.13130908;
1348  fasx4 = 2.8843198;
1349  fasx6 = 0.37448087;
1350  g22 = 5.7686396;
1351  g32 = 0.95240898;
1352  g44 = 1.8014998;
1353  g52 = 1.0508330;
1354  g54 = 4.4108898;
1355  rptim = 4.37526908801129966e-3; // this equates to 7.29211514668855e-5 rad/sec
1356  stepp = 720.0;
1357  stepn = -720.0;
1358  step2 = 259200.0;
1359 
1360  /* ----------- calculate deep space resonance effects ----------- */
1361  dndt = 0.0;
1362  theta = fmod(gsto + tc * rptim, twopi);
1363  em = em + dedt * t;
1364 
1365  inclm = inclm + didt * t;
1366  argpm = argpm + domdt * t;
1367  nodem = nodem + dnodt * t;
1368  mm = mm + dmdt * t;
1369 
1370  // sgp4fix for negative inclinations
1371  // the following if statement should be commented out
1372  // if (inclm < 0.0)
1373  // {
1374  // inclm = -inclm;
1375  // argpm = argpm - pi;
1376  // nodem = nodem + pi;
1377  // }
1378 
1379  /* - update resonances : numerical (euler-maclaurin) integration - */
1380  /* ------------------------- epoch restart ---------------------- */
1381  // sgp4fix for propagator problems
1382  // the following integration works for negative time steps and periods
1383  // the specific changes are unknown because the original code was so convoluted
1384 
1385  // sgp4fix take out atime = 0.0 and fix for faster operation
1386  ft = 0.0;
1387  if (irez != 0)
1388  {
1389  // sgp4fix streamline check
1390  if ((atime == 0.0) || (t * atime <= 0.0) || (fabs(t) < fabs(atime)))
1391  {
1392  atime = 0.0;
1393  xni = no;
1394  xli = xlamo;
1395  }
1396  // sgp4fix move check outside loop
1397  if (t > 0.0)
1398  delt = stepp;
1399  else
1400  delt = stepn;
1401 
1402  iretn = 381; // added for do loop
1403  while (iretn == 381)
1404  {
1405  /* ------------------- dot terms calculated ------------- */
1406  /* ----------- near - synchronous resonance terms ------- */
1407  if (irez != 2)
1408  {
1409  xndt = del1 * sin(xli - fasx2) + del2 * sin(2.0 * (xli - fasx4)) +
1410  del3 * sin(3.0 * (xli - fasx6));
1411  xldot = xni + xfact;
1412  xnddt = del1 * cos(xli - fasx2) + 2.0 * del2 * cos(2.0 * (xli - fasx4)) +
1413  3.0 * del3 * cos(3.0 * (xli - fasx6));
1414  xnddt = xnddt * xldot;
1415  }
1416  else
1417  {
1418  /* --------- near - half-day resonance terms -------- */
1419  xomi = argpo + argpdot * atime;
1420  x2omi = xomi + xomi;
1421  x2li = xli + xli;
1422  xndt = d2201 * sin(x2omi + xli - g22) + d2211 * sin(xli - g22) +
1423  d3210 * sin(xomi + xli - g32) + d3222 * sin(-xomi + xli - g32) +
1424  d4410 * sin(x2omi + x2li - g44) + d4422 * sin(x2li - g44) +
1425  d5220 * sin(xomi + xli - g52) + d5232 * sin(-xomi + xli - g52) +
1426  d5421 * sin(xomi + x2li - g54) + d5433 * sin(-xomi + x2li - g54);
1427  xldot = xni + xfact;
1428  xnddt = d2201 * cos(x2omi + xli - g22) + d2211 * cos(xli - g22) +
1429  d3210 * cos(xomi + xli - g32) + d3222 * cos(-xomi + xli - g32) +
1430  d5220 * cos(xomi + xli - g52) + d5232 * cos(-xomi + xli - g52) +
1431  2.0 * (d4410 * cos(x2omi + x2li - g44) + d4422 * cos(x2li - g44) +
1432  d5421 * cos(xomi + x2li - g54) + d5433 * cos(-xomi + x2li - g54));
1433  xnddt = xnddt * xldot;
1434  }
1435 
1436  /* ----------------------- integrator ------------------- */
1437  // sgp4fix move end checks to end of routine
1438  if (fabs(t - atime) >= stepp)
1439  {
1440  iretn = 381;
1441  }
1442  else // exit here
1443  {
1444  ft = t - atime;
1445  iretn = 0;
1446  }
1447 
1448  if (iretn == 381)
1449  {
1450  xli = xli + xldot * delt + xndt * step2;
1451  xni = xni + xndt * delt + xnddt * step2;
1452  atime = atime + delt;
1453  }
1454  } // while iretn = 381
1455 
1456  nm = xni + xndt * ft + xnddt * ft * ft * 0.5;
1457  xl = xli + xldot * ft + xndt * ft * ft * 0.5;
1458  if (irez != 1)
1459  {
1460  mm = xl - 2.0 * nodem + 2.0 * theta;
1461  dndt = nm - no;
1462  }
1463  else
1464  {
1465  mm = xl - nodem - argpm + theta;
1466  dndt = nm - no;
1467  }
1468  nm = no + dndt;
1469  }
1470 } // end dsspace
1471 
1472 /*-----------------------------------------------------------------------------
1473  *
1474  * procedure initl
1475  *
1476  * this procedure initializes the spg4 propagator. all the initialization is
1477  * consolidated here instead of having multiple loops inside other routines.
1478  *
1479  * author : david vallado 719-573-2600 28 jun 2005
1480  *
1481  * inputs :
1482  * ecco - eccentricity 0.0 - 1.0
1483  * epoch - epoch time in days from jan 0, 1950. 0 hr
1484  * inclo - inclination of satellite
1485  * no - mean motion of satellite
1486  * satn - satellite number
1487  *
1488  * outputs :
1489  * ainv - 1.0 / a
1490  * ao - semi major axis
1491  * con41 -
1492  * con42 - 1.0 - 5.0 cos(i)
1493  * cosio - cosine of inclination
1494  * cosio2 - cosio squared
1495  * eccsq - eccentricity squared
1496  * method - flag for deep space 'd', 'n'
1497  * omeosq - 1.0 - ecco * ecco
1498  * posq - semi-parameter squared
1499  * rp - radius of perigee
1500  * rteosq - square root of (1.0 - ecco*ecco)
1501  * sinio - sine of inclination
1502  * gsto - gst at time of observation rad
1503  * no - mean motion of satellite
1504  *
1505  * locals :
1506  * ak -
1507  * d1 -
1508  * del -
1509  * adel -
1510  * po -
1511  *
1512  * coupling :
1513  * getgravconst
1514  * gstime - find greenwich sidereal time from the julian date
1515  *
1516  * references :
1517  * hoots, roehrich, norad spacetrack report #3 1980
1518  * hoots, norad spacetrack report #6 1986
1519  * hoots, schumacher and glover 2004
1520  * vallado, crawford, hujsak, kelso 2006
1521  * ----------------------------------------------------------------------------*/
1522 
1523 static void
1524 initl(long int satn,
1525  gravconsttype whichconst,
1526  double ecco,
1527  double epoch,
1528  double inclo,
1529  double& no,
1530  char& method,
1531  double& ainv,
1532  double& ao,
1533  double& con41,
1534  double& con42,
1535  double& cosio,
1536  double& cosio2,
1537  double& eccsq,
1538  double& omeosq,
1539  double& posq,
1540  double& rp,
1541  double& rteosq,
1542  double& sinio,
1543  double& gsto,
1544  char opsmode)
1545 {
1546  /* --------------------- local variables ------------------------ */
1547  double ak, d1, del, adel, po, x2o3, j2, xke, tumin, mu, radiusearthkm, j3, j4, j3oj2;
1548 
1549  // ensure that all variables are initialized to make the compiler happy
1550  ak = d1 = del = adel = po = x2o3 = j2 = xke = tumin = mu = radiusearthkm = j3 = j4 = j3oj2 =
1551  0.0;
1552 
1553  // sgp4fix use old way of finding gst
1554  double ds70;
1555  double ts70, tfrac, c1, thgr70, fk5r, c1p2p;
1556  const double twopi = 2.0 * pi;
1557 
1558  /* ----------------------- earth constants ---------------------- */
1559  // sgp4fix identify constants and allow alternate values
1560  getgravconst(whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2);
1561  x2o3 = 2.0 / 3.0;
1562 
1563  /* ------------- calculate auxillary epoch quantities ---------- */
1564  eccsq = ecco * ecco;
1565  omeosq = 1.0 - eccsq;
1566  rteosq = sqrt(omeosq);
1567  cosio = cos(inclo);
1568  cosio2 = cosio * cosio;
1569 
1570  /* ------------------ un-kozai the mean motion ----------------- */
1571  ak = pow(xke / no, x2o3);
1572  d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq);
1573  del = d1 / (ak * ak);
1574  adel = ak * (1.0 - del * del - del * (1.0 / 3.0 + 134.0 * del * del / 81.0));
1575  del = d1 / (adel * adel);
1576  no = no / (1.0 + del);
1577 
1578  ao = pow(xke / no, x2o3);
1579  sinio = sin(inclo);
1580  po = ao * omeosq;
1581  con42 = 1.0 - 5.0 * cosio2;
1582  con41 = -con42 - cosio2 - cosio2;
1583  ainv = 1.0 / ao;
1584  posq = po * po;
1585  rp = ao * (1.0 - ecco);
1586  method = 'n';
1587 
1588  // sgp4fix modern approach to finding sidereal time
1589  if (opsmode == 'a')
1590  {
1591  // sgp4fix use old way of finding gst
1592  // count integer number of days from 0 jan 1970
1593  ts70 = epoch - 7305.0;
1594  ds70 = floor(ts70 + 1.0e-8);
1595  tfrac = ts70 - ds70;
1596  // find greenwich location at epoch
1597  c1 = 1.72027916940703639e-2;
1598  thgr70 = 1.7321343856509374;
1599  fk5r = 5.07551419432269442e-15;
1600  c1p2p = c1 + twopi;
1601  gsto = fmod(thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r, twopi);
1602  if (gsto < 0.0)
1603  gsto = gsto + twopi;
1604  }
1605  else
1606  gsto = gstime(epoch + 2433281.5);
1607 
1608 } // end initl
1609 
1610 /*-----------------------------------------------------------------------------
1611  *
1612  * procedure sgp4init
1613  *
1614  * this procedure initializes variables for sgp4.
1615  *
1616  * author : david vallado 719-573-2600 28 jun 2005
1617  *
1618  * inputs :
1619  * opsmode - mode of operation afspc or improved 'a', 'i'
1620  * whichconst - which set of constants to use 72, 84
1621  * satn - satellite number
1622  * bstar - sgp4 type drag coefficient kg/m2er
1623  * ecco - eccentricity
1624  * epoch - epoch time in days from jan 0, 1950. 0 hr
1625  * argpo - argument of perigee (output if ds)
1626  * inclo - inclination
1627  * mo - mean anomaly (output if ds)
1628  * no - mean motion
1629  * nodeo - right ascension of ascending node
1630  *
1631  * outputs :
1632  * satrec - common values for subsequent calls
1633  * return code - non-zero on error.
1634  * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er
1635  * 2 - mean motion less than 0.0
1636  * 3 - pert elements, ecc < 0.0 or ecc > 1.0
1637  * 4 - semi-latus rectum < 0.0
1638  * 5 - epoch elements are sub-orbital
1639  * 6 - satellite has decayed
1640  *
1641  * locals :
1642  * cnodm , snodm , cosim , sinim , cosomm , sinomm
1643  * cc1sq , cc2 , cc3
1644  * coef , coef1
1645  * cosio4 -
1646  * day -
1647  * dndt -
1648  * em - eccentricity
1649  * emsq - eccentricity squared
1650  * eeta -
1651  * etasq -
1652  * gam -
1653  * argpm - argument of perigee
1654  * nodem -
1655  * inclm - inclination
1656  * mm - mean anomaly
1657  * nm - mean motion
1658  * perige - perigee
1659  * pinvsq -
1660  * psisq -
1661  * qzms24 -
1662  * rtemsq -
1663  * s1, s2, s3, s4, s5, s6, s7 -
1664  * sfour -
1665  * ss1, ss2, ss3, ss4, ss5, ss6, ss7 -
1666  * sz1, sz2, sz3
1667  * sz11, sz12, sz13, sz21, sz22, sz23, sz31, sz32, sz33 -
1668  * tc -
1669  * temp -
1670  * temp1, temp2, temp3 -
1671  * tsi -
1672  * xpidot -
1673  * xhdot1 -
1674  * z1, z2, z3 -
1675  * z11, z12, z13, z21, z22, z23, z31, z32, z33 -
1676  *
1677  * coupling :
1678  * getgravconst-
1679  * initl -
1680  * dscom -
1681  * dpper -
1682  * dsinit -
1683  * sgp4 -
1684  *
1685  * references :
1686  * hoots, roehrich, norad spacetrack report #3 1980
1687  * hoots, norad spacetrack report #6 1986
1688  * hoots, schumacher and glover 2004
1689  * vallado, crawford, hujsak, kelso 2006
1690  * ----------------------------------------------------------------------------*/
1691 
1692 bool
1693 sgp4init(gravconsttype whichconst,
1694  char opsmode,
1695  const long int satn,
1696  const double epoch,
1697  const double xbstar,
1698  const double xecco,
1699  const double xargpo,
1700  const double xinclo,
1701  const double xmo,
1702  const double xno,
1703  const double xnodeo,
1704  elsetrec& satrec)
1705 {
1706  /* --------------------- local variables ------------------------ */
1707  double ao, ainv, con42, cosio, sinio, cosio2, eccsq, omeosq, posq, rp, rteosq, cnodm, snodm,
1708  cosim, sinim, cosomm, sinomm, cc1sq, cc2, cc3, coef, coef1, cosio4, day, dndt, em, emsq,
1709  eeta, etasq, gam, argpm, nodem, inclm, mm, nm, perige, pinvsq, psisq, qzms24, rtemsq, s1,
1710  s2, s3, s4, s5, s6, s7, sfour, ss1, ss2, ss3, ss4, ss5, ss6, ss7, sz1, sz2, sz3, sz11, sz12,
1711  sz13, sz21, sz22, sz23, sz31, sz32, sz33, tc, temp, temp1, temp2, temp3, tsi, xpidot,
1712  xhdot1, z1, z2, z3, z11, z12, z13, z21, z22, z23, z31, z32, z33, qzms2t, ss, j2, j3oj2, j4,
1713  x2o3, r[3], v[3], tumin, mu, radiusearthkm, xke, j3, delmotemp, qzms2ttemp, qzms24temp;
1714 
1715  // ensure that all variables are initialized to make the compiler happy
1716  ao = ainv = con42 = cosio = sinio = cosio2 = eccsq = omeosq = posq = rp = rteosq = cnodm =
1717  snodm = cosim = sinim = cosomm = sinomm = cc1sq = cc2 = cc3 = coef = coef1 = cosio4 = day =
1718  dndt = em = emsq = eeta = etasq = gam = argpm = nodem = inclm = mm = nm = perige =
1719  pinvsq = psisq = qzms24 = rtemsq = s1 = s2 = s3 = s4 = s5 = s6 = s7 = sfour = ss1 =
1720  ss2 = ss3 = ss4 = ss5 = ss6 = ss7 = sz1 = sz2 = sz3 = sz11 = sz12 = sz13 =
1721  sz21 = sz22 = sz23 = sz31 = sz32 = sz33 = tc = temp = temp1 = temp2 =
1722  temp3 = tsi = xpidot = xhdot1 = z1 = z2 = z3 = z11 = z12 = z13 = z21 =
1723  z22 = z23 = z31 = z32 = z33 = qzms2t = ss = j2 = j3oj2 = j4 = x2o3 =
1724  tumin = mu = radiusearthkm = xke = j3 = delmotemp = qzms2ttemp =
1725  qzms24temp = 0.0;
1726 
1727  /* ------------------------ initialization --------------------- */
1728  // sgp4fix divisor for divide by zero check on inclination
1729  // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to
1730  // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency
1731  const double temp4 = 1.5e-12;
1732 
1733  /* ----------- set all near earth variables to zero ------------ */
1734  satrec.isimp = 0;
1735  satrec.method = 'n';
1736  satrec.aycof = 0.0;
1737  satrec.con41 = 0.0;
1738  satrec.cc1 = 0.0;
1739  satrec.cc4 = 0.0;
1740  satrec.cc5 = 0.0;
1741  satrec.d2 = 0.0;
1742  satrec.d3 = 0.0;
1743  satrec.d4 = 0.0;
1744  satrec.delmo = 0.0;
1745  satrec.eta = 0.0;
1746  satrec.argpdot = 0.0;
1747  satrec.omgcof = 0.0;
1748  satrec.sinmao = 0.0;
1749  satrec.t = 0.0;
1750  satrec.t2cof = 0.0;
1751  satrec.t3cof = 0.0;
1752  satrec.t4cof = 0.0;
1753  satrec.t5cof = 0.0;
1754  satrec.x1mth2 = 0.0;
1755  satrec.x7thm1 = 0.0;
1756  satrec.mdot = 0.0;
1757  satrec.nodedot = 0.0;
1758  satrec.xlcof = 0.0;
1759  satrec.xmcof = 0.0;
1760  satrec.nodecf = 0.0;
1761 
1762  /* ----------- set all deep space variables to zero ------------ */
1763  satrec.irez = 0;
1764  satrec.d2201 = 0.0;
1765  satrec.d2211 = 0.0;
1766  satrec.d3210 = 0.0;
1767  satrec.d3222 = 0.0;
1768  satrec.d4410 = 0.0;
1769  satrec.d4422 = 0.0;
1770  satrec.d5220 = 0.0;
1771  satrec.d5232 = 0.0;
1772  satrec.d5421 = 0.0;
1773  satrec.d5433 = 0.0;
1774  satrec.dedt = 0.0;
1775  satrec.del1 = 0.0;
1776  satrec.del2 = 0.0;
1777  satrec.del3 = 0.0;
1778  satrec.didt = 0.0;
1779  satrec.dmdt = 0.0;
1780  satrec.dnodt = 0.0;
1781  satrec.domdt = 0.0;
1782  satrec.e3 = 0.0;
1783  satrec.ee2 = 0.0;
1784  satrec.peo = 0.0;
1785  satrec.pgho = 0.0;
1786  satrec.pho = 0.0;
1787  satrec.pinco = 0.0;
1788  satrec.plo = 0.0;
1789  satrec.se2 = 0.0;
1790  satrec.se3 = 0.0;
1791  satrec.sgh2 = 0.0;
1792  satrec.sgh3 = 0.0;
1793  satrec.sgh4 = 0.0;
1794  satrec.sh2 = 0.0;
1795  satrec.sh3 = 0.0;
1796  satrec.si2 = 0.0;
1797  satrec.si3 = 0.0;
1798  satrec.sl2 = 0.0;
1799  satrec.sl3 = 0.0;
1800  satrec.sl4 = 0.0;
1801  satrec.gsto = 0.0;
1802  satrec.xfact = 0.0;
1803  satrec.xgh2 = 0.0;
1804  satrec.xgh3 = 0.0;
1805  satrec.xgh4 = 0.0;
1806  satrec.xh2 = 0.0;
1807  satrec.xh3 = 0.0;
1808  satrec.xi2 = 0.0;
1809  satrec.xi3 = 0.0;
1810  satrec.xl2 = 0.0;
1811  satrec.xl3 = 0.0;
1812  satrec.xl4 = 0.0;
1813  satrec.xlamo = 0.0;
1814  satrec.zmol = 0.0;
1815  satrec.zmos = 0.0;
1816  satrec.atime = 0.0;
1817  satrec.xli = 0.0;
1818  satrec.xni = 0.0;
1819 
1820  // sgp4fix - note the following variables are also passed directly via satrec.
1821  // it is possible to streamline the sgp4init call by deleting the "x"
1822  // variables, but the user would need to set the satrec.* values first. we
1823  // include the additional assignments in case twoline2rv is not used.
1824  satrec.bstar = xbstar;
1825  satrec.ecco = xecco;
1826  satrec.argpo = xargpo;
1827  satrec.inclo = xinclo;
1828  satrec.mo = xmo;
1829  satrec.no = xno;
1830  satrec.nodeo = xnodeo;
1831 
1832  // sgp4fix add opsmode
1833  satrec.operationmode = opsmode;
1834 
1835  /* ------------------------ earth constants ----------------------- */
1836  // sgp4fix identify constants and allow alternate values
1837  getgravconst(whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2);
1838  ss = 78.0 / radiusearthkm + 1.0;
1839  // sgp4fix use multiply for speed instead of pow
1840  qzms2ttemp = (120.0 - 78.0) / radiusearthkm;
1841  qzms2t = qzms2ttemp * qzms2ttemp * qzms2ttemp * qzms2ttemp;
1842  x2o3 = 2.0 / 3.0;
1843 
1844  satrec.init = 'y';
1845  satrec.t = 0.0;
1846 
1847  initl(satn,
1848  whichconst,
1849  satrec.ecco,
1850  epoch,
1851  satrec.inclo,
1852  satrec.no,
1853  satrec.method,
1854  ainv,
1855  ao,
1856  satrec.con41,
1857  con42,
1858  cosio,
1859  cosio2,
1860  eccsq,
1861  omeosq,
1862  posq,
1863  rp,
1864  rteosq,
1865  sinio,
1866  satrec.gsto,
1867  satrec.operationmode);
1868  satrec.error = 0;
1869 
1870  // sgp4fix remove this check as it is unnecessary
1871  // the mrt check in sgp4 handles decaying satellite cases even if the starting
1872  // condition is below the surface of te earth
1873  // if (rp < 1.0)
1874  // {
1875  // printf("# *** satn%d epoch elts sub-orbital ***\n", satn);
1876  // satrec.error = 5;
1877  // }
1878 
1879  if ((omeosq >= 0.0) || (satrec.no >= 0.0))
1880  {
1881  satrec.isimp = 0;
1882  if (rp < (220.0 / radiusearthkm + 1.0))
1883  satrec.isimp = 1;
1884  sfour = ss;
1885  qzms24 = qzms2t;
1886  perige = (rp - 1.0) * radiusearthkm;
1887 
1888  /* - for perigees below 156 km, s and qoms2t are altered - */
1889  if (perige < 156.0)
1890  {
1891  sfour = perige - 78.0;
1892  if (perige < 98.0)
1893  sfour = 20.0;
1894  // sgp4fix use multiply for speed instead of pow
1895  qzms24temp = (120.0 - sfour) / radiusearthkm;
1896  qzms24 = qzms24temp * qzms24temp * qzms24temp * qzms24temp;
1897  sfour = sfour / radiusearthkm + 1.0;
1898  }
1899  pinvsq = 1.0 / posq;
1900 
1901  tsi = 1.0 / (ao - sfour);
1902  satrec.eta = ao * satrec.ecco * tsi;
1903  etasq = satrec.eta * satrec.eta;
1904  eeta = satrec.ecco * satrec.eta;
1905  psisq = fabs(1.0 - etasq);
1906  coef = qzms24 * pow(tsi, 4.0);
1907  coef1 = coef / pow(psisq, 3.5);
1908  cc2 = coef1 * satrec.no *
1909  (ao * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) +
1910  0.375 * j2 * tsi / psisq * satrec.con41 * (8.0 + 3.0 * etasq * (8.0 + etasq)));
1911  satrec.cc1 = satrec.bstar * cc2;
1912  cc3 = 0.0;
1913  if (satrec.ecco > 1.0e-4)
1914  cc3 = -2.0 * coef * tsi * j3oj2 * satrec.no * sinio / satrec.ecco;
1915  satrec.x1mth2 = 1.0 - cosio2;
1916  satrec.cc4 = 2.0 * satrec.no * coef1 * ao * omeosq *
1917  (satrec.eta * (2.0 + 0.5 * etasq) + satrec.ecco * (0.5 + 2.0 * etasq) -
1918  j2 * tsi / (ao * psisq) *
1919  (-3.0 * satrec.con41 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta)) +
1920  0.75 * satrec.x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) *
1921  cos(2.0 * satrec.argpo)));
1922  satrec.cc5 = 2.0 * coef1 * ao * omeosq * (1.0 + 2.75 * (etasq + eeta) + eeta * etasq);
1923  cosio4 = cosio2 * cosio2;
1924  temp1 = 1.5 * j2 * pinvsq * satrec.no;
1925  temp2 = 0.5 * temp1 * j2 * pinvsq;
1926  temp3 = -0.46875 * j4 * pinvsq * pinvsq * satrec.no;
1927  satrec.mdot = satrec.no + 0.5 * temp1 * rteosq * satrec.con41 +
1928  0.0625 * temp2 * rteosq * (13.0 - 78.0 * cosio2 + 137.0 * cosio4);
1929  satrec.argpdot = -0.5 * temp1 * con42 +
1930  0.0625 * temp2 * (7.0 - 114.0 * cosio2 + 395.0 * cosio4) +
1931  temp3 * (3.0 - 36.0 * cosio2 + 49.0 * cosio4);
1932  xhdot1 = -temp1 * cosio;
1933  satrec.nodedot =
1934  xhdot1 +
1935  (0.5 * temp2 * (4.0 - 19.0 * cosio2) + 2.0 * temp3 * (3.0 - 7.0 * cosio2)) * cosio;
1936  xpidot = satrec.argpdot + satrec.nodedot;
1937  satrec.omgcof = satrec.bstar * cc3 * cos(satrec.argpo);
1938  satrec.xmcof = 0.0;
1939  if (satrec.ecco > 1.0e-4)
1940  satrec.xmcof = -x2o3 * coef * satrec.bstar / eeta;
1941  satrec.nodecf = 3.5 * omeosq * xhdot1 * satrec.cc1;
1942  satrec.t2cof = 1.5 * satrec.cc1;
1943  // sgp4fix for divide by zero with xinco = 180 deg
1944  if (fabs(cosio + 1.0) > 1.5e-12)
1945  satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / (1.0 + cosio);
1946  else
1947  satrec.xlcof = -0.25 * j3oj2 * sinio * (3.0 + 5.0 * cosio) / temp4;
1948  satrec.aycof = -0.5 * j3oj2 * sinio;
1949  // sgp4fix use multiply for speed instead of pow
1950  delmotemp = 1.0 + satrec.eta * cos(satrec.mo);
1951  satrec.delmo = delmotemp * delmotemp * delmotemp;
1952  satrec.sinmao = sin(satrec.mo);
1953  satrec.x7thm1 = 7.0 * cosio2 - 1.0;
1954 
1955  /* --------------- deep space initialization ------------- */
1956  if ((2 * pi / satrec.no) >= 225.0)
1957  {
1958  satrec.method = 'd';
1959  satrec.isimp = 1;
1960  tc = 0.0;
1961  inclm = satrec.inclo;
1962 
1963  dscom(epoch,
1964  satrec.ecco,
1965  satrec.argpo,
1966  tc,
1967  satrec.inclo,
1968  satrec.nodeo,
1969  satrec.no,
1970  snodm,
1971  cnodm,
1972  sinim,
1973  cosim,
1974  sinomm,
1975  cosomm,
1976  day,
1977  satrec.e3,
1978  satrec.ee2,
1979  em,
1980  emsq,
1981  gam,
1982  satrec.peo,
1983  satrec.pgho,
1984  satrec.pho,
1985  satrec.pinco,
1986  satrec.plo,
1987  rtemsq,
1988  satrec.se2,
1989  satrec.se3,
1990  satrec.sgh2,
1991  satrec.sgh3,
1992  satrec.sgh4,
1993  satrec.sh2,
1994  satrec.sh3,
1995  satrec.si2,
1996  satrec.si3,
1997  satrec.sl2,
1998  satrec.sl3,
1999  satrec.sl4,
2000  s1,
2001  s2,
2002  s3,
2003  s4,
2004  s5,
2005  s6,
2006  s7,
2007  ss1,
2008  ss2,
2009  ss3,
2010  ss4,
2011  ss5,
2012  ss6,
2013  ss7,
2014  sz1,
2015  sz2,
2016  sz3,
2017  sz11,
2018  sz12,
2019  sz13,
2020  sz21,
2021  sz22,
2022  sz23,
2023  sz31,
2024  sz32,
2025  sz33,
2026  satrec.xgh2,
2027  satrec.xgh3,
2028  satrec.xgh4,
2029  satrec.xh2,
2030  satrec.xh3,
2031  satrec.xi2,
2032  satrec.xi3,
2033  satrec.xl2,
2034  satrec.xl3,
2035  satrec.xl4,
2036  nm,
2037  z1,
2038  z2,
2039  z3,
2040  z11,
2041  z12,
2042  z13,
2043  z21,
2044  z22,
2045  z23,
2046  z31,
2047  z32,
2048  z33,
2049  satrec.zmol,
2050  satrec.zmos);
2051  dpper(satrec.e3,
2052  satrec.ee2,
2053  satrec.peo,
2054  satrec.pgho,
2055  satrec.pho,
2056  satrec.pinco,
2057  satrec.plo,
2058  satrec.se2,
2059  satrec.se3,
2060  satrec.sgh2,
2061  satrec.sgh3,
2062  satrec.sgh4,
2063  satrec.sh2,
2064  satrec.sh3,
2065  satrec.si2,
2066  satrec.si3,
2067  satrec.sl2,
2068  satrec.sl3,
2069  satrec.sl4,
2070  satrec.t,
2071  satrec.xgh2,
2072  satrec.xgh3,
2073  satrec.xgh4,
2074  satrec.xh2,
2075  satrec.xh3,
2076  satrec.xi2,
2077  satrec.xi3,
2078  satrec.xl2,
2079  satrec.xl3,
2080  satrec.xl4,
2081  satrec.zmol,
2082  satrec.zmos,
2083  inclm,
2084  satrec.init,
2085  satrec.ecco,
2086  satrec.inclo,
2087  satrec.nodeo,
2088  satrec.argpo,
2089  satrec.mo,
2090  satrec.operationmode);
2091 
2092  argpm = 0.0;
2093  nodem = 0.0;
2094  mm = 0.0;
2095 
2096  dsinit(whichconst,
2097  cosim,
2098  emsq,
2099  satrec.argpo,
2100  s1,
2101  s2,
2102  s3,
2103  s4,
2104  s5,
2105  sinim,
2106  ss1,
2107  ss2,
2108  ss3,
2109  ss4,
2110  ss5,
2111  sz1,
2112  sz3,
2113  sz11,
2114  sz13,
2115  sz21,
2116  sz23,
2117  sz31,
2118  sz33,
2119  satrec.t,
2120  tc,
2121  satrec.gsto,
2122  satrec.mo,
2123  satrec.mdot,
2124  satrec.no,
2125  satrec.nodeo,
2126  satrec.nodedot,
2127  xpidot,
2128  z1,
2129  z3,
2130  z11,
2131  z13,
2132  z21,
2133  z23,
2134  z31,
2135  z33,
2136  satrec.ecco,
2137  eccsq,
2138  em,
2139  argpm,
2140  inclm,
2141  mm,
2142  nm,
2143  nodem,
2144  satrec.irez,
2145  satrec.atime,
2146  satrec.d2201,
2147  satrec.d2211,
2148  satrec.d3210,
2149  satrec.d3222,
2150  satrec.d4410,
2151  satrec.d4422,
2152  satrec.d5220,
2153  satrec.d5232,
2154  satrec.d5421,
2155  satrec.d5433,
2156  satrec.dedt,
2157  satrec.didt,
2158  satrec.dmdt,
2159  dndt,
2160  satrec.dnodt,
2161  satrec.domdt,
2162  satrec.del1,
2163  satrec.del2,
2164  satrec.del3,
2165  satrec.xfact,
2166  satrec.xlamo,
2167  satrec.xli,
2168  satrec.xni);
2169  }
2170 
2171  /* ----------- set variables if not deep space ----------- */
2172  if (satrec.isimp != 1)
2173  {
2174  cc1sq = satrec.cc1 * satrec.cc1;
2175  satrec.d2 = 4.0 * ao * tsi * cc1sq;
2176  temp = satrec.d2 * tsi * satrec.cc1 / 3.0;
2177  satrec.d3 = (17.0 * ao + sfour) * temp;
2178  satrec.d4 = 0.5 * temp * ao * tsi * (221.0 * ao + 31.0 * sfour) * satrec.cc1;
2179  satrec.t3cof = satrec.d2 + 2.0 * cc1sq;
2180  satrec.t4cof =
2181  0.25 * (3.0 * satrec.d3 + satrec.cc1 * (12.0 * satrec.d2 + 10.0 * cc1sq));
2182  satrec.t5cof =
2183  0.2 * (3.0 * satrec.d4 + 12.0 * satrec.cc1 * satrec.d3 +
2184  6.0 * satrec.d2 * satrec.d2 + 15.0 * cc1sq * (2.0 * satrec.d2 + cc1sq));
2185  }
2186  } // if omeosq = 0 ...
2187 
2188  /* finally propogate to zero epoch to initialize all others. */
2189  // sgp4fix take out check to let satellites process until they are actually below earth surface
2190  // if(satrec.error == 0)
2191  sgp4(whichconst, satrec, 0.0, r, v);
2192 
2193  satrec.init = 'n';
2194 
2195  return true;
2196 } // end sgp4init
2197 
2198 /*-----------------------------------------------------------------------------
2199  *
2200  * procedure sgp4
2201  *
2202  * this procedure is the sgp4 prediction model from space command. this is an
2203  * updated and combined version of sgp4 and sdp4, which were originally
2204  * published separately in spacetrack report #3. this version follows the
2205  * methodology from the aiaa paper (2006) describing the history and
2206  * development of the code.
2207  *
2208  * author : david vallado 719-573-2600 28 jun 2005
2209  *
2210  * inputs :
2211  * satrec - initialised structure from sgp4init() call.
2212  * tsince - time eince epoch (minutes)
2213  *
2214  * outputs :
2215  * r - position vector km
2216  * v - velocity km/sec
2217  * return code - non-zero on error.
2218  * 1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er
2219  * 2 - mean motion less than 0.0
2220  * 3 - pert elements, ecc < 0.0 or ecc > 1.0
2221  * 4 - semi-latus rectum < 0.0
2222  * 5 - epoch elements are sub-orbital
2223  * 6 - satellite has decayed
2224  *
2225  * locals :
2226  * am -
2227  * axnl, aynl -
2228  * betal -
2229  * cosim , sinim , cosomm , sinomm , cnod , snod , cos2u ,
2230  * sin2u , coseo1 , sineo1 , cosi , sini , cosip , sinip ,
2231  * cosisq , cossu , sinsu , cosu , sinu
2232  * delm -
2233  * delomg -
2234  * dndt -
2235  * eccm -
2236  * emsq -
2237  * ecose -
2238  * el2 -
2239  * eo1 -
2240  * eccp -
2241  * esine -
2242  * argpm -
2243  * argpp -
2244  * omgadf -
2245  * pl -
2246  * r -
2247  * rtemsq -
2248  * rdotl -
2249  * rl -
2250  * rvdot -
2251  * rvdotl -
2252  * su -
2253  * t2 , t3 , t4 , tc
2254  * tem5, temp , temp1 , temp2 , tempa , tempe , templ
2255  * u , ux , uy , uz , vx , vy , vz
2256  * inclm - inclination
2257  * mm - mean anomaly
2258  * nm - mean motion
2259  * nodem - right asc of ascending node
2260  * xinc -
2261  * xincp -
2262  * xl -
2263  * xlm -
2264  * mp -
2265  * xmdf -
2266  * xmx -
2267  * xmy -
2268  * nodedf -
2269  * xnode -
2270  * nodep -
2271  * np -
2272  *
2273  * coupling :
2274  * getgravconst-
2275  * dpper
2276  * dpspace
2277  *
2278  * references :
2279  * hoots, roehrich, norad spacetrack report #3 1980
2280  * hoots, norad spacetrack report #6 1986
2281  * hoots, schumacher and glover 2004
2282  * vallado, crawford, hujsak, kelso 2006
2283  * ----------------------------------------------------------------------------*/
2284 
2285 bool
2286 sgp4(gravconsttype whichconst, elsetrec& satrec, double tsince, double r[3], double v[3])
2287 {
2288  double am, axnl, aynl, betal, cosim, cnod, cos2u, coseo1, cosi, cosip, cosisq, cossu, cosu,
2289  delm, delomg, em, emsq, ecose, el2, eo1, ep, esine, argpm, argpp, argpdf, pl,
2290  mrt = 0.0, mvt, rdotl, rl, rvdot, rvdotl, sinim, sin2u, sineo1, sini, sinip, sinsu, sinu,
2291  snod, su, t2, t3, t4, tem5, temp, temp1, temp2, tempa, tempe, templ, u, ux, uy, uz, vx, vy,
2292  vz, inclm, mm, nm, nodem, xinc, xincp, xl, xlm, mp, xmdf, xmx, xmy, nodedf, xnode, nodep,
2293  tc, dndt, twopi, x2o3, j2, j3, tumin, j4, xke, j3oj2, radiusearthkm, mu, vkmpersec,
2294  delmtemp;
2295  int ktr;
2296 
2297  // ensure that all variables are initialized to make the compiler happy
2298  am = axnl = aynl = betal = cosim = cnod = cos2u = coseo1 = cosi = cosip = cosisq = cossu =
2299  cosu = delm = delomg = em = emsq = ecose = el2 = eo1 = ep = esine = argpm = argpp = argpdf =
2300  pl = mvt = rdotl = rl = rvdot = rvdotl = sinim = sin2u = sineo1 = sini = sinip = sinsu =
2301  sinu = snod = su = t2 = t3 = t4 = tem5 = temp = temp1 = temp2 = tempa = tempe =
2302  templ = u = ux = uy = uz = vx = vy = vz = inclm = mm = nm = nodem = xinc =
2303  xincp = xl = xlm = mp = xmdf = xmx = xmy = nodedf = xnode = nodep = tc =
2304  dndt = twopi = x2o3 = j2 = j3 = tumin = j4 = xke = j3oj2 =
2305  radiusearthkm = mu = vkmpersec = delmtemp = 0.0;
2306 
2307  /* ------------------ set mathematical constants --------------- */
2308  // sgp4fix divisor for divide by zero check on inclination
2309  // the old check used 1.0 + cos(pi-1.0e-9), but then compared it to
2310  // 1.5 e-12, so the threshold was changed to 1.5e-12 for consistency
2311  const double temp4 = 1.5e-12;
2312  twopi = 2.0 * pi;
2313  x2o3 = 2.0 / 3.0;
2314  // sgp4fix identify constants and allow alternate values
2315  getgravconst(whichconst, tumin, mu, radiusearthkm, xke, j2, j3, j4, j3oj2);
2316  vkmpersec = radiusearthkm * xke / 60.0;
2317 
2318  /* --------------------- clear sgp4 error flag ----------------- */
2319  satrec.t = tsince;
2320  satrec.error = 0;
2321 
2322  /* ------- update for secular gravity and atmospheric drag ----- */
2323  xmdf = satrec.mo + satrec.mdot * satrec.t;
2324  argpdf = satrec.argpo + satrec.argpdot * satrec.t;
2325  nodedf = satrec.nodeo + satrec.nodedot * satrec.t;
2326  argpm = argpdf;
2327  mm = xmdf;
2328  t2 = satrec.t * satrec.t;
2329  nodem = nodedf + satrec.nodecf * t2;
2330  tempa = 1.0 - satrec.cc1 * satrec.t;
2331  tempe = satrec.bstar * satrec.cc4 * satrec.t;
2332  templ = satrec.t2cof * t2;
2333 
2334  if (satrec.isimp != 1)
2335  {
2336  delomg = satrec.omgcof * satrec.t;
2337  // sgp4fix use mutliply for speed instead of pow
2338  delmtemp = 1.0 + satrec.eta * cos(xmdf);
2339  delm = satrec.xmcof * (delmtemp * delmtemp * delmtemp - satrec.delmo);
2340  temp = delomg + delm;
2341  mm = xmdf + temp;
2342  argpm = argpdf - temp;
2343  t3 = t2 * satrec.t;
2344  t4 = t3 * satrec.t;
2345  tempa = tempa - satrec.d2 * t2 - satrec.d3 * t3 - satrec.d4 * t4;
2346  tempe = tempe + satrec.bstar * satrec.cc5 * (sin(mm) - satrec.sinmao);
2347  templ = templ + satrec.t3cof * t3 + t4 * (satrec.t4cof + satrec.t * satrec.t5cof);
2348  }
2349 
2350  nm = satrec.no;
2351  em = satrec.ecco;
2352  inclm = satrec.inclo;
2353  if (satrec.method == 'd')
2354  {
2355  tc = satrec.t;
2356  dspace(satrec.irez,
2357  satrec.d2201,
2358  satrec.d2211,
2359  satrec.d3210,
2360  satrec.d3222,
2361  satrec.d4410,
2362  satrec.d4422,
2363  satrec.d5220,
2364  satrec.d5232,
2365  satrec.d5421,
2366  satrec.d5433,
2367  satrec.dedt,
2368  satrec.del1,
2369  satrec.del2,
2370  satrec.del3,
2371  satrec.didt,
2372  satrec.dmdt,
2373  satrec.dnodt,
2374  satrec.domdt,
2375  satrec.argpo,
2376  satrec.argpdot,
2377  satrec.t,
2378  tc,
2379  satrec.gsto,
2380  satrec.xfact,
2381  satrec.xlamo,
2382  satrec.no,
2383  satrec.atime,
2384  em,
2385  argpm,
2386  inclm,
2387  satrec.xli,
2388  mm,
2389  satrec.xni,
2390  nodem,
2391  dndt,
2392  nm);
2393  } // if method = d
2394 
2395  if (nm <= 0.0)
2396  {
2397  satrec.error = 2;
2398  // sgp4fix add return
2399  return false;
2400  }
2401  am = pow((xke / nm), x2o3) * tempa * tempa;
2402  nm = xke / pow(am, 1.5);
2403  em = em - tempe;
2404 
2405  // fix tolerance for error recognition
2406  // sgp4fix am is fixed from the previous nm check
2407  if ((em >= 1.0) || (em < -0.001) /* || (am < 0.95)*/)
2408  {
2409  satrec.error = 1;
2410  // sgp4fix to return if there is an error in eccentricity
2411  return false;
2412  }
2413  // sgp4fix fix tolerance to avoid a divide by zero
2414  if (em < 1.0e-6)
2415  em = 1.0e-6;
2416  mm = mm + satrec.no * templ;
2417  xlm = mm + argpm + nodem;
2418  emsq = em * em;
2419  temp = 1.0 - emsq;
2420 
2421  nodem = fmod(nodem, twopi);
2422  argpm = fmod(argpm, twopi);
2423  xlm = fmod(xlm, twopi);
2424  mm = fmod(xlm - argpm - nodem, twopi);
2425 
2426  /* ----------------- compute extra mean quantities ------------- */
2427  sinim = sin(inclm);
2428  cosim = cos(inclm);
2429 
2430  /* -------------------- add lunar-solar periodics -------------- */
2431  ep = em;
2432  xincp = inclm;
2433  argpp = argpm;
2434  nodep = nodem;
2435  mp = mm;
2436  sinip = sinim;
2437  cosip = cosim;
2438  if (satrec.method == 'd')
2439  {
2440  dpper(satrec.e3,
2441  satrec.ee2,
2442  satrec.peo,
2443  satrec.pgho,
2444  satrec.pho,
2445  satrec.pinco,
2446  satrec.plo,
2447  satrec.se2,
2448  satrec.se3,
2449  satrec.sgh2,
2450  satrec.sgh3,
2451  satrec.sgh4,
2452  satrec.sh2,
2453  satrec.sh3,
2454  satrec.si2,
2455  satrec.si3,
2456  satrec.sl2,
2457  satrec.sl3,
2458  satrec.sl4,
2459  satrec.t,
2460  satrec.xgh2,
2461  satrec.xgh3,
2462  satrec.xgh4,
2463  satrec.xh2,
2464  satrec.xh3,
2465  satrec.xi2,
2466  satrec.xi3,
2467  satrec.xl2,
2468  satrec.xl3,
2469  satrec.xl4,
2470  satrec.zmol,
2471  satrec.zmos,
2472  satrec.inclo,
2473  'n',
2474  ep,
2475  xincp,
2476  nodep,
2477  argpp,
2478  mp,
2479  satrec.operationmode);
2480  if (xincp < 0.0)
2481  {
2482  xincp = -xincp;
2483  nodep = nodep + pi;
2484  argpp = argpp - pi;
2485  }
2486  if ((ep < 0.0) || (ep > 1.0))
2487  {
2488  satrec.error = 3;
2489  // sgp4fix add return
2490  return false;
2491  }
2492  } // if method = d
2493 
2494  /* -------------------- long period periodics ------------------ */
2495  if (satrec.method == 'd')
2496  {
2497  sinip = sin(xincp);
2498  cosip = cos(xincp);
2499  satrec.aycof = -0.5 * j3oj2 * sinip;
2500  // sgp4fix for divide by zero for xincp = 180 deg
2501  if (fabs(cosip + 1.0) > 1.5e-12)
2502  satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / (1.0 + cosip);
2503  else
2504  satrec.xlcof = -0.25 * j3oj2 * sinip * (3.0 + 5.0 * cosip) / temp4;
2505  }
2506  axnl = ep * cos(argpp);
2507  temp = 1.0 / (am * (1.0 - ep * ep));
2508  aynl = ep * sin(argpp) + temp * satrec.aycof;
2509  xl = mp + argpp + nodep + temp * satrec.xlcof * axnl;
2510 
2511  /* --------------------- solve kepler's equation --------------- */
2512  u = fmod(xl - nodep, twopi);
2513  eo1 = u;
2514  tem5 = 9999.9;
2515  ktr = 1;
2516  // sgp4fix for kepler iteration
2517  // the following iteration needs better limits on corrections
2518  while ((fabs(tem5) >= 1.0e-12) && (ktr <= 10))
2519  {
2520  sineo1 = sin(eo1);
2521  coseo1 = cos(eo1);
2522  tem5 = 1.0 - coseo1 * axnl - sineo1 * aynl;
2523  tem5 = (u - aynl * coseo1 + axnl * sineo1 - eo1) / tem5;
2524  if (fabs(tem5) >= 0.95)
2525  tem5 = tem5 > 0.0 ? 0.95 : -0.95;
2526  eo1 = eo1 + tem5;
2527  ktr = ktr + 1;
2528  }
2529 
2530  /* ------------- short period preliminary quantities ----------- */
2531  ecose = axnl * coseo1 + aynl * sineo1;
2532  esine = axnl * sineo1 - aynl * coseo1;
2533  el2 = axnl * axnl + aynl * aynl;
2534  pl = am * (1.0 - el2);
2535  if (pl < 0.0)
2536  {
2537  satrec.error = 4;
2538  // sgp4fix add return
2539  return false;
2540  }
2541  else
2542  {
2543  rl = am * (1.0 - ecose);
2544  rdotl = sqrt(am) * esine / rl;
2545  rvdotl = sqrt(pl) / rl;
2546  betal = sqrt(1.0 - el2);
2547  temp = esine / (1.0 + betal);
2548  sinu = am / rl * (sineo1 - aynl - axnl * temp);
2549  cosu = am / rl * (coseo1 - axnl + aynl * temp);
2550  su = atan2(sinu, cosu);
2551  sin2u = (cosu + cosu) * sinu;
2552  cos2u = 1.0 - 2.0 * sinu * sinu;
2553  temp = 1.0 / pl;
2554  temp1 = 0.5 * j2 * temp;
2555  temp2 = temp1 * temp;
2556 
2557  /* -------------- update for short period periodics ------------ */
2558  if (satrec.method == 'd')
2559  {
2560  cosisq = cosip * cosip;
2561  satrec.con41 = 3.0 * cosisq - 1.0;
2562  satrec.x1mth2 = 1.0 - cosisq;
2563  satrec.x7thm1 = 7.0 * cosisq - 1.0;
2564  }
2565  mrt = rl * (1.0 - 1.5 * temp2 * betal * satrec.con41) + 0.5 * temp1 * satrec.x1mth2 * cos2u;
2566  su = su - 0.25 * temp2 * satrec.x7thm1 * sin2u;
2567  xnode = nodep + 1.5 * temp2 * cosip * sin2u;
2568  xinc = xincp + 1.5 * temp2 * cosip * sinip * cos2u;
2569  mvt = rdotl - nm * temp1 * satrec.x1mth2 * sin2u / xke;
2570  rvdot = rvdotl + nm * temp1 * (satrec.x1mth2 * cos2u + 1.5 * satrec.con41) / xke;
2571 
2572  /* --------------------- orientation vectors ------------------- */
2573  sinsu = sin(su);
2574  cossu = cos(su);
2575  snod = sin(xnode);
2576  cnod = cos(xnode);
2577  sini = sin(xinc);
2578  cosi = cos(xinc);
2579  xmx = -snod * cosi;
2580  xmy = cnod * cosi;
2581  ux = xmx * sinsu + cnod * cossu;
2582  uy = xmy * sinsu + snod * cossu;
2583  uz = sini * sinsu;
2584  vx = xmx * cossu - cnod * sinsu;
2585  vy = xmy * cossu - snod * sinsu;
2586  vz = sini * cossu;
2587 
2588  /* --------- position and velocity (in km and km/sec) ---------- */
2589  r[0] = (mrt * ux) * radiusearthkm;
2590  r[1] = (mrt * uy) * radiusearthkm;
2591  r[2] = (mrt * uz) * radiusearthkm;
2592  v[0] = (mvt * ux + rvdot * vx) * vkmpersec;
2593  v[1] = (mvt * uy + rvdot * vy) * vkmpersec;
2594  v[2] = (mvt * uz + rvdot * vz) * vkmpersec;
2595  } // if pl > 0
2596 
2597  // sgp4fix for decaying satellites
2598  if (mrt < 1.0)
2599  {
2600  satrec.error = 6;
2601  return false;
2602  }
2603 
2604  return true;
2605 } // end sgp4
2606 
2607 /* -----------------------------------------------------------------------------
2608  *
2609  * function gstime
2610  *
2611  * this function finds the greenwich sidereal time.
2612  *
2613  * author : david vallado 719-573-2600 1 mar 2001
2614  *
2615  * inputs description range / units
2616  * jdut1 - julian date in ut1 days from 4713 bc
2617  *
2618  * outputs :
2619  * gstime - greenwich sidereal time 0 to 2pi rad
2620  *
2621  * locals :
2622  * temp - temporary variable for doubles rad
2623  * tut1 - julian centuries from the
2624  * jan 1, 2000 12 h epoch (ut1)
2625  *
2626  * coupling :
2627  * none
2628  *
2629  * references :
2630  * vallado 2004, 191, eq 3-45
2631  * --------------------------------------------------------------------------- */
2632 
2633 double
2634 gstime(double jdut1)
2635 {
2636  const double twopi = 2.0 * pi;
2637  const double deg2rad = pi / 180.0;
2638  double temp, tut1;
2639 
2640  tut1 = (jdut1 - 2451545.0) / 36525.0;
2641  temp = -6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 +
2642  (876600.0 * 3600 + 8640184.812866) * tut1 + 67310.54841; // sec
2643  temp = fmod(temp * deg2rad / 240.0, twopi); // 360/86400 = 1/240, to deg, to rad
2644 
2645  // ------------------------ check quadrants ---------------------
2646  if (temp < 0.0)
2647  temp += twopi;
2648 
2649  return temp;
2650 } // end gstime
2651 
2652 /* -----------------------------------------------------------------------------
2653  *
2654  * function getgravconst
2655  *
2656  * this function gets constants for the propagator. note that mu is identified to
2657  * facilitiate comparisons with newer models. the common useage is wgs72.
2658  *
2659  * author : david vallado 719-573-2600 21 jul 2006
2660  *
2661  * inputs :
2662  * whichconst - which set of constants to use wgs72old, wgs72, wgs84
2663  *
2664  * outputs :
2665  * tumin - minutes in one time unit
2666  * mu - earth gravitational parameter
2667  * radiusearthkm - radius of the earth in km
2668  * xke - reciprocal of tumin
2669  * j2, j3, j4 - un-normalized zonal harmonic values
2670  * j3oj2 - j3 divided by j2
2671  *
2672  * locals :
2673  *
2674  * coupling :
2675  * none
2676  *
2677  * references :
2678  * norad spacetrack report #3
2679  * vallado, crawford, hujsak, kelso 2006
2680  * --------------------------------------------------------------------------- */
2681 
2682 void
2683 getgravconst(gravconsttype whichconst,
2684  double& tumin,
2685  double& mu,
2686  double& radiusearthkm,
2687  double& xke,
2688  double& j2,
2689  double& j3,
2690  double& j4,
2691  double& j3oj2)
2692 {
2693  switch (whichconst)
2694  {
2695  // -- wgs-72 low precision str#3 constants --
2696  case wgs72old:
2697  mu = 398600.79964; // in km3 / s2
2698  radiusearthkm = 6378.135; // km
2699  xke = 0.0743669161;
2700  tumin = 1.0 / xke;
2701  j2 = 0.001082616;
2702  j3 = -0.00000253881;
2703  j4 = -0.00000165597;
2704  j3oj2 = j3 / j2;
2705  break;
2706  // ------------ wgs-72 constants ------------
2707  case wgs72:
2708  mu = 398600.8; // in km3 / s2
2709  radiusearthkm = 6378.135; // km
2710  xke = 60.0 / sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu);
2711  tumin = 1.0 / xke;
2712  j2 = 0.001082616;
2713  j3 = -0.00000253881;
2714  j4 = -0.00000165597;
2715  j3oj2 = j3 / j2;
2716  break;
2717  // ------------ wgs-84 constants ------------
2718  case wgs84:
2719  mu = 398600.5; // in km3 / s2
2720  radiusearthkm = 6378.137; // km
2721  xke = 60.0 / sqrt(radiusearthkm * radiusearthkm * radiusearthkm / mu);
2722  tumin = 1.0 / xke;
2723  j2 = 0.00108262998905;
2724  j3 = -0.00000253215306;
2725  j4 = -0.00000161098761;
2726  j3oj2 = j3 / j2;
2727  break;
2728  default:
2729  fprintf(stderr, "unknown gravity option (%d)\n", whichconst);
2730  break;
2731  }
2732 } // end getgravconst
dspace
static void dspace(int irez, double d2201, double d2211, double d3210, double d3222, double d4410, double d4422, double d5220, double d5232, double d5421, double d5433, double dedt, double del1, double del2, double del3, double didt, double dmdt, double dnodt, double domdt, double argpo, double argpdot, double t, double tc, double gsto, double xfact, double xlamo, double no, double &atime, double &em, double &argpm, double &inclm, double &xli, double &mm, double &xni, double &nodem, double &dndt, double &nm)
Definition: satellite-sgp4unit.cc:1304
initl
static void initl(long int satn, gravconsttype whichconst, double ecco, double epoch, double inclo, double &no, char &method, double &ainv, double &ao, double &con41, double &con42, double &cosio, double &cosio2, double &eccsq, double &omeosq, double &posq, double &rp, double &rteosq, double &sinio, double &gsto, char opsmode)
Definition: satellite-sgp4unit.cc:1524
sgp4
bool sgp4(gravconsttype whichconst, elsetrec &satrec, double tsince, double r[3], double v[3])
Definition: satellite-sgp4unit.cc:2286
dsinit
static void dsinit(gravconsttype whichconst, double cosim, double emsq, double argpo, double s1, double s2, double s3, double s4, double s5, double sinim, double ss1, double ss2, double ss3, double ss4, double ss5, double sz1, double sz3, double sz11, double sz13, double sz21, double sz23, double sz31, double sz33, double t, double tc, double gsto, double mo, double mdot, double no, double nodeo, double nodedot, double xpidot, double z1, double z3, double z11, double z13, double z21, double z23, double z31, double z33, double ecco, double eccsq, double &em, double &argpm, double &inclm, double &mm, double &nm, double &nodem, int &irez, double &atime, double &d2201, double &d2211, double &d3210, double &d3222, double &d4410, double &d4422, double &d5220, double &d5232, double &d5421, double &d5433, double &dedt, double &didt, double &dmdt, double &dndt, double &dnodt, double &domdt, double &del1, double &del2, double &del3, double &xfact, double &xlamo, double &xli, double &xni)
Definition: satellite-sgp4unit.cc:950
getgravconst
void getgravconst(gravconsttype whichconst, double &tumin, double &mu, double &radiusearthkm, double &xke, double &j2, double &j3, double &j4, double &j3oj2)
Definition: satellite-sgp4unit.cc:2683
dpper
static void dpper(double e3, double ee2, double peo, double pgho, double pho, double pinco, double plo, double se2, double se3, double sgh2, double sgh3, double sgh4, double sh2, double sh3, double si2, double si3, double sl2, double sl3, double sl4, double t, double xgh2, double xgh3, double xgh4, double xh2, double xh3, double xi2, double xi3, double xl2, double xl3, double xl4, double zmol, double zmos, double inclo, char init, double &ep, double &inclp, double &nodep, double &argpp, double &mp, char opsmode)
Definition: satellite-sgp4unit.cc:387
dscom
static void dscom(double epoch, double ep, double argpp, double tc, double inclp, double nodep, double np, double &snodm, double &cnodm, double &sinim, double &cosim, double &sinomm, double &cosomm, double &day, double &e3, double &ee2, double &em, double &emsq, double &gam, double &peo, double &pgho, double &pho, double &pinco, double &plo, double &rtemsq, double &se2, double &se3, double &sgh2, double &sgh3, double &sgh4, double &sh2, double &sh3, double &si2, double &si3, double &sl2, double &sl3, double &sl4, double &s1, double &s2, double &s3, double &s4, double &s5, double &s6, double &s7, double &ss1, double &ss2, double &ss3, double &ss4, double &ss5, double &ss6, double &ss7, double &sz1, double &sz2, double &sz3, double &sz11, double &sz12, double &sz13, double &sz21, double &sz22, double &sz23, double &sz31, double &sz32, double &sz33, double &xgh2, double &xgh3, double &xgh4, double &xh2, double &xh3, double &xi2, double &xi3, double &xl2, double &xl3, double &xl4, double &nm, double &z1, double &z2, double &z3, double &z11, double &z12, double &z13, double &z21, double &z22, double &z23, double &z31, double &z32, double &z33, double &zmol, double &zmos)
Definition: satellite-sgp4unit.cc:608
gstime
double gstime(double jdut1)
Definition: satellite-sgp4unit.cc:2634
sgp4init
bool sgp4init(gravconsttype whichconst, char opsmode, const long int satn, const double epoch, const double xbstar, const double xecco, const double xargpo, const double xinclo, const double xmo, const double xno, const double xnodeo, elsetrec &satrec)
Definition: satellite-sgp4unit.cc:1693
pi
#define pi
Definition: satellite-sgp4unit.h:59
gravconsttype
gravconsttype
Definition: satellite-sgp4unit.h:63
wgs72
@ wgs72
Definition: satellite-sgp4unit.h:65
wgs72old
@ wgs72old
Definition: satellite-sgp4unit.h:64
wgs84
@ wgs84
Definition: satellite-sgp4unit.h:66