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