/[PAMELA software]/yodaUtility/sgp4/cNoradSDP4.cpp
ViewVC logotype

Annotation of /yodaUtility/sgp4/cNoradSDP4.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Sun Apr 30 11:08:15 2006 UTC (18 years, 7 months ago) by kusanagi
Branch point for: MAIN
Initial revision

1 kusanagi 1.1 //
2     // cNoradSDP4.cpp
3     //
4     // NORAD SDP4 implementation. See historical note in cNoradBase.cpp
5     // Copyright (c) 2003 Michael F. Henry
6     //
7     // mfh 12/07/2003
8     //
9     #include "stdafx.h"
10     #include "cNoradSDP4.h"
11     #include "cTle.h"
12     #include "coord.h"
13     #include "cOrbit.h"
14     #include "cVector.h"
15    
16     const double zns = 1.19459E-5; const double c1ss = 2.9864797E-6;
17     const double zes = 0.01675; const double znl = 1.5835218E-4;
18     const double c1l = 4.7968065E-7; const double zel = 0.05490;
19     const double zcosis = 0.91744867; const double zsinis = 0.39785416;
20     const double zsings = -0.98088458; const double zcosgs = 0.1945905;
21     const double q22 = 1.7891679E-6; const double q31 = 2.1460748E-6;
22     const double q33 = 2.2123015E-7; const double g22 = 5.7686396;
23     const double g32 = 0.95240898; const double g44 = 1.8014998;
24     const double g52 = 1.0508330; const double g54 = 4.4108898;
25     const double root22 = 1.7891679E-6; const double root32 = 3.7393792E-7;
26     const double root44 = 7.3636953E-9; const double root52 = 1.1428639E-7;
27     const double root54 = 2.1765803E-9; const double thdt = 4.3752691E-3;
28    
29     //////////////////////////////////////////////////////////////////////////////
30     cNoradSDP4::cNoradSDP4(const cOrbit &orbit) :
31     cNoradBase(orbit)
32     {
33     m_sing = sin(m_Orbit.ArgPerigee());
34     m_cosg = cos(m_Orbit.ArgPerigee());
35    
36     dp_savtsn = 0.0;
37     dp_zmos = 0.0;
38     dp_se2 = 0.0;
39     dp_se3 = 0.0;
40     dp_si2 = 0.0;
41     dp_si3 = 0.0;
42     dp_sl2 = 0.0;
43     dp_sl3 = 0.0;
44     dp_sl4 = 0.0;
45     dp_sghs = 0.0;
46     dp_sgh2 = 0.0;
47     dp_sgh3 = 0.0;
48     dp_sgh4 = 0.0;
49     dp_sh2 = 0.0;
50     dp_sh3 = 0.0;
51     dp_zmol = 0.0;
52     dp_ee2 = 0.0;
53     dp_e3 = 0.0;
54     dp_xi2 = 0.0;
55     dp_xi3 = 0.0;
56     dp_xl2 = 0.0;
57     dp_xl3 = 0.0;
58     dp_xl4 = 0.0;
59     dp_xgh2 = 0.0;
60     dp_xgh3 = 0.0;
61     dp_xgh4 = 0.0;
62     dp_xh2 = 0.0;
63     dp_xh3 = 0.0;
64     dp_xqncl = 0.0;
65    
66     dp_thgr = 0.0;
67     dp_omegaq = 0.0;
68     dp_sse = 0.0;
69     dp_ssi = 0.0;
70     dp_ssl = 0.0;
71     dp_ssh = 0.0;
72     dp_ssg = 0.0;
73     dp_d2201 = 0.0;
74     dp_d2211 = 0.0;
75     dp_d3210 = 0.0;
76     dp_d3222 = 0.0;
77     dp_d4410 = 0.0;
78     dp_d4422 = 0.0;
79     dp_d5220 = 0.0;
80     dp_d5232 = 0.0;
81     dp_d5421 = 0.0;
82     dp_d5433 = 0.0;
83     dp_xlamo = 0.0;
84     dp_del1 = 0.0;
85     dp_del2 = 0.0;
86     dp_del3 = 0.0;
87     dp_fasx2 = 0.0;
88     dp_fasx4 = 0.0;
89     dp_fasx6 = 0.0;
90     dp_xfact = 0.0;
91     dp_xli = 0.0;
92     dp_xni = 0.0;
93     dp_atime = 0.0;
94     dp_stepp = 0.0;
95     dp_stepn = 0.0;
96     dp_step2 = 0.0;
97    
98     dp_iresfl = false;
99     dp_isynfl = false;
100    
101     }
102    
103     cNoradSDP4::~cNoradSDP4(void)
104     {
105     }
106    
107     //////////////////////////////////////////////////////////////////////////////
108     bool cNoradSDP4::DeepInit(double *eosq, double *sinio, double *cosio,
109     double *betao, double *aodp, double *theta2,
110     double *sing, double *cosg, double *betao2,
111     double *xmdot, double *omgdot, double *xnodott)
112     {
113     eqsq = *eosq;
114     siniq = *sinio;
115     cosiq = *cosio;
116     rteqsq = *betao;
117     ao = *aodp;
118     cosq2 = *theta2;
119     sinomo = *sing;
120     cosomo = *cosg;
121     bsq = *betao2;
122     xlldot = *xmdot;
123     omgdt = *omgdot;
124     xnodot = *xnodott;
125    
126     // Deep space initialization
127     cJulian jd = m_Orbit.Epoch();
128    
129     dp_thgr = jd.toGMST();
130    
131     double eq = m_Orbit.Eccentricity();
132     double aqnv = 1.0 / ao;
133    
134     dp_xqncl = m_Orbit.Inclination();
135    
136     double xmao = m_Orbit.mnAnomaly();
137     double xpidot = omgdt + xnodot;
138     double sinq = sin(m_Orbit.RAAN());
139     double cosq = cos(m_Orbit.RAAN());
140    
141     dp_omegaq = m_Orbit.ArgPerigee();
142    
143     // Initialize lunar solar terms
144     double day = jd.FromJan1_12h_1900();
145    
146     if (day != dpi_day)
147     {
148     dpi_day = day;
149     dpi_xnodce = 4.5236020 - 9.2422029E-4 * day;
150     dpi_stem = sin(dpi_xnodce);
151     dpi_ctem = cos(dpi_xnodce);
152     dpi_zcosil = 0.91375164 - 0.03568096 * dpi_ctem;
153     dpi_zsinil = sqrt(1.0 - dpi_zcosil * dpi_zcosil);
154     dpi_zsinhl = 0.089683511 *dpi_stem / dpi_zsinil;
155     dpi_zcoshl = sqrt(1.0 - dpi_zsinhl * dpi_zsinhl);
156     dpi_c = 4.7199672 + 0.22997150 * day;
157     dpi_gam = 5.8351514 + 0.0019443680 * day;
158     dp_zmol = Fmod2p(dpi_c - dpi_gam);
159     dpi_zx = 0.39785416 * dpi_stem / dpi_zsinil;
160     dpi_zy = dpi_zcoshl * dpi_ctem + 0.91744867 * dpi_zsinhl * dpi_stem;
161     dpi_zx = AcTan(dpi_zx,dpi_zy) + dpi_gam - dpi_xnodce;
162     dpi_zcosgl = cos(dpi_zx);
163     dpi_zsingl = sin(dpi_zx);
164     dp_zmos = 6.2565837 + 0.017201977 * day;
165     dp_zmos = Fmod2p(dp_zmos);
166     }
167    
168     dp_savtsn = 1.0e20;
169    
170     double zcosg = zcosgs;
171     double zsing = zsings;
172     double zcosi = zcosis;
173     double zsini = zsinis;
174     double zcosh = cosq;
175     double zsinh = sinq;
176     double cc = c1ss;
177     double zn = zns;
178     double ze = zes;
179     double zmo = dp_zmos;
180     double xnoi = 1.0 / m_xnodp;
181    
182     double a1; double a3; double a7; double a8; double a9; double a10;
183     double a2; double a4; double a5; double a6; double x1; double x2;
184     double x3; double x4; double x5; double x6; double x7; double x8;
185     double z31; double z32; double z33; double z1; double z2; double z3;
186     double z11; double z12; double z13; double z21; double z22; double z23;
187     double s3; double s2; double s4; double s1; double s5; double s6;
188     double s7;
189     double se = 0.0; double si = 0.0; double sl = 0.0;
190     double sgh = 0.0; double sh = 0.0;
191    
192     // Apply the solar and lunar terms on the first pass, then re-apply the
193     // solar terms again on the second pass.
194    
195     for (int pass = 1; pass <= 2; pass++)
196     {
197     // Do solar terms
198     a1 = zcosg * zcosh + zsing * zcosi * zsinh;
199     a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
200     a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
201     a8 = zsing * zsini;
202     a9 = zsing * zsinh + zcosg * zcosi * zcosh;
203     a10 = zcosg * zsini;
204     a2 = cosiq * a7 + siniq * a8;
205     a4 = cosiq * a9 + siniq * a10;
206     a5 = -siniq * a7 + cosiq * a8;
207     a6 = -siniq * a9 + cosiq * a10;
208     x1 = a1 * cosomo + a2 * sinomo;
209     x2 = a3 * cosomo + a4 * sinomo;
210     x3 = -a1 * sinomo + a2 * cosomo;
211     x4 = -a3 * sinomo + a4 * cosomo;
212     x5 = a5 * sinomo;
213     x6 = a6 * sinomo;
214     x7 = a5 * cosomo;
215     x8 = a6 * cosomo;
216     z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3;
217     z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4;
218     z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4;
219     z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * eqsq;
220     z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * eqsq;
221     z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * eqsq;
222     z11 = -6.0 * a1 * a5 + eqsq*(-24.0 * x1 * x7 - 6.0 * x3 * x5);
223     z12 = -6.0 * (a1 * a6 + a3 * a5) +
224     eqsq * (-24.0 * (x2 * x7 + x1 * x8) - 6.0 * (x3 * x6 + x4 * x5));
225     z13 = -6.0 * a3 * a6 + eqsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6);
226     z21 = 6.0 * a2 * a5 + eqsq * (24.0 * x1 * x5 - 6.0 * x3 * x7);
227     z22 = 6.0*(a4 * a5 + a2 * a6) +
228     eqsq * (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8));
229     z23 = 6.0 * a4 * a6 + eqsq*(24.0 * x2 * x6 - 6.0 * x4 * x8);
230     z1 = z1 + z1 + bsq * z31;
231     z2 = z2 + z2 + bsq * z32;
232     z3 = z3 + z3 + bsq * z33;
233     s3 = cc * xnoi;
234     s2 = -0.5 * s3/rteqsq;
235     s4 = s3 * rteqsq;
236     s1 = -15.0 * eq * s4;
237     s5 = x1 * x3 + x2 * x4;
238     s6 = x2 * x3 + x1 * x4;
239     s7 = x2 * x4 - x1 * x3;
240     se = s1 * zn * s5;
241     si = s2 * zn * (z11 + z13);
242     sl = -zn * s3 * (z1 + z3 - 14.0 - 6.0 * eqsq);
243     sgh = s4 * zn * (z31 + z33 - 6.0);
244     sh = -zn * s2 * (z21 + z23);
245    
246     if (dp_xqncl < 5.2359877E-2)
247     sh = 0.0;
248    
249     dp_ee2 = 2.0 * s1 * s6;
250     dp_e3 = 2.0 * s1 * s7;
251     dp_xi2 = 2.0 * s2 * z12;
252     dp_xi3 = 2.0 * s2 * (z13 - z11);
253     dp_xl2 = -2.0 * s3 * z2;
254     dp_xl3 = -2.0 * s3 * (z3 - z1);
255     dp_xl4 = -2.0 * s3 * (-21.0 - 9.0 * eqsq) * ze;
256     dp_xgh2 = 2.0 * s4 * z32;
257     dp_xgh3 = 2.0 * s4 * (z33 - z31);
258     dp_xgh4 = -18.0 * s4 * ze;
259     dp_xh2 = -2.0 * s2 * z22;
260     dp_xh3 = -2.0 * s2 * (z23 - z21);
261    
262     if (pass == 1)
263     {
264     // Do lunar terms
265     dp_sse = se;
266     dp_ssi = si;
267     dp_ssl = sl;
268     dp_ssh = sh / siniq;
269     dp_ssg = sgh - cosiq * dp_ssh;
270     dp_se2 = dp_ee2;
271     dp_si2 = dp_xi2;
272     dp_sl2 = dp_xl2;
273     dp_sgh2 = dp_xgh2;
274     dp_sh2 = dp_xh2;
275     dp_se3 = dp_e3;
276     dp_si3 = dp_xi3;
277     dp_sl3 = dp_xl3;
278     dp_sgh3 = dp_xgh3;
279     dp_sh3 = dp_xh3;
280     dp_sl4 = dp_xl4;
281     dp_sgh4 = dp_xgh4;
282     zcosg = dpi_zcosgl;
283     zsing = dpi_zsingl;
284     zcosi = dpi_zcosil;
285     zsini = dpi_zsinil;
286     zcosh = dpi_zcoshl * cosq + dpi_zsinhl * sinq;
287     zsinh = sinq * dpi_zcoshl - cosq * dpi_zsinhl;
288     zn = znl;
289     cc = c1l;
290     ze = zel;
291     zmo = dp_zmol;
292     }
293     }
294    
295     dp_sse = dp_sse + se;
296     dp_ssi = dp_ssi + si;
297     dp_ssl = dp_ssl + sl;
298     dp_ssg = dp_ssg + sgh - cosiq / siniq * sh;
299     dp_ssh = dp_ssh + sh / siniq;
300    
301     // Geopotential resonance initialization for 12 hour orbits
302     dp_iresfl = false;
303     dp_isynfl = false;
304    
305     bool bInitOnExit = true;
306     double g310;
307     double f220;
308     double bfact = 0.0;
309    
310     if ((m_xnodp >= 0.0052359877) || (m_xnodp <= 0.0034906585))
311     {
312     if ((m_xnodp < 8.26E-3) || (m_xnodp > 9.24E-3) || (eq < 0.5))
313     {
314     bInitOnExit = false;
315     }
316     else
317     {
318     dp_iresfl = true;
319    
320     double eoc = eq * eqsq;
321     double g201 = -0.306 - (eq - 0.64) * 0.440;
322    
323     double g211; double g322;
324     double g410; double g422;
325     double g520;
326    
327     if (eq <= 0.65)
328     {
329     g211 = 3.616 - 13.247 * eq + 16.290 * eqsq;
330     g310 = -19.302 + 117.390 * eq - 228.419 * eqsq + 156.591 * eoc;
331     g322 = -18.9068 + 109.7927 * eq - 214.6334 * eqsq + 146.5816 * eoc;
332     g410 = -41.122 + 242.694 * eq - 471.094 * eqsq + 313.953 * eoc;
333     g422 = -146.407 + 841.880 * eq - 1629.014 * eqsq + 1083.435 * eoc;
334     g520 = -532.114 + 3017.977 * eq - 5740.0 * eqsq + 3708.276 * eoc;
335     }
336     else
337     {
338     g211 = -72.099 + 331.819 * eq - 508.738 * eqsq + 266.724 * eoc;
339     g310 = -346.844 + 1582.851 * eq - 2415.925 * eqsq + 1246.113 * eoc;
340     g322 = -342.585 + 1554.908 * eq - 2366.899 * eqsq + 1215.972 * eoc;
341     g410 = -1052.797 + 4758.686 * eq - 7193.992 * eqsq + 3651.957 * eoc;
342     g422 = -3581.69 + 16178.11 * eq - 24462.77 * eqsq + 12422.52 * eoc;
343    
344     if (eq <= 0.715)
345     g520 = 1464.74 - 4664.75 * eq + 3763.64 * eqsq;
346     else
347     g520 = -5149.66 + 29936.92 * eq - 54087.36 * eqsq + 31324.56 * eoc;
348     }
349    
350     double g533;
351     double g521;
352     double g532;
353    
354     if (eq < 0.7)
355     {
356     g533 = -919.2277 + 4988.61 * eq - 9064.77 * eqsq + 5542.21 * eoc;
357     g521 = -822.71072 + 4568.6173 * eq - 8491.4146 * eqsq + 5337.524 * eoc;
358     g532 = -853.666 + 4690.25 * eq - 8624.77 * eqsq + 5341.4 * eoc;
359     }
360     else
361     {
362     g533 = -37995.78 + 161616.52 * eq - 229838.2 * eqsq + 109377.94 * eoc;
363     g521 = -51752.104 + 218913.95 * eq - 309468.16 * eqsq + 146349.42 * eoc;
364     g532 = -40023.88 + 170470.89 * eq - 242699.48 * eqsq + 115605.82 * eoc;
365     }
366    
367     double sini2 = siniq * siniq;
368     f220 = 0.75*(1.0 + 2.0 * cosiq + cosq2);
369     double f221 = 1.5 * sini2;
370     double f321 = 1.875 * siniq*(1.0 - 2.0 * cosiq - 3.0 * cosq2);
371     double f322 = -1.875 * siniq*(1.0 + 2.0 * cosiq - 3.0 * cosq2);
372     double f441 = 35.0 * sini2 * f220;
373     double f442 = 39.3750 * sini2 * sini2;
374     double f522 = 9.84375 * siniq*(sini2*(1.0 - 2.0 * cosiq - 5.0 * cosq2) +
375     0.33333333*(-2.0 + 4.0 * cosiq + 6.0 * cosq2));
376     double f523 = siniq*(4.92187512 * sini2*(-2.0 - 4.0 * cosiq + 10.0 * cosq2) +
377     6.56250012 * (1.0 + 2.0 * cosiq - 3.0 * cosq2));
378     double f542 = 29.53125 * siniq*(2.0 - 8.0 * cosiq + cosq2 * (-12.0 + 8.0 * cosiq + 10.0 * cosq2));
379     double f543 = 29.53125 * siniq*(-2.0 - 8.0 * cosiq + cosq2 * (12.0 + 8.0 * cosiq - 10.0 * cosq2));
380     double xno2 = m_xnodp * m_xnodp;
381     double ainv2 = aqnv * aqnv;
382     double temp1 = 3.0 * xno2 * ainv2;
383     double temp = temp1 * root22;
384    
385     dp_d2201 = temp * f220 * g201;
386     dp_d2211 = temp * f221 * g211;
387     temp1 = temp1 * aqnv;
388     temp = temp1 * root32;
389     dp_d3210 = temp * f321 * g310;
390     dp_d3222 = temp * f322 * g322;
391     temp1 = temp1 * aqnv;
392     temp = 2.0 * temp1 * root44;
393     dp_d4410 = temp * f441 * g410;
394     dp_d4422 = temp * f442 * g422;
395     temp1 = temp1 * aqnv;
396     temp = temp1 * root52;
397     dp_d5220 = temp * f522 * g520;
398     dp_d5232 = temp * f523 * g532;
399     temp = 2.0 * temp1 * root54;
400     dp_d5421 = temp * f542 * g521;
401     dp_d5433 = temp * f543 * g533;
402     dp_xlamo = xmao + m_Orbit.RAAN() + m_Orbit.RAAN() - dp_thgr - dp_thgr;
403     bfact = xlldot + xnodot + xnodot - thdt - thdt;
404     bfact = bfact + dp_ssl + dp_ssh + dp_ssh;
405     }
406     }
407     else
408     {
409     // Synchronous resonance terms initialization
410     dp_iresfl = true;
411     dp_isynfl = true;
412     double g200 = 1.0 + eqsq * (-2.5 + 0.8125 * eqsq);
413     g310 = 1.0 + 2.0 * eqsq;
414     double g300 = 1.0 + eqsq * (-6.0 + 6.60937 * eqsq);
415     f220 = 0.75 * (1.0 + cosiq) * (1.0 + cosiq);
416     double f311 = 0.9375 * siniq * siniq * (1.0 + 3 * cosiq) - 0.75 * (1.0 + cosiq);
417     double f330 = 1.0 + cosiq;
418     f330 = 1.875 * f330 * f330 * f330;
419     dp_del1 = 3.0 * m_xnodp * m_xnodp * aqnv * aqnv;
420     dp_del2 = 2.0 * dp_del1 * f220 * g200 * q22;
421     dp_del3 = 3.0 * dp_del1 * f330 * g300 * q33 * aqnv;
422     dp_del1 = dp_del1 * f311 * g310 * q31 * aqnv;
423     dp_fasx2 = 0.13130908;
424     dp_fasx4 = 2.8843198;
425     dp_fasx6 = 0.37448087;
426     dp_xlamo = xmao + m_Orbit.RAAN() + m_Orbit.ArgPerigee() - dp_thgr;
427     bfact = xlldot + xpidot - thdt;
428     bfact = bfact + dp_ssl + dp_ssg + dp_ssh;
429     }
430    
431     if (bInitOnExit)
432     {
433     dp_xfact = bfact - m_xnodp;
434    
435     // Initialize integrator
436     dp_xli = dp_xlamo;
437     dp_xni = m_xnodp;
438     dp_atime = 0.0;
439     dp_stepp = 720.0;
440     dp_stepn = -720.0;
441     dp_step2 = 259200.0;
442     }
443    
444     *eosq = eqsq;
445     *sinio = siniq;
446     *cosio = cosiq;
447     *betao = rteqsq;
448     *aodp = ao;
449     *theta2 = cosq2;
450     *sing = sinomo;
451     *cosg = cosomo;
452     *betao2 = bsq;
453     *xmdot = xlldot;
454     *omgdot = omgdt;
455     *xnodott = xnodot;
456    
457     return true;
458     }
459    
460     //////////////////////////////////////////////////////////////////////////////
461     bool cNoradSDP4::DeepCalcDotTerms(double *pxndot, double *pxnddt, double *pxldot)
462     {
463     // Dot terms calculated
464     if (dp_isynfl)
465     {
466     *pxndot = dp_del1 * sin(dp_xli - dp_fasx2) +
467     dp_del2 * sin(2.0 * (dp_xli - dp_fasx4)) +
468     dp_del3 * sin(3.0 * (dp_xli - dp_fasx6));
469     *pxnddt = dp_del1 * cos(dp_xli - dp_fasx2) +
470     2.0 * dp_del2 * cos(2.0 * (dp_xli - dp_fasx4)) +
471     3.0 * dp_del3 * cos(3.0 * (dp_xli - dp_fasx6));
472     }
473     else
474     {
475     double xomi = dp_omegaq + omgdt * dp_atime;
476     double x2omi = xomi + xomi;
477     double x2li = dp_xli + dp_xli;
478    
479     *pxndot = dp_d2201 * sin(x2omi + dp_xli - g22) +
480     dp_d2211 * sin(dp_xli - g22) +
481     dp_d3210 * sin(xomi + dp_xli - g32) +
482     dp_d3222 * sin(-xomi + dp_xli - g32) +
483     dp_d4410 * sin(x2omi + x2li - g44) +
484     dp_d4422 * sin(x2li - g44) +
485     dp_d5220 * sin(xomi + dp_xli - g52) +
486     dp_d5232 * sin(-xomi + dp_xli - g52) +
487     dp_d5421 * sin(xomi + x2li - g54) +
488     dp_d5433 * sin(-xomi + x2li - g54);
489    
490     *pxnddt = dp_d2201 * cos(x2omi + dp_xli - g22) +
491     dp_d2211 * cos(dp_xli - g22) +
492     dp_d3210 * cos(xomi + dp_xli - g32) +
493     dp_d3222 * cos(-xomi + dp_xli - g32) +
494     dp_d5220 * cos(xomi + dp_xli - g52) +
495     dp_d5232 * cos(-xomi + dp_xli - g52) +
496     2.0 * (dp_d4410 * cos(x2omi + x2li - g44) +
497     dp_d4422 * cos(x2li - g44) +
498     dp_d5421 * cos(xomi + x2li - g54) +
499     dp_d5433 * cos(-xomi + x2li - g54));
500     }
501    
502     *pxldot = dp_xni + dp_xfact;
503     *pxnddt = (*pxnddt) * (*pxldot);
504    
505     return true;
506     }
507    
508     //////////////////////////////////////////////////////////////////////////////
509     void cNoradSDP4::DeepCalcIntegrator(double *pxndot, double *pxnddt,
510     double *pxldot, const double &delt)
511     {
512     DeepCalcDotTerms(pxndot, pxnddt, pxldot);
513    
514     dp_xli = dp_xli + (*pxldot) * delt + (*pxndot) * dp_step2;
515     dp_xni = dp_xni + (*pxndot) * delt + (*pxnddt) * dp_step2;
516     dp_atime = dp_atime + delt;
517     }
518    
519     //////////////////////////////////////////////////////////////////////////////
520     bool cNoradSDP4::DeepSecular(double *xmdf, double *omgadf, double *xnode,
521     double *emm, double *xincc, double *xnn,
522     double *tsince)
523     {
524     xll = *xmdf;
525     omgasm = *omgadf;
526     xnodes = *xnode;
527     xn = *xnn;
528     t = *tsince;
529    
530     // Deep space secular effects
531     xll = xll + dp_ssl * t;
532     omgasm = omgasm + dp_ssg * t;
533     xnodes = xnodes + dp_ssh * t;
534     _em = m_Orbit.Eccentricity() + dp_sse * t;
535     xinc = m_Orbit.Inclination() + dp_ssi * t;
536    
537     if (xinc < 0.0)
538     {
539     xinc = -xinc;
540     xnodes = xnodes + PI;
541     omgasm = omgasm - PI;
542     }
543    
544     double xnddt = 0.0;
545     double xndot = 0.0;
546     double xldot = 0.0;
547     double ft = 0.0;
548     double delt = 0.0;
549    
550     bool fDone = false;
551    
552     if (dp_iresfl)
553     {
554     while (!fDone)
555     {
556     if ((dp_atime == 0.0) ||
557     ((t >= 0.0) && (dp_atime < 0.0)) ||
558     ((t < 0.0) && (dp_atime >= 0.0)))
559     {
560     if (t < 0)
561     delt = dp_stepn;
562     else
563     delt = dp_stepp;
564    
565     // Epoch restart
566     dp_atime = 0.0;
567     dp_xni = m_xnodp;
568     dp_xli = dp_xlamo;
569    
570     fDone = true;
571     }
572     else
573     {
574     if (fabs(t) < fabs(dp_atime))
575     {
576     delt = dp_stepp;
577    
578     if (t >= 0.0)
579     delt = dp_stepn;
580    
581     DeepCalcIntegrator(&xndot, &xnddt, &xldot, delt);
582     }
583     else
584     {
585     delt = dp_stepn;
586    
587     if (t > 0.0)
588     delt = dp_stepp;
589    
590     fDone = true;
591     }
592     }
593     }
594    
595     while (fabs(t - dp_atime) >= dp_stepp)
596     {
597     DeepCalcIntegrator(&xndot, &xnddt, &xldot, delt);
598     }
599    
600     ft = t - dp_atime;
601    
602     DeepCalcDotTerms(&xndot, &xnddt, &xldot);
603    
604     xn = dp_xni + xndot * ft + xnddt * ft * ft * 0.5;
605    
606     double xl = dp_xli + xldot * ft + xndot * ft * ft * 0.5;
607     double temp = -xnodes + dp_thgr + t * thdt;
608    
609     xll = xl - omgasm + temp;
610    
611     if (!dp_isynfl)
612     xll = xl + temp + temp;
613     }
614    
615     *xmdf = xll;
616     *omgadf = omgasm;
617     *xnode = xnodes;
618     *emm = _em;
619     *xincc = xinc;
620     *xnn = xn;
621     *tsince = t;
622    
623     return true;
624     }
625    
626     //////////////////////////////////////////////////////////////////////////////
627     bool cNoradSDP4::DeepPeriodics(double *e, double *xincc,
628     double *omgadf, double *xnode,
629     double *xmam)
630     {
631     _em = *e;
632     xinc = *xincc;
633     omgasm = *omgadf;
634     xnodes = *xnode;
635     xll = *xmam;
636    
637     // Lunar-solar periodics
638     double sinis = sin(xinc);
639     double cosis = cos(xinc);
640    
641     double sghs = 0.0;
642     double shs = 0.0;
643     double sh1 = 0.0;
644     double pe = 0.0;
645     double pinc = 0.0;
646     double pl = 0.0;
647     double sghl = 0.0;
648    
649     if (fabs(dp_savtsn - t) >= 30.0)
650     {
651     dp_savtsn = t;
652    
653     double zm = dp_zmos + zns * t;
654     double zf = zm + 2.0 * zes * sin(zm);
655     double sinzf = sin(zf);
656     double f2 = 0.5 * sinzf * sinzf - 0.25;
657     double f3 = -0.5 * sinzf * cos(zf);
658     double ses = dp_se2 * f2 + dp_se3 * f3;
659     double sis = dp_si2 * f2 + dp_si3 * f3;
660     double sls = dp_sl2 * f2 + dp_sl3 * f3 + dp_sl4 * sinzf;
661    
662     sghs = dp_sgh2 * f2 + dp_sgh3 * f3 + dp_sgh4 * sinzf;
663     shs = dp_sh2 * f2 + dp_sh3 * f3;
664     zm = dp_zmol + znl * t;
665     zf = zm + 2.0 * zel * sin(zm);
666     sinzf = sin(zf);
667     f2 = 0.5 * sinzf * sinzf - 0.25;
668     f3 = -0.5 * sinzf * cos(zf);
669    
670     double sel = dp_ee2 * f2 + dp_e3 * f3;
671     double sil = dp_xi2 * f2 + dp_xi3 * f3;
672     double sll = dp_xl2 * f2 + dp_xl3 * f3 + dp_xl4 * sinzf;
673    
674     sghl = dp_xgh2 * f2 + dp_xgh3 * f3 + dp_xgh4 * sinzf;
675     sh1 = dp_xh2 * f2 + dp_xh3 * f3;
676     pe = ses + sel;
677     pinc = sis + sil;
678     pl = sls + sll;
679     }
680    
681     double pgh = sghs + sghl;
682     double ph = shs + sh1;
683     xinc = xinc + pinc;
684     _em = _em + pe;
685    
686     if (dp_xqncl >= 0.2)
687     {
688     // Apply periodics directly
689     ph = ph / siniq;
690     pgh = pgh - cosiq * ph;
691     omgasm = omgasm + pgh;
692     xnodes = xnodes + ph;
693     xll = xll + pl;
694     }
695     else
696     {
697     // Apply periodics with Lyddane modification
698     double sinok = sin(xnodes);
699     double cosok = cos(xnodes);
700     double alfdp = sinis * sinok;
701     double betdp = sinis * cosok;
702     double dalf = ph * cosok + pinc * cosis * sinok;
703     double dbet = -ph * sinok + pinc * cosis * cosok;
704    
705     alfdp = alfdp + dalf;
706     betdp = betdp + dbet;
707    
708     double xls = xll + omgasm + cosis * xnodes;
709     double dls = pl + pgh - pinc * xnodes * sinis;
710    
711     xls = xls + dls;
712     xnodes = AcTan(alfdp, betdp);
713     xll = xll + pl;
714     omgasm = xls - xll - cos(xinc) * xnodes;
715     }
716    
717     *e = _em;
718     *xincc = xinc;
719     *omgadf = omgasm;
720     *xnode = xnodes;
721     *xmam = xll;
722    
723     return true;
724     }
725    
726     //////////////////////////////////////////////////////////////////////////////
727     // getPosition()
728     // This procedure returns the ECI position and velocity for the satellite
729     // in the orbit at the given number of minutes since the TLE epoch time
730     // using the NORAD Simplified General Perturbation 4, "deep space" orbit
731     // model.
732     //
733     // tsince - Time in minutes since the TLE epoch (GMT).
734     // pECI - pointer to location to store the ECI data.
735     // To convert the returned ECI position vector to km,
736     // multiply each component by:
737     // (XKMPER_WGS72 / AE).
738     // To convert the returned ECI velocity vector to km/sec,
739     // multiply each component by:
740     // (XKMPER_WGS72 / AE) * (MIN_PER_DAY / 86400).
741     bool cNoradSDP4::getPosition(double tsince, cEci &eci)
742     {
743     DeepInit(&m_eosq, &m_sinio, &m_cosio, &m_betao, &m_aodp, &m_theta2,
744     &m_sing, &m_cosg, &m_betao2, &m_xmdot, &m_omgdot, &m_xnodot);
745    
746     // Update for secular gravity and atmospheric drag
747     double xmdf = m_Orbit.mnAnomaly() + m_xmdot * tsince;
748     double omgadf = m_Orbit.ArgPerigee() + m_omgdot * tsince;
749     double xnoddf = m_Orbit.RAAN() + m_xnodot * tsince;
750     double tsq = tsince * tsince;
751     double xnode = xnoddf + m_xnodcf * tsq;
752     double tempa = 1.0 - m_c1 * tsince;
753     double tempe = m_Orbit.BStar() * m_c4 * tsince;
754     double templ = m_t2cof * tsq;
755     double xn = m_xnodp;
756     double em;
757     double xinc;
758    
759     DeepSecular(&xmdf, &omgadf, &xnode, &em, &xinc, &xn, &tsince);
760    
761     double a = pow(XKE / xn, TWOTHRD) * sqr(tempa);
762     double e = em - tempe;
763     double xmam = xmdf + m_xnodp * templ;
764    
765     DeepPeriodics(&e, &xinc, &omgadf, &xnode, &xmam);
766    
767     double xl = xmam + omgadf + xnode;
768    
769     xn = XKE / pow(a, 1.5);
770    
771     return FinalPosition(xinc, omgadf, e, a, xl, xnode, xn, tsince, eci);
772     }

  ViewVC Help
Powered by ViewVC 1.1.23