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