1   /* Copyright 2002-2024 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation.analytical.tle;
18  
19  import org.hipparchus.util.FastMath;
20  import org.hipparchus.util.MathUtils;
21  import org.hipparchus.util.SinCos;
22  import org.orekit.annotation.DefaultDataContext;
23  import org.orekit.attitudes.AttitudeProvider;
24  import org.orekit.data.DataContext;
25  import org.orekit.frames.Frame;
26  import org.orekit.time.DateTimeComponents;
27  import org.orekit.utils.Constants;
28  
29  
30  /** This class contains the methods that compute deep space perturbation terms.
31   * <p>
32   * The user should not bother in this class since it is handled internaly by the
33   * {@link TLEPropagator}.
34   * </p>
35   * <p>This implementation is largely inspired from the paper and source code <a
36   * href="https://www.celestrak.com/publications/AIAA/2006-6753/">Revisiting Spacetrack
37   * Report #3</a> and is fully compliant with its results and tests cases.</p>
38   * @author Felix R. Hoots, Ronald L. Roehrich, December 1980 (original fortran)
39   * @author David A. Vallado, Paul Crawford, Richard Hujsak, T.S. Kelso (C++ translation and improvements)
40   * @author Fabien Maussion (java translation)
41   */
42  public class DeepSDP4 extends SDP4 {
43  
44      // CHECKSTYLE: stop JavadocVariable check
45  
46      /** Integration step (seconds). */
47      private static final double SECULAR_INTEGRATION_STEP  = 720.0;
48  
49      /** Intermediate values. */
50      private double thgr;
51      private double xnq;
52      private double omegaq;
53      private double zcosil;
54      private double zsinil;
55      private double zsinhl;
56      private double zcoshl;
57      private double zmol;
58      private double zcosgl;
59      private double zsingl;
60      private double zmos;
61      private double savtsn;
62  
63      private double ee2;
64      private double e3;
65      private double xi2;
66      private double xi3;
67      private double xl2;
68      private double xl3;
69      private double xl4;
70      private double xgh2;
71      private double xgh3;
72      private double xgh4;
73      private double xh2;
74      private double xh3;
75  
76      private double d2201;
77      private double d2211;
78      private double d3210;
79      private double d3222;
80      private double d4410;
81      private double d4422;
82      private double d5220;
83      private double d5232;
84      private double d5421;
85      private double d5433;
86      private double xlamo;
87  
88      private double sse;
89      private double ssi;
90      private double ssl;
91      private double ssh;
92      private double ssg;
93      private double se2;
94      private double si2;
95      private double sl2;
96      private double sgh2;
97      private double sh2;
98      private double se3;
99      private double si3;
100     private double sl3;
101     private double sgh3;
102     private double sh3;
103     private double sl4;
104     private double sgh4;
105 
106     private double del1;
107     private double del2;
108     private double del3;
109     private double xfact;
110     private double xli;
111     private double xni;
112     private double atime;
113 
114     private double pe;
115     private double pinc;
116     private double pl;
117     private double pgh;
118     private double ph;
119 
120     private double[] derivs;
121 
122     // CHECKSTYLE: resume JavadocVariable check
123 
124     /** Flag for resonant orbits. */
125     private boolean resonant;
126 
127     /** Flag for synchronous orbits. */
128     private boolean synchronous;
129 
130     /** Flag for compliance with Dundee modifications. */
131     private boolean isDundeeCompliant = true;
132 
133     /** Constructor for a unique initial TLE.
134      *
135      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
136      *
137      * @param initialTLE the TLE to propagate.
138      * @param attitudeProvider provider for attitude computation
139      * @param mass spacecraft mass (kg)
140      * @see #DeepSDP4(TLE, AttitudeProvider, double, Frame)
141      */
142     @DefaultDataContext
143     public DeepSDP4(final TLE initialTLE, final AttitudeProvider attitudeProvider,
144                     final double mass) {
145         this(initialTLE, attitudeProvider, mass,
146                 DataContext.getDefault().getFrames().getTEME());
147     }
148 
149     /** Constructor for a unique initial TLE.
150      * @param initialTLE the TLE to propagate.
151      * @param attitudeProvider provider for attitude computation
152      * @param mass spacecraft mass (kg)
153      * @param teme the TEME frame to use for propagation.
154      * @since 10.1
155      */
156     public DeepSDP4(final TLE initialTLE,
157                     final AttitudeProvider attitudeProvider,
158                     final double mass,
159                     final Frame teme) {
160         super(initialTLE, attitudeProvider, mass, teme);
161     }
162 
163     /** Computes luni - solar terms from initial coordinates and epoch.
164      */
165     protected void luniSolarTermsComputation() {
166 
167         final SinCos scg  = FastMath.sinCos(tle.getPerigeeArgument());
168         final double sing = scg.sin();
169         final double cosg = scg.cos();
170 
171         final SinCos scq  = FastMath.sinCos(tle.getRaan());
172         final double sinq = scq.sin();
173         final double cosq = scq.cos();
174         final double aqnv = 1.0 / a0dp;
175 
176         // Compute julian days since 1900
177         final double daysSince1900 = (tle.getDate()
178                 .getComponents(utc)
179                 .offsetFrom(DateTimeComponents.JULIAN_EPOCH)) /
180                 Constants.JULIAN_DAY - 2415020;
181 
182         double cc = TLEConstants.C1SS;
183         double ze = TLEConstants.ZES;
184         double zn = TLEConstants.ZNS;
185         double zsinh = sinq;
186         double zcosh = cosq;
187 
188         thgr = thetaG(tle.getDate());
189         xnq = xn0dp;
190         omegaq = tle.getPerigeeArgument();
191 
192         final double xnodce = 4.5236020 - 9.2422029e-4 * daysSince1900;
193         final SinCos scTem  = FastMath.sinCos(xnodce);
194         final double stem   = scTem.sin();
195         final double ctem   = scTem.cos();
196         final double c_minus_gam = 0.228027132 * daysSince1900 - 1.1151842;
197         final double gam = 5.8351514 + 0.0019443680 * daysSince1900;
198 
199         zcosil = 0.91375164 - 0.03568096 * ctem;
200         zsinil = FastMath.sqrt(1.0 - zcosil * zcosil);
201         zsinhl = 0.089683511 * stem / zsinil;
202         zcoshl = FastMath.sqrt(1.0 - zsinhl * zsinhl);
203         zmol = MathUtils.normalizeAngle(c_minus_gam, FastMath.PI);
204 
205         double zx = 0.39785416 * stem / zsinil;
206         final double zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
207         zx = FastMath.atan2( zx, zy) + gam - xnodce;
208         final SinCos scZx = FastMath.sinCos(zx);
209         zcosgl = scZx.cos();
210         zsingl = scZx.sin();
211         zmos = MathUtils.normalizeAngle(6.2565837 + 0.017201977 * daysSince1900, FastMath.PI);
212 
213         // Do solar terms
214         savtsn = 1e20;
215 
216         double zcosi =  0.91744867;
217         double zsini =  0.39785416;
218         double zsing = -0.98088458;
219         double zcosg =  0.1945905;
220 
221         double se = 0;
222         double sgh = 0;
223         double sh = 0;
224         double si = 0;
225         double sl = 0;
226 
227         // There was previously some convoluted logic here, but it boils
228         // down to this:  we compute the solar terms,  then the lunar terms.
229         // On a second pass,  we recompute the solar terms, taking advantage
230         // of the improved data that resulted from computing lunar terms.
231         for (int iteration = 0; iteration < 2; ++iteration) {
232             final double a1 = zcosg * zcosh + zsing * zcosi * zsinh;
233             final double a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
234             final double a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
235             final double a8 = zsing * zsini;
236             final double a9 = zsing * zsinh + zcosg * zcosi * zcosh;
237             final double a10 = zcosg * zsini;
238             final double a2 = cosi0 * a7 + sini0 * a8;
239             final double a4 = cosi0 * a9 + sini0 * a10;
240             final double a5 = -sini0 * a7 + cosi0 * a8;
241             final double a6 = -sini0 * a9 + cosi0 * a10;
242             final double x1 = a1 * cosg + a2 * sing;
243             final double x2 = a3 * cosg + a4 * sing;
244             final double x3 = -a1 * sing + a2 * cosg;
245             final double x4 = -a3 * sing + a4 * cosg;
246             final double x5 = a5 * sing;
247             final double x6 = a6 * sing;
248             final double x7 = a5 * cosg;
249             final double x8 = a6 * cosg;
250             final double z31 = 12 * x1 * x1 - 3 * x3 * x3;
251             final double z32 = 24 * x1 * x2 - 6 * x3 * x4;
252             final double z33 = 12 * x2 * x2 - 3 * x4 * x4;
253             final double z11 = -6 * a1 * a5 + e0sq * (-24 * x1 * x7 - 6 * x3 * x5);
254             final double z12 = -6 * (a1 * a6 + a3 * a5) +
255                                e0sq * (-24 * (x2 * x7 + x1 * x8) - 6 * (x3 * x6 + x4 * x5));
256             final double z13 = -6 * a3 * a6 + e0sq * (-24 * x2 * x8 - 6 * x4 * x6);
257             final double z21 = 6 * a2 * a5 + e0sq * (24 * x1 * x5 - 6 * x3 * x7);
258             final double z22 = 6 * (a4 * a5 + a2 * a6) +
259                                e0sq * (24 * (x2 * x5 + x1 * x6) - 6 * (x4 * x7 + x3 * x8));
260             final double z23 = 6 * a4 * a6 + e0sq * (24 * x2 * x6 - 6 * x4 * x8);
261             final double s3 = cc / xnq;
262             final double s2 = -0.5 * s3 / beta0;
263             final double s4 = s3 * beta0;
264             final double s1 = -15 * tle.getE() * s4;
265             final double s5 = x1 * x3 + x2 * x4;
266             final double s6 = x2 * x3 + x1 * x4;
267             final double s7 = x2 * x4 - x1 * x3;
268             double z1 = 3 * (a1 * a1 + a2 * a2) + z31 * e0sq;
269             double z2 = 6 * (a1 * a3 + a2 * a4) + z32 * e0sq;
270             double z3 = 3 * (a3 * a3 + a4 * a4) + z33 * e0sq;
271 
272             z1 = z1 + z1 + beta02 * z31;
273             z2 = z2 + z2 + beta02 * z32;
274             z3 = z3 + z3 + beta02 * z33;
275             se = s1 * zn * s5;
276             si = s2 * zn * (z11 + z13);
277             sl = -zn * s3 * (z1 + z3 - 14 - 6 * e0sq);
278             sgh = s4 * zn * (z31 + z33 - 6);
279             if (tle.getI() < (FastMath.PI / 60.0)) {
280                 // inclination smaller than 3 degrees
281                 sh = 0;
282             } else {
283                 sh = -zn * s2 * (z21 + z23);
284             }
285             ee2  =  2 * s1 * s6;
286             e3   =  2 * s1 * s7;
287             xi2  =  2 * s2 * z12;
288             xi3  =  2 * s2 * (z13 - z11);
289             xl2  = -2 * s3 * z2;
290             xl3  = -2 * s3 * (z3 - z1);
291             xl4  = -2 * s3 * (-21 - 9 * e0sq) * ze;
292             xgh2 =  2 * s4 * z32;
293             xgh3 =  2 * s4 * (z33 - z31);
294             xgh4 = -18 * s4 * ze;
295             xh2  = -2 * s2 * z22;
296             xh3  = -2 * s2 * (z23 - z21);
297 
298             if (iteration == 0) { // we compute lunar terms only on the first pass:
299                 sse = se;
300                 ssi = si;
301                 ssl = sl;
302                 ssh = (tle.getI() < (FastMath.PI / 60.0)) ? 0 : sh / sini0;
303                 ssg = sgh - cosi0 * ssh;
304                 se2 = ee2;
305                 si2 = xi2;
306                 sl2 = xl2;
307                 sgh2 = xgh2;
308                 sh2 = xh2;
309                 se3 = e3;
310                 si3 = xi3;
311                 sl3 = xl3;
312                 sgh3 = xgh3;
313                 sh3 = xh3;
314                 sl4 = xl4;
315                 sgh4 = xgh4;
316                 zcosg = zcosgl;
317                 zsing = zsingl;
318                 zcosi = zcosil;
319                 zsini = zsinil;
320                 zcosh = zcoshl * cosq + zsinhl * sinq;
321                 zsinh = sinq * zcoshl - cosq * zsinhl;
322                 zn = TLEConstants.ZNL;
323                 cc = TLEConstants.C1L;
324                 ze = TLEConstants.ZEL;
325             }
326         } // end of solar - lunar - solar terms computation
327 
328         sse += se;
329         ssi += si;
330         ssl += sl;
331         ssg += sgh - ((tle.getI() < (FastMath.PI / 60.0)) ? 0 : (cosi0 / sini0 * sh));
332         ssh += (tle.getI() < (FastMath.PI / 60.0)) ? 0 : sh / sini0;
333 
334 
335 
336         //        Start the resonant-synchronous tests and initialization
337 
338         double bfact = 0;
339 
340         // if mean motion is 1.893053 to 2.117652 revs/day, and eccentricity >= 0.5,
341         // start of the 12-hour orbit, e > 0.5 section
342         if (xnq >= 0.00826 && xnq <= 0.00924 && tle.getE() >= 0.5) {
343 
344             final double g201 = -0.306 - (tle.getE() - 0.64) * 0.440;
345             final double eoc = tle.getE() * e0sq;
346             final double sini2 = sini0 * sini0;
347             final double f220 = 0.75 * (1 + 2 * cosi0 + theta2);
348             final double f221 = 1.5 * sini2;
349             final double f321 =  1.875 * sini0 * (1 - 2 * cosi0 - 3 * theta2);
350             final double f322 = -1.875 * sini0 * (1 + 2 * cosi0 - 3 * theta2);
351             final double f441 = 35 * sini2 * f220;
352             final double f442 = 39.3750 * sini2 * sini2;
353             final double f522 = 9.84375 * sini0 * (sini2 * (1 - 2 * cosi0 - 5 * theta2) +
354                                                    0.33333333 * (-2 + 4 * cosi0 + 6 * theta2));
355             final double f523 = sini0 * (4.92187512 * sini2 * (-2 - 4 * cosi0 + 10 * theta2) +
356                                          6.56250012 * (1 + 2 * cosi0 - 3 * theta2));
357             final double f542 = 29.53125 * sini0 * (2 - 8 * cosi0 + theta2 * (-12 + 8 * cosi0 + 10 * theta2));
358             final double f543 = 29.53125 * sini0 * (-2 - 8 * cosi0 + theta2 * (12 + 8 * cosi0 - 10 * theta2));
359             final double g211;
360             final double g310;
361             final double g322;
362             final double g410;
363             final double g422;
364             final double g520;
365 
366             resonant = true;       // it is resonant...
367             synchronous = false;     // but it's not synchronous
368 
369             // Geopotential resonance initialization for 12 hour orbits :
370             if (tle.getE() <= 0.65) {
371                 g211 =    3.616  -   13.247  * tle.getE() +   16.290  * e0sq;
372                 g310 =  -19.302  +  117.390  * tle.getE() -  228.419  * e0sq +  156.591  * eoc;
373                 g322 =  -18.9068 +  109.7927 * tle.getE() -  214.6334 * e0sq +  146.5816 * eoc;
374                 g410 =  -41.122  +  242.694  * tle.getE() -  471.094  * e0sq +  313.953  * eoc;
375                 g422 = -146.407  +  841.880  * tle.getE() - 1629.014  * e0sq + 1083.435  * eoc;
376                 g520 = -532.114  + 3017.977  * tle.getE() - 5740.032  * e0sq + 3708.276  * eoc;
377             } else  {
378                 g211 =   -72.099 +   331.819 * tle.getE() -   508.738 * e0sq +   266.724 * eoc;
379                 g310 =  -346.844 +  1582.851 * tle.getE() -  2415.925 * e0sq +  1246.113 * eoc;
380                 g322 =  -342.585 +  1554.908 * tle.getE() -  2366.899 * e0sq +  1215.972 * eoc;
381                 g410 = -1052.797 +  4758.686 * tle.getE() -  7193.992 * e0sq +  3651.957 * eoc;
382                 g422 = -3581.69  + 16178.11  * tle.getE() - 24462.77  * e0sq + 12422.52  * eoc;
383                 if (tle.getE() <= 0.715) {
384                     g520 = 1464.74 - 4664.75 * tle.getE() + 3763.64 * e0sq;
385                 } else {
386                     g520 = -5149.66 + 29936.92 * tle.getE() - 54087.36 * e0sq + 31324.56 * eoc;
387                 }
388             }
389 
390             final double g533;
391             final double g521;
392             final double g532;
393             if (tle.getE() < 0.7) {
394                 g533 = -919.2277  + 4988.61   * tle.getE() - 9064.77   * e0sq + 5542.21  * eoc;
395                 g521 = -822.71072 + 4568.6173 * tle.getE() - 8491.4146 * e0sq + 5337.524 * eoc;
396                 g532 = -853.666   + 4690.25   * tle.getE() - 8624.77   * e0sq + 5341.4   * eoc;
397             } else {
398                 g533 = -37995.78  + 161616.52 * tle.getE() - 229838.2  * e0sq + 109377.94 * eoc;
399                 g521 = -51752.104 + 218913.95 * tle.getE() - 309468.16 * e0sq + 146349.42 * eoc;
400                 g532 = -40023.88  + 170470.89 * tle.getE() - 242699.48 * e0sq + 115605.82 * eoc;
401             }
402 
403             double temp1 = 3 * xnq * xnq * aqnv * aqnv;
404             double temp = temp1 * TLEConstants.ROOT22;
405             d2201 = temp * f220 * g201;
406             d2211 = temp * f221 * g211;
407             temp1 *= aqnv;
408             temp = temp1 * TLEConstants.ROOT32;
409             d3210 = temp * f321 * g310;
410             d3222 = temp * f322 * g322;
411             temp1 *= aqnv;
412             temp = 2 * temp1 * TLEConstants.ROOT44;
413             d4410 = temp * f441 * g410;
414             d4422 = temp * f442 * g422;
415             temp1 *= aqnv;
416             temp = temp1 * TLEConstants.ROOT52;
417             d5220 = temp * f522 * g520;
418             d5232 = temp * f523 * g532;
419             temp = 2 * temp1 * TLEConstants.ROOT54;
420             d5421 = temp * f542 * g521;
421             d5433 = temp * f543 * g533;
422             xlamo = tle.getMeanAnomaly() + tle.getRaan() + tle.getRaan() - thgr - thgr;
423             bfact = xmdot + xnodot + xnodot - TLEConstants.THDT - TLEConstants.THDT;
424             bfact += ssl + ssh + ssh;
425         } else if (xnq < 0.0052359877 && xnq > 0.0034906585) {
426             // if mean motion is .8 to 1.2 revs/day : (geosynch)
427 
428             final double cosio_plus_1 = 1.0 + cosi0;
429             final double g200 = 1 + e0sq * (-2.5 + 0.8125  * e0sq);
430             final double g300 = 1 + e0sq * (-6   + 6.60937 * e0sq);
431             final double f311 = 0.9375 * sini0 * sini0 * (1 + 3 * cosi0) - 0.75 * cosio_plus_1;
432             final double g310 = 1 + 2 * e0sq;
433             final double f220 = 0.75 * cosio_plus_1 * cosio_plus_1;
434             final double f330 = 2.5 * f220 * cosio_plus_1;
435 
436             resonant = true;
437             synchronous = true;
438 
439             // Synchronous resonance terms initialization
440             del1 = 3 * xnq * xnq * aqnv * aqnv;
441             del2 = 2 * del1 * f220 * g200 * TLEConstants.Q22;
442             del3 = 3 * del1 * f330 * g300 * TLEConstants.Q33 * aqnv;
443             del1 = del1 * f311 * g310 * TLEConstants.Q31 * aqnv;
444             xlamo = tle.getMeanAnomaly() + tle.getRaan() + tle.getPerigeeArgument() - thgr;
445             bfact = xmdot + omgdot + xnodot - TLEConstants.THDT;
446             bfact = bfact + ssl + ssg + ssh;
447         } else {
448             // it's neither a high-e 12-hours orbit nor a geosynchronous:
449             resonant = false;
450             synchronous = false;
451         }
452 
453         if (resonant) {
454             xfact = bfact - xnq;
455 
456             // Initialize integrator
457             xli   = xlamo;
458             xni   = xnq;
459             atime = 0;
460         }
461         derivs = new double[2];
462     }
463 
464     /** Computes secular terms from current coordinates and epoch.
465      * @param t offset from initial epoch (minutes)
466      */
467     protected void deepSecularEffects(final double t)  {
468 
469         xll    += ssl * t;
470         omgadf += ssg * t;
471         xnode  += ssh * t;
472         em      = tle.getE() + sse * t;
473         xinc    = tle.getI() + ssi * t;
474 
475         if (resonant) {
476             // If we're closer to t = 0 than to the currently-stored data
477             // from the previous call to this function,  then we're
478             // better off "restarting",  going back to the initial data.
479             // The Dundee code rigs things up to _always_ take 720-minute
480             // steps from epoch to end time,  except for the final step.
481             // Easiest way to arrange similar behavior in this code is
482             // just to always do a restart,  if we're in Dundee-compliant
483             // mode.
484             if (FastMath.abs(t) < FastMath.abs(t - atime) || isDundeeCompliant)  {
485                 // Epoch restart
486                 atime = 0;
487                 xni = xnq;
488                 xli = xlamo;
489             }
490             boolean lastIntegrationStep = false;
491             // if |step|>|step max| then do one step at step max
492             while (!lastIntegrationStep) {
493                 double delt = t - atime;
494                 if (delt > SECULAR_INTEGRATION_STEP) {
495                     delt = SECULAR_INTEGRATION_STEP;
496                 } else if (delt < -SECULAR_INTEGRATION_STEP) {
497                     delt = -SECULAR_INTEGRATION_STEP;
498                 } else {
499                     lastIntegrationStep = true;
500                 }
501 
502                 computeSecularDerivs();
503 
504                 final double xldot = xni + xfact;
505 
506                 double xlpow = 1.;
507                 xli += delt * xldot;
508                 xni += delt * derivs[0];
509                 double delt_factor = delt;
510                 xlpow *= xldot;
511                 derivs[1] *= xlpow;
512                 delt_factor *= delt / 2;
513                 xli += delt_factor * derivs[0];
514                 xni += delt_factor * derivs[1];
515                 atime += delt;
516             }
517             xn = xni;
518             final double temp = -xnode + thgr + t * TLEConstants.THDT;
519             xll = xli + temp + (synchronous ? -omgadf : temp);
520         }
521     }
522 
523     /** Computes periodic terms from current coordinates and epoch.
524      * @param t offset from initial epoch (min)
525      */
526     protected void deepPeriodicEffects(final double t)  {
527 
528         // If the time didn't change by more than 30 minutes,
529         // there's no good reason to recompute the perturbations;
530         // they don't change enough over so short a time span.
531         // However,  the Dundee code _always_ recomputes,  so if
532         // we're attempting to replicate its results,  we've gotta
533         // recompute everything,  too.
534         if (FastMath.abs(savtsn - t) >= 30.0 || isDundeeCompliant)  {
535 
536             savtsn = t;
537 
538             // Update solar perturbations for time T
539             double zm = zmos + TLEConstants.ZNS * t;
540             double zf = zm + 2 * TLEConstants.ZES * FastMath.sin(zm);
541             SinCos sczf = FastMath.sinCos(zf);
542             double sinzf = sczf.sin();
543             double f2 = 0.5 * sinzf * sinzf - 0.25;
544             double f3 = -0.5 * sinzf * sczf.cos();
545             final double ses = se2 * f2 + se3 * f3;
546             final double sis = si2 * f2 + si3 * f3;
547             final double sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
548             final double sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
549             final double shs = sh2 * f2 + sh3 * f3;
550 
551             // Update lunar perturbations for time T
552             zm = zmol + TLEConstants.ZNL * t;
553             zf = zm + 2 * TLEConstants.ZEL * FastMath.sin(zm);
554             sczf = FastMath.sinCos(zf);
555             sinzf = sczf.sin();
556             f2 =  0.5 * sinzf * sinzf - 0.25;
557             f3 = -0.5 * sinzf * sczf.cos();
558             final double sel = ee2 * f2 + e3 * f3;
559             final double sil = xi2 * f2 + xi3 * f3;
560             final double sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf;
561             final double sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
562             final double sh1 = xh2 * f2 + xh3 * f3;
563 
564             // Sum the solar and lunar contributions
565             pe   = ses  + sel;
566             pinc = sis  + sil;
567             pl   = sls  + sll;
568             pgh  = sghs + sghl;
569             ph   = shs  + sh1;
570         }
571 
572         xinc += pinc;
573 
574         final SinCos scis = FastMath.sinCos(xinc);
575         final double sinis = scis.sin();
576         final double cosis = scis.cos();
577 
578         /* Add solar/lunar perturbation correction to eccentricity: */
579         em     += pe;
580         xll    += pl;
581         omgadf += pgh;
582         xinc    = MathUtils.normalizeAngle(xinc, 0);
583 
584         if (FastMath.abs(xinc) >= 0.2) {
585             // Apply periodics directly
586             final double temp_val = ph / sinis;
587             omgadf -= cosis * temp_val;
588             xnode += temp_val;
589         } else {
590             // Apply periodics with Lyddane modification
591             final SinCos scok  = FastMath.sinCos(xnode);
592             final double sinok = scok.sin();
593             final double cosok = scok.cos();
594             final double alfdp =  ph * cosok + (pinc * cosis + sinis) * sinok;
595             final double betdp = -ph * sinok + (pinc * cosis + sinis) * cosok;
596             final double delta_xnode = MathUtils.normalizeAngle(FastMath.atan2(alfdp, betdp) - xnode, 0);
597             final double dls = -xnode * sinis * pinc;
598             omgadf += dls - cosis * delta_xnode;
599             xnode  += delta_xnode;
600         }
601     }
602 
603     /** Computes internal secular derivs. */
604     private void computeSecularDerivs() {
605 
606         final SinCos sc_li  = FastMath.sinCos(xli);
607         final double sin_li = sc_li.sin();
608         final double cos_li = sc_li.cos();
609         final double sin_2li = 2. * sin_li * cos_li;
610         final double cos_2li = 2. * cos_li * cos_li - 1.;
611 
612         // Dot terms calculated :
613         if (synchronous)  {
614             final double sin_3li = sin_2li * cos_li + cos_2li * sin_li;
615             final double cos_3li = cos_2li * cos_li - sin_2li * sin_li;
616             final double term1a = del1 * (sin_li  * TLEConstants.C_FASX2  - cos_li  * TLEConstants.S_FASX2);
617             final double term2a = del2 * (sin_2li * TLEConstants.C_2FASX4 - cos_2li * TLEConstants.S_2FASX4);
618             final double term3a = del3 * (sin_3li * TLEConstants.C_3FASX6 - cos_3li * TLEConstants.S_3FASX6);
619             final double term1b = del1 * (cos_li  * TLEConstants.C_FASX2  + sin_li  * TLEConstants.S_FASX2);
620             final double term2b = 2.0 * del2 * (cos_2li * TLEConstants.C_2FASX4 + sin_2li * TLEConstants.S_2FASX4);
621             final double term3b = 3.0 * del3 * (cos_3li * TLEConstants.C_3FASX6 + sin_3li * TLEConstants.S_3FASX6);
622             derivs[0] = term1a + term2a + term3a;
623             derivs[1] = term1b + term2b + term3b;
624         } else {
625             // orbit is a 12-hour resonant one
626             final double xomi = omegaq + omgdot * atime;
627             final SinCos sc_omi  = FastMath.sinCos(xomi);
628             final double sin_omi = sc_omi.sin();
629             final double cos_omi = sc_omi.cos();
630             final double sin_li_m_omi = sin_li * cos_omi - sin_omi * cos_li;
631             final double sin_li_p_omi = sin_li * cos_omi + sin_omi * cos_li;
632             final double cos_li_m_omi = cos_li * cos_omi + sin_omi * sin_li;
633             final double cos_li_p_omi = cos_li * cos_omi - sin_omi * sin_li;
634             final double sin_2omi = 2. * sin_omi * cos_omi;
635             final double cos_2omi = 2. * cos_omi * cos_omi - 1.;
636             final double sin_2li_m_omi = sin_2li * cos_omi - sin_omi * cos_2li;
637             final double sin_2li_p_omi = sin_2li * cos_omi + sin_omi * cos_2li;
638             final double cos_2li_m_omi = cos_2li * cos_omi + sin_omi * sin_2li;
639             final double cos_2li_p_omi = cos_2li * cos_omi - sin_omi * sin_2li;
640             final double sin_2li_p_2omi = sin_2li * cos_2omi + sin_2omi * cos_2li;
641             final double cos_2li_p_2omi = cos_2li * cos_2omi - sin_2omi * sin_2li;
642             final double sin_2omi_p_li = sin_li * cos_2omi + sin_2omi * cos_li;
643             final double cos_2omi_p_li = cos_li * cos_2omi - sin_2omi * sin_li;
644             final double term1a = d2201 * (sin_2omi_p_li  * TLEConstants.C_G22 - cos_2omi_p_li  * TLEConstants.S_G22) +
645                                   d2211 * (sin_li         * TLEConstants.C_G22 - cos_li         * TLEConstants.S_G22) +
646                                   d3210 * (sin_li_p_omi   * TLEConstants.C_G32 - cos_li_p_omi   * TLEConstants.S_G32) +
647                                   d3222 * (sin_li_m_omi   * TLEConstants.C_G32 - cos_li_m_omi   * TLEConstants.S_G32) +
648                                   d5220 * (sin_li_p_omi   * TLEConstants.C_G52 - cos_li_p_omi   * TLEConstants.S_G52) +
649                                   d5232 * (sin_li_m_omi   * TLEConstants.C_G52 - cos_li_m_omi   * TLEConstants.S_G52);
650             final double term2a = d4410 * (sin_2li_p_2omi * TLEConstants.C_G44 - cos_2li_p_2omi * TLEConstants.S_G44) +
651                                   d4422 * (sin_2li        * TLEConstants.C_G44 - cos_2li        * TLEConstants.S_G44) +
652                                   d5421 * (sin_2li_p_omi  * TLEConstants.C_G54 - cos_2li_p_omi  * TLEConstants.S_G54) +
653                                   d5433 * (sin_2li_m_omi  * TLEConstants.C_G54 - cos_2li_m_omi  * TLEConstants.S_G54);
654             final double term1b = d2201 * (cos_2omi_p_li  * TLEConstants.C_G22 + sin_2omi_p_li  * TLEConstants.S_G22) +
655                                   d2211 * (cos_li         * TLEConstants.C_G22 + sin_li         * TLEConstants.S_G22) +
656                                   d3210 * (cos_li_p_omi   * TLEConstants.C_G32 + sin_li_p_omi   * TLEConstants.S_G32) +
657                                   d3222 * (cos_li_m_omi   * TLEConstants.C_G32 + sin_li_m_omi   * TLEConstants.S_G32) +
658                                   d5220 * (cos_li_p_omi   * TLEConstants.C_G52 + sin_li_p_omi   * TLEConstants.S_G52) +
659                                   d5232 * (cos_li_m_omi   * TLEConstants.C_G52 + sin_li_m_omi   * TLEConstants.S_G52);
660             final double term2b = 2.0 * (d4410 * (cos_2li_p_2omi * TLEConstants.C_G44 + sin_2li_p_2omi * TLEConstants.S_G44) +
661                                          d4422 * (cos_2li        * TLEConstants.C_G44 + sin_2li        * TLEConstants.S_G44) +
662                                          d5421 * (cos_2li_p_omi  * TLEConstants.C_G54 + sin_2li_p_omi  * TLEConstants.S_G54) +
663                                          d5433 * (cos_2li_m_omi  * TLEConstants.C_G54 + sin_2li_m_omi  * TLEConstants.S_G54));
664 
665             derivs[0] = term1a + term2a;
666             derivs[1] = term1b + term2b;
667 
668         }
669     }
670 
671 }