/* Copyright 2002-2025 CS GROUP
    * Licensed to CS GROUP (CS) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * CS licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *   http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.orekit.propagation.analytical.tle;
  
  import org.hipparchus.CalculusFieldElement;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.FieldSinCos;
  import org.hipparchus.util.MathArrays;
  import org.hipparchus.util.MathUtils;
  import org.hipparchus.util.SinCos;
  import org.orekit.annotation.DefaultDataContext;
  import org.orekit.attitudes.AttitudeProvider;
  import org.orekit.data.DataContext;
  import org.orekit.frames.Frame;
  import org.orekit.time.DateTimeComponents;
  import org.orekit.utils.Constants;
  
  
  /** This class contains the methods that compute deep space perturbation terms.
   * <p>
   * The user should not bother in this class since it is handled internaly by the
   * {@link TLEPropagator}.
   * </p>
   * <p>This implementation is largely inspired from the paper and source code <a
   * href="https://www.celestrak.com/publications/AIAA/2006-6753/">Revisiting Spacetrack
   * Report #3</a> and is fully compliant with its results and tests cases.</p>
   * @author Felix R. Hoots, Ronald L. Roehrich, December 1980 (original fortran)
   * @author David A. Vallado, Paul Crawford, Richard Hujsak, T.S. Kelso (C++ translation and improvements)
   * @author Fabien Maussion (java translation)
   * @author Thomas Paulet (field translation)
   * @since 11.0
   * @param <T> type of the field elements
   */
  public class FieldDeepSDP4<T extends CalculusFieldElement<T>> extends FieldSDP4<T> {
  
      // CHECKSTYLE: stop JavadocVariable check
  
      /** Integration step (seconds). */
      private static final double SECULAR_INTEGRATION_STEP  = 720.0;
  
      /** Intermediate values. */
      private double thgr;
      private T xnq;
      private T omegaq;
      private double zcosil;
      private double zsinil;
      private double zsinhl;
      private double zcoshl;
      private double zmol;
      private double zcosgl;
      private double zsingl;
      private double zmos;
      private T savtsn;
  
      private T ee2;
      private T e3;
      private T xi2;
      private T xi3;
      private T xl2;
      private T xl3;
      private T xl4;
      private T xgh2;
      private T xgh3;
      private T xgh4;
      private T xh2;
      private T xh3;
  
      private T d2201;
      private T d2211;
      private T d3210;
      private T d3222;
      private T d4410;
      private T d4422;
      private T d5220;
      private T d5232;
      private T d5421;
      private T d5433;
      private T xlamo;
  
      private T sse;
      private T ssi;
      private T ssl;
      private T ssh;
      private T ssg;
      private T se2;
     private T si2;
     private T sl2;
     private T sgh2;
     private T sh2;
     private T se3;
     private T si3;
     private T sl3;
     private T sgh3;
     private T sh3;
     private T sl4;
     private T sgh4;
 
     private T del1;
     private T del2;
     private T del3;
     private T xfact;
     private T xli;
     private T xni;
     private T atime;
 
     private T pe;
     private T pinc;
     private T pl;
     private T pgh;
     private T ph;
 
     private T[] derivs;
 
     // CHECKSTYLE: resume JavadocVariable check
 
     /** Flag for resonant orbits. */
     private boolean resonant;
 
     /** Flag for synchronous orbits. */
     private boolean synchronous;
 
     /** Flag for compliance with Dundee modifications. */
     private boolean isDundeeCompliant = true;
 
     /** Constructor for a unique initial TLE.
      *
      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
      *
      * @param initialTLE the TLE to propagate.
      * @param attitudeProvider provider for attitude computation
      * @param mass spacecraft mass (kg)
      * @param parameters SGP4 and SDP4 model parameters
      * @see #FieldDeepSDP4(FieldTLE, AttitudeProvider, CalculusFieldElement, Frame, CalculusFieldElement[])
      */
     @DefaultDataContext
     public FieldDeepSDP4(final FieldTLE<T> initialTLE, final AttitudeProvider attitudeProvider,
                     final T mass, final T[] parameters) {
         this(initialTLE, attitudeProvider, mass,
                 DataContext.getDefault().getFrames().getTEME(), parameters);
     }
 
     /** Constructor for a unique initial TLE.
      * @param initialTLE the TLE to propagate.
      * @param attitudeProvider provider for attitude computation
      * @param mass spacecraft mass (kg)
      * @param teme the TEME frame to use for propagation.
      * @param parameters SGP4 and SDP4 model parameters
      */
     public FieldDeepSDP4(final FieldTLE<T> initialTLE,
                          final AttitudeProvider attitudeProvider,
                          final T mass,
                          final Frame teme,
                          final T[] parameters) {
         super(initialTLE, attitudeProvider, mass, teme, parameters);
     }
 
     /** Computes luni - solar terms from initial coordinates and epoch.
      */
     protected void luniSolarTermsComputation() {
 
         final T zero = tle.getPerigeeArgument().getField().getZero();
         final T pi   = zero.getPi();
 
         final FieldSinCos<T> scg  = FastMath.sinCos(tle.getPerigeeArgument());
         final T sing = scg.sin();
         final T cosg = scg.cos();
 
         final FieldSinCos<T> scq  = FastMath.sinCos(tle.getRaan());
         final T sinq = scq.sin();
         final T cosq = scq.cos();
         final T aqnv = a0dp.reciprocal();
 
         // Compute julian days since 1900
         final double daysSince1900 = (tle.getDate()
                 .getComponents(utc)
                 .offsetFrom(DateTimeComponents.JULIAN_EPOCH)) /
                 Constants.JULIAN_DAY - 2415020;
 
         double cc = TLEConstants.C1SS;
         double ze = TLEConstants.ZES;
         double zn = TLEConstants.ZNS;
         T zsinh = sinq;
         T zcosh = cosq;
 
         thgr = thetaG(tle.getDate());
         xnq = xn0dp;
         omegaq = tle.getPerigeeArgument();
 
         final double xnodce = 4.5236020 - 9.2422029e-4 * daysSince1900;
         final SinCos scTem  = FastMath.sinCos(xnodce);
         final double stem = scTem.sin();
         final double ctem = scTem.cos();
         final double c_minus_gam = 0.228027132 * daysSince1900 - 1.1151842;
         final double gam = 5.8351514 + 0.0019443680 * daysSince1900;
 
         zcosil = 0.91375164 - 0.03568096 * ctem;
         zsinil = FastMath.sqrt(1.0 - zcosil * zcosil);
         zsinhl = 0.089683511 * stem / zsinil;
         zcoshl = FastMath.sqrt(1.0 - zsinhl * zsinhl);
         zmol = MathUtils.normalizeAngle(c_minus_gam, pi.getReal());
 
         double zx = 0.39785416 * stem / zsinil;
         final double zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
         zx = FastMath.atan2( zx, zy) + gam - xnodce;
         final SinCos scZx = FastMath.sinCos(zx);
         zcosgl = scZx.cos();
         zsingl = scZx.sin();
         zmos = MathUtils.normalizeAngle(6.2565837 + 0.017201977 * daysSince1900, pi.getReal());
 
         // Do solar terms
         savtsn = zero.newInstance(1e20);
 
         T zcosi = zero.newInstance(0.91744867);
         T zsini = zero.newInstance(0.39785416);
         T zsing = zero.newInstance(-0.98088458);
         T zcosg = zero.newInstance(0.1945905);
 
         T se =  zero;
         T sgh = zero;
         T sh =  zero;
         T si =  zero;
         T sl =  zero;
 
         // There was previously some convoluted logic here, but it boils
         // down to this:  we compute the solar terms,  then the lunar terms.
         // On a second pass,  we recompute the solar terms, taking advantage
         // of the improved data that resulted from computing lunar terms.
         for (int iteration = 0; iteration < 2; ++iteration) {
             final T a1  = zcosh.multiply(zcosg).add(zsinh.multiply(zsing).multiply(zcosi));
             final T a3  = zcosh.multiply(zsing.negate()).add(zsinh.multiply(zcosg).multiply(zcosi));
             final T a7  = zsinh.negate().multiply(zcosg).add(zcosh.multiply(zcosi).multiply(zsing));
             final T a8  = zsing.multiply(zsini);
             final T a9  = zsinh.multiply(zsing).add(zcosh.multiply(zcosi).multiply(zcosg));
             final T a10 = zcosg.multiply(zsini);
             final T a2  = cosi0.multiply(a7).add(sini0.multiply(a8));
             final T a4  = cosi0.multiply(a9).add(sini0.multiply(a10));
             final T a5  = sini0.negate().multiply(a7).add(cosi0.multiply(a8));
             final T a6  = sini0.negate().multiply(a9).add(cosi0.multiply(a10));
             final T x1  = a1.multiply(cosg).add(a2.multiply(sing));
             final T x2  = a3.multiply(cosg).add(a4.multiply(sing));
             final T x3  = a1.negate().multiply(sing).add(a2.multiply(cosg));
             final T x4  = a3.negate().multiply(sing).add(a4.multiply(cosg));
             final T x5  = a5.multiply(sing);
             final T x6  = a6.multiply(sing);
             final T x7  = a5.multiply(cosg);
             final T x8  = a6.multiply(cosg);
             final T z31 = x1.square().multiply(12).subtract(x3.square().multiply(3));
             final T z32 = x1.multiply(x2).multiply(24).subtract(x3.multiply(x4).multiply(6));
             final T z33 = x2.square().multiply(12).subtract(x4.square().multiply(3));
             final T z11 = a1.multiply(-6).multiply(a5).add(e0sq.multiply(x1.multiply(x7).multiply(-24).add(x3.multiply(x5).multiply(-6))));
             final T z12 = a1.multiply(a6).add(a3.multiply(a5)).multiply(-6).add(
                                 e0sq.multiply(x2.multiply(x7).add(x1.multiply(x8)).multiply(-24).add(
                                 x3.multiply(x6).add(x4.multiply(x5)).multiply(-6))));
             final T z13 = a3.multiply(a6).multiply(-6).add(e0sq.multiply(
                                x2.multiply(x8).multiply(-24).subtract(x4.multiply(x6).multiply(6))));
             final T z21 = a2.multiply(a5).multiply(6).add(e0sq.multiply(
                                x1.multiply(x5).multiply(24).subtract(x3.multiply(x7).multiply(6))));
             final T z22 = a4.multiply(a5).add(a2.multiply(a6)).multiply(6).add(
                                e0sq.multiply(x2.multiply(x5).add(x1.multiply(x6)).multiply(24).subtract(
                                x4.multiply(x7).add(x3.multiply(x8)).multiply(6))));
             final T z23 = a4.multiply(a6).multiply(6).add(e0sq.multiply(x2.multiply(x6).multiply(24).subtract(x4.multiply(x8).multiply(6))));
             final T s3  = xnq.reciprocal().multiply(cc);
             final T s2  = beta0.reciprocal().multiply(s3.multiply(-0.5));
             final T s4  = s3.multiply(beta0);
             final T s1  = tle.getE().multiply(s4).multiply(-15);
             final T s5  = x1.multiply(x3).add(x2.multiply(x4));
             final T s6  = x2.multiply(x3).add(x1.multiply(x4));
             final T s7  = x2.multiply(x4).subtract(x1.multiply(x3));
             T z1 = a1.square().add(a2.square()).multiply(3).add(z31.multiply(e0sq));
             T z2 = a1.multiply(a3).add(a2.multiply(a4)).multiply(6).add(z32.multiply(e0sq));
             T z3 = a3.square().add(a4.square()).multiply(3).add(z33.multiply(e0sq));
 
             z1 = z1.add(z1).add(beta02.multiply(z31));
             z2 = z2.add(z2).add(beta02.multiply(z32));
             z3 = z3.add(z3).add(beta02.multiply(z33));
             se = s1.multiply(zn).multiply(s5);
             si = s2.multiply(zn).multiply(z11.add(z13));
             sl = s3.multiply(-zn).multiply(z1.add(z3).subtract(14).subtract(e0sq.multiply(6)));
             sgh = s4.multiply(zn).multiply(z31.add(z33).subtract(6));
             if (tle.getI().getReal() < pi.divide(60.0).getReal()) {
                 // inclination smaller than 3 degrees
                 sh = zero;
             } else {
                 sh = s2.multiply(-zn).multiply(z21.add(z23));
             }
             ee2  = s1.multiply(s6).multiply(2);
             e3   = s1.multiply(s7).multiply(2);
             xi2  = s2.multiply(z12).multiply(2);
             xi3  = s2.multiply(z13.subtract(z11)).multiply(2);
             xl2  = s3.multiply(z2).multiply(-2);
             xl3  = s3.multiply(z3.subtract(z1)).multiply(-2);
             xl4  = s3.multiply(e0sq.multiply(-9).add(-21)).multiply(ze).multiply(-2);
             xgh2 = s4.multiply(z32).multiply(2);
             xgh3 = s4.multiply(z33.subtract(z31)).multiply(2);
             xgh4 = s4.multiply(ze).multiply(-18);
             xh2  = s2.multiply(z22).multiply(-2);
             xh3  = s2.multiply(z23.subtract(z21)).multiply(-2);
 
             if (iteration == 0) { // we compute lunar terms only on the first pass:
                 sse = se;
                 ssi = si;
                 ssl = sl;
                 ssh = (tle.getI().getReal() < pi.divide(60.0).getReal()) ? zero : sh.divide(sini0);
                 ssg = sgh.subtract(cosi0.multiply(ssh));
                 se2 = ee2;
                 si2 = xi2;
                 sl2 = xl2;
                 sgh2 = xgh2;
                 sh2 = xh2;
                 se3 = e3;
                 si3 = xi3;
                 sl3 = xl3;
                 sgh3 = xgh3;
                 sh3 = xh3;
                 sl4 = xl4;
                 sgh4 = xgh4;
                 zcosg = zero.newInstance(zcosgl);
                 zsing = zero.newInstance(zsingl);
                 zcosi = zero.newInstance(zcosil);
                 zsini = zero.newInstance(zsinil);
                 zcosh = cosq.multiply(zcoshl).add(sinq.multiply(zsinhl));
                 zsinh = sinq.multiply(zcoshl).subtract(cosq.multiply(zsinhl));
                 zn = TLEConstants.ZNL;
                 cc = TLEConstants.C1L;
                 ze = TLEConstants.ZEL;
             }
         } // end of solar - lunar - solar terms computation
 
         sse = sse.add(se);
         ssi = ssi.add(si);
         ssl = ssl.add(sl);
         ssg = ssg.add(sgh).subtract((tle.getI().getReal() < pi.divide(60.0).getReal()) ? zero : (cosi0.divide(sini0).multiply(sh)));
         ssh = ssh.add((tle.getI().getReal() < pi.divide(60.0).getReal()) ? zero : sh.divide(sini0));
 
 
 
         //        Start the resonant-synchronous tests and initialization
 
         T bfact = zero;
 
         // if mean motion is 1.893053 to 2.117652 revs/day, and eccentricity >= 0.5,
         // start of the 12-hour orbit, e > 0.5 section
         if (xnq.getReal() >= 0.00826 && xnq.getReal() <= 0.00924 && tle.getE().getReal() >= 0.5) {
 
             final T g201  = tle.getE().subtract(0.64).negate().multiply(0.440).add(-0.306);
             final T eoc   = tle.getE().multiply(e0sq);
             final T sini2 = sini0.multiply(sini0);
             final T f220  = cosi0.multiply(2).add(theta2).add(1).multiply(0.75);
             final T f221  = sini2.multiply(1.5);
             final T f321  = sini0.multiply(1.875).multiply(cosi0.multiply(2).negate().subtract(theta2.multiply(3)).add(1));
             final T f322  = sini0.multiply(-1.875).multiply(cosi0.multiply(2).subtract(theta2.multiply(3)).add(1));
             final T f441  = sini2.multiply(f220).multiply(35);
             final T f442  = sini2.multiply(sini2).multiply(39.3750);
             final T f522  = sini0.multiply(9.84375).multiply(sini2.multiply(cosi0.multiply(-2).add(theta2.multiply(-5)).add(1.0)).add(
                                     cosi0.multiply(4.0).add(theta2.multiply(6.0)).add(-2).multiply(0.33333333)));
             final T f523  = sini0.multiply(sini2.multiply(cosi0.multiply(-4).add(theta2.multiply(10)).add(-2)).multiply(4.92187512).add(
                                     cosi0.multiply(2).subtract(theta2.multiply(3)).add(1).multiply(6.56250012)));
             final T f542  = sini0.multiply(29.53125).multiply(cosi0.multiply(-8).add(2).add(
                                     theta2.multiply(cosi0.multiply(8).add(theta2.multiply(10)).add(-12))));
             final T f543  = sini0.multiply(29.53125).multiply(cosi0.multiply(-8).add(-2).add(
                                     theta2.multiply(cosi0.multiply(8).subtract(theta2.multiply(10)).add(12))));
             final T g211;
             final T g310;
             final T g322;
             final T g410;
             final T g422;
             final T g520;
 
             resonant = true;       // it is resonant...
             synchronous = false;     // but it's not synchronous
 
             // Geopotential resonance initialization for 12 hour orbits :
             if (tle.getE().getReal() <= 0.65) {
                 g211 = tle.getE().multiply( -13.247).add(  e0sq.multiply(   16.290)).add(                                  3.616);
                 g310 = tle.getE().multiply( 117.390).add(  e0sq.multiply( -228.419)).add(  eoc.multiply( 156.591)).add(  -19.302);
                 g322 = tle.getE().multiply(109.7927).add(  e0sq.multiply(-214.6334)).add(  eoc.multiply(146.5816)).add( -18.9068);
                 g410 = tle.getE().multiply( 242.694).add(  e0sq.multiply( -471.094)).add(  eoc.multiply( 313.953)).add(  -41.122);
                 g422 = tle.getE().multiply( 841.880).add(  e0sq.multiply(-1629.014)).add(  eoc.multiply(1083.435)).add( -146.407);
                 g520 = tle.getE().multiply(3017.977).add(  e0sq.multiply(-5740.032)).add(  eoc.multiply(3708.276)).add( -532.114);
             } else  {
                 g211 = tle.getE().multiply( 331.819).add(  e0sq.multiply( -508.738)).add(  eoc.multiply( 266.724)).add(  -72.099);
                 g310 = tle.getE().multiply(1582.851).add(  e0sq.multiply(-2415.925)).add(  eoc.multiply(1246.113)).add( -346.844);
                 g322 = tle.getE().multiply(1554.908).add(  e0sq.multiply(-2366.899)).add(  eoc.multiply(1215.972)).add( -342.585);
                 g410 = tle.getE().multiply(4758.686).add(  e0sq.multiply(-7193.992)).add(  eoc.multiply(3651.957)).add(-1052.797);
                 g422 = tle.getE().multiply(16178.11).add(  e0sq.multiply(-24462.77)).add(  eoc.multiply(12422.52)).add( -3581.69);
                 if (tle.getE().getReal() <= 0.715) {
                     g520 = tle.getE().multiply(-4664.75).add(  e0sq.multiply(  3763.64)).add(                                1464.74);
                 } else {
                     g520 = tle.getE().multiply(29936.92).add(  e0sq.multiply(-54087.36)).add(  eoc.multiply(31324.56)).add( -5149.66);
                 }
             }
 
             final T g533;
             final T g521;
             final T g532;
             if (tle.getE().getReal() < 0.7) {
                 g533 = tle.getE().multiply(  4988.61).add(  e0sq.multiply(  -9064.77)).add(  eoc.multiply(  5542.21)).add(  -919.2277);
                 g521 = tle.getE().multiply(4568.6173).add(  e0sq.multiply(-8491.4146)).add(  eoc.multiply( 5337.524)).add( -822.71072);
                 g532 = tle.getE().multiply(  4690.25).add(  e0sq.multiply(  -8624.77)).add(  eoc.multiply(   5341.4)).add(   -853.666);
             } else {
                 g533 = tle.getE().multiply(161616.52).add(  e0sq.multiply( -229838.2)).add(  eoc.multiply(109377.94)).add(  -37995.78);
                 g521 = tle.getE().multiply(218913.95).add(  e0sq.multiply(-309468.16)).add(  eoc.multiply(146349.42)).add( -51752.104);
                 g532 = tle.getE().multiply(170470.89).add(  e0sq.multiply(-242699.48)).add(  eoc.multiply(115605.82)).add(  -40023.88);
             }
 
             T temp1 = xnq.multiply(xnq).multiply(aqnv).multiply(aqnv).multiply(3);
             T temp  = temp1.multiply(TLEConstants.ROOT22);
             d2201   = temp.multiply(f220).multiply(g201);
             d2211   = temp.multiply(f221).multiply(g211);
             temp1   = temp1.multiply(aqnv);
             temp    = temp1.multiply(TLEConstants.ROOT32);
             d3210   = temp.multiply(f321).multiply(g310);
             d3222   = temp.multiply(f322).multiply(g322);
             temp1   = temp1.multiply(aqnv);
             temp    = temp1.multiply(2 * TLEConstants.ROOT44);
             d4410   = temp.multiply(f441).multiply(g410);
             d4422   = temp.multiply(f442).multiply(g422);
             temp1   = temp1.multiply(aqnv);
             temp    = temp1.multiply(TLEConstants.ROOT52);
             d5220   = temp.multiply(f522).multiply(g520);
             d5232   = temp.multiply(f523).multiply(g532);
             temp    = temp1.multiply(2 * TLEConstants.ROOT54);
             d5421   = temp.multiply(f542).multiply(g521);
             d5433   = temp.multiply(f543).multiply(g533);
             xlamo   = tle.getMeanAnomaly().add(tle.getRaan()).add(tle.getRaan()).subtract(thgr + thgr);
             bfact   = xmdot.add(xnodot).add(xnodot).subtract(TLEConstants.THDT + TLEConstants.THDT);
             bfact   = bfact.add(ssl).add(ssh).add(ssh);
         } else if (xnq.getReal() < 0.0052359877 && xnq.getReal() > 0.0034906585) {
             // if mean motion is .8 to 1.2 revs/day : (geosynch)
 
             final T cosio_plus_1 = cosi0.add(1.0);
             final T g200 = e0sq.multiply(e0sq.multiply(0.8125).add(-2.5)).add(1);
             final T g300 = e0sq.multiply(e0sq.multiply(6.60937).add(-6)).add(1);
             final T f311 = sini0.multiply(0.9375).multiply(sini0.multiply(cosi0.multiply(3).add(1))).subtract(cosio_plus_1.multiply(0.75));
             final T g310 = e0sq.multiply(2).add(1);
             final T f220 = cosio_plus_1.multiply(cosio_plus_1).multiply(0.75);
             final T f330 = f220.multiply(cosio_plus_1).multiply(2.5);
 
             resonant = true;
             synchronous = true;
 
             // Synchronous resonance terms initialization
             del1 = xnq.multiply(xnq).multiply(aqnv).multiply(aqnv).multiply(3);
             del2 = del1.multiply(f220).multiply(g200).multiply(2 * TLEConstants.Q22);
             del3 = del1.multiply(f330).multiply(g300).multiply(aqnv).multiply(3 * TLEConstants.Q33);
             del1 = del1.multiply(f311).multiply(g310).multiply(TLEConstants.Q31).multiply(aqnv);
             xlamo = tle.getMeanAnomaly().add(tle.getRaan()).add(tle.getPerigeeArgument()).subtract(thgr);
             bfact = xmdot.add(omgdot).add(xnodot).subtract(TLEConstants.THDT);
             bfact = bfact.add(ssl).add(ssg).add(ssh);
         } else {
             // it's neither a high-e 12-hours orbit nor a geosynchronous:
             resonant = false;
             synchronous = false;
         }
 
         if (resonant) {
             xfact = bfact.subtract(xnq);
 
             // Initialize integrator
             xli   = xlamo;
             xni   = xnq;
             atime = zero;
         }
         derivs = MathArrays.buildArray(xnq.getField(), 2);
     }
 
     /** Computes secular terms from current coordinates and epoch.
      * @param t offset from initial epoch (minutes)
      */
     protected void deepSecularEffects(final T t)  {
 
         xll     = xll.add(ssl.multiply(t));
         omgadf  = omgadf.add(ssg.multiply(t));
         xnode   = xnode.add(ssh.multiply(t));
         em      = tle.getE().add(sse.multiply(t));
         xinc    = tle.getI().add(ssi.multiply(t));
 
         if (resonant) {
             // If we're closer to t = 0 than to the currently-stored data
             // from the previous call to this function,  then we're
             // better off "restarting",  going back to the initial data.
             // The Dundee code rigs things up to _always_ take 720-minute
             // steps from epoch to end time,  except for the final step.
             // Easiest way to arrange similar behavior in this code is
             // just to always do a restart,  if we're in Dundee-compliant
             // mode.
             if (FastMath.abs(t).getReal() < FastMath.abs(t.subtract(atime)).getReal() || isDundeeCompliant)  {
                 // Epoch restart
                 atime = t.getField().getZero();
                 xni = xnq;
                 xli = xlamo;
             }
             boolean lastIntegrationStep = false;
             // if |step|>|step max| then do one step at step max
             while (!lastIntegrationStep) {
                 double delt = t.subtract(atime).getReal();
                 if (delt > SECULAR_INTEGRATION_STEP) {
                     delt = SECULAR_INTEGRATION_STEP;
                 } else if (delt < -SECULAR_INTEGRATION_STEP) {
                     delt = -SECULAR_INTEGRATION_STEP;
                 } else {
                     lastIntegrationStep = true;
                 }
 
                 computeSecularDerivs();
 
                 final T xldot = xni.add(xfact);
 
                 T xlpow = t.getField().getOne();
                 xli = xli.add(xldot.multiply(delt));
                 xni = xni.add(derivs[0].multiply(delt));
                 double delt_factor = delt;
                 xlpow = xlpow.multiply(xldot);
                 derivs[1] = derivs[1].multiply(xlpow);
                 delt_factor *= delt / 2;
                 xli = xli.add(derivs[0].multiply(delt_factor));
                 xni = xni.add(derivs[1].multiply(delt_factor));
                 atime = atime.add(delt);
             }
             xn = xni;
             final T temp = xnode.negate().add(thgr).add(t.multiply(TLEConstants.THDT));
             xll = xli.add(temp).add(synchronous ? omgadf.negate() : temp);
         }
     }
 
     /** Computes periodic terms from current coordinates and epoch.
      * @param t offset from initial epoch (min)
      */
     protected void deepPeriodicEffects(final T t)  {
 
         // If the time didn't change by more than 30 minutes,
         // there's no good reason to recompute the perturbations;
         // they don't change enough over so short a time span.
         // However,  the Dundee code _always_ recomputes,  so if
         // we're attempting to replicate its results,  we've gotta
         // recompute everything,  too.
         if (FastMath.abs(savtsn.subtract(t).getReal()) >= 30.0 || isDundeeCompliant)  {
 
             savtsn = t;
 
             // Update solar perturbations for time T
             T zm = t.multiply(TLEConstants.ZNS).add(zmos);
             T zf = zm.add(FastMath.sin(zm).multiply(2 * TLEConstants.ZES));
             FieldSinCos<T> sczf = FastMath.sinCos(zf);
             T sinzf = sczf.sin();
             T f2 = sinzf.multiply(sinzf).multiply(0.5).subtract(0.25);
             T f3 = sinzf.multiply(sczf.cos()).multiply(-0.5);
             final T ses = se2.multiply(f2).add(se3.multiply(f3));
             final T sis = si2.multiply(f2).add(si3.multiply(f3));
             final T sls = sl2.multiply(f2).add(sl3.multiply(f3)).add(sl4.multiply(sinzf));
             final T sghs = sgh2.multiply(f2).add(sgh3.multiply(f3)).add(sgh4.multiply(sinzf));
             final T shs = sh2.multiply(f2).add(sh3.multiply(f3));
 
             // Update lunar perturbations for time T
             zm = t.multiply(TLEConstants.ZNL).add(zmol);
             zf = zm.add(FastMath.sin(zm).multiply(2 * TLEConstants.ZEL));
             sczf = FastMath.sinCos(zf);
             sinzf = sczf.sin();
             f2 =  sinzf.multiply(sinzf).multiply(0.5).subtract(0.25);
             f3 = sinzf.multiply(sczf.cos()).multiply(-0.5);
             final T sel = ee2.multiply(f2).add(e3.multiply(f3));
             final T sil = xi2.multiply(f2).add(xi3.multiply(f3));
             final T sll = xl2.multiply(f2).add(xl3.multiply(f3)).add(xl4.multiply(sinzf));
             final T sghl = xgh2.multiply(f2).add(xgh3.multiply(f3)).add(xgh4.multiply(sinzf));
             final T sh1 = xh2.multiply(f2).add(xh3.multiply(f3));
 
             // Sum the solar and lunar contributions
             pe   = ses.add(sel);
             pinc = sis.add(sil);
             pl   = sls.add(sll);
             pgh  = sghs.add(sghl);
             ph   = shs.add(sh1);
         }
 
         xinc = xinc.add(pinc);
 
         final FieldSinCos<T> scis = FastMath.sinCos(xinc);
         final T sinis = scis.sin();
         final T cosis = scis.cos();
 
         /* Add solar/lunar perturbation correction to eccentricity: */
         em     = em.add(pe);
         xll    = xll.add(pl);
         omgadf = omgadf.add(pgh);
         xinc   = MathUtils.normalizeAngle(xinc, t.getField().getZero());
 
         if (FastMath.abs(xinc).getReal() >= 0.2) {
             // Apply periodics directly
             final T temp_val = ph.divide(sinis);
             omgadf = omgadf.subtract(cosis.multiply(temp_val));
             xnode  = xnode.add(temp_val);
         } else {
             // Apply periodics with Lyddane modification
             final FieldSinCos<T> scok = FastMath.sinCos(xnode);
             final T sinok = scok.sin();
             final T cosok = scok.cos();
             final T alfdp =  ph.multiply(cosok).add((pinc.multiply(cosis).add(sinis)).multiply(sinok));
             final T betdp = ph.negate().multiply(sinok).add((pinc.multiply(cosis).add(sinis)).multiply(cosok));
             final T delta_xnode = MathUtils.normalizeAngle(FastMath.atan2(alfdp, betdp).subtract(xnode), t.getField().getZero());
             final T dls = xnode.negate().multiply(sinis).multiply(pinc);
             omgadf = omgadf.add(dls.subtract(cosis.multiply(delta_xnode)));
             xnode  = xnode.add(delta_xnode);
         }
     }
 
     /** Computes internal secular derivs. */
     private void computeSecularDerivs() {
 
         final FieldSinCos<T> sc_li  = FastMath.sinCos(xli);
         final T sin_li = sc_li.sin();
         final T cos_li = sc_li.cos();
         final T sin_2li = sin_li.multiply(cos_li).multiply(2.);
         final T cos_2li = cos_li.multiply(cos_li).multiply(2.).subtract(1.);
 
         // Dot terms calculated :
         if (synchronous)  {
             final T sin_3li = sin_2li.multiply(cos_li).add(cos_2li.multiply(sin_li));
             final T cos_3li = cos_2li.multiply(cos_li).subtract(sin_2li.multiply(sin_li));
             final T term1a = del1.multiply(sin_li .multiply(TLEConstants.C_FASX2) .subtract(cos_li .multiply(TLEConstants.S_FASX2 )));
             final T term2a = del2.multiply(sin_2li.multiply(TLEConstants.C_2FASX4).subtract(cos_2li.multiply(TLEConstants.S_2FASX4)));
             final T term3a = del3.multiply(sin_3li.multiply(TLEConstants.C_3FASX6).subtract(cos_3li.multiply(TLEConstants.S_3FASX6)));
             final T term1b = del1.multiply(cos_li .multiply(TLEConstants.C_FASX2)      .add(sin_li .multiply(TLEConstants.S_FASX2 )));
             final T term2b = del2.multiply(cos_2li.multiply(TLEConstants.C_2FASX4)     .add(sin_2li.multiply(TLEConstants.S_2FASX4))).multiply(2.0);
             final T term3b = del3.multiply(cos_3li.multiply(TLEConstants.C_3FASX6)     .add(sin_3li.multiply(TLEConstants.S_3FASX6))).multiply(3.0);
             derivs[0] = term1a.add(term2a).add(term3a);
             derivs[1] = term1b.add(term2b).add(term3b);
         } else {
             // orbit is a 12-hour resonant one
             final T xomi = omegaq.add(omgdot.multiply(atime));
             final FieldSinCos<T> sc_omi  = FastMath.sinCos(xomi);
             final T sin_omi = sc_omi.sin();
             final T cos_omi = sc_omi.cos();
             final T sin_li_m_omi = sin_li.multiply(cos_omi).subtract(sin_omi.multiply(cos_li));
             final T sin_li_p_omi = sin_li.multiply(cos_omi).add(     sin_omi.multiply(cos_li));
             final T cos_li_m_omi = cos_li.multiply(cos_omi).add(     sin_omi.multiply(sin_li));
             final T cos_li_p_omi = cos_li.multiply(cos_omi).subtract(sin_omi.multiply(sin_li));
             final T sin_2omi = sin_omi.multiply(cos_omi).multiply(2.0);
             final T cos_2omi = cos_omi.multiply(cos_omi).multiply(2.0).subtract(1.0);
             final T sin_2li_m_omi  = sin_2li.multiply(cos_omi ).subtract(sin_omi .multiply(cos_2li));
             final T sin_2li_p_omi  = sin_2li.multiply(cos_omi ).add(     sin_omi .multiply(cos_2li));
             final T cos_2li_m_omi  = cos_2li.multiply(cos_omi ).add(     sin_omi .multiply(sin_2li));
             final T cos_2li_p_omi  = cos_2li.multiply(cos_omi ).subtract(sin_omi .multiply(sin_2li));
             final T sin_2li_p_2omi = sin_2li.multiply(cos_2omi).add(     sin_2omi.multiply(cos_2li));
             final T cos_2li_p_2omi = cos_2li.multiply(cos_2omi).subtract(sin_2omi.multiply(sin_2li));
             final T sin_2omi_p_li  = sin_li .multiply(cos_2omi).add(     sin_2omi.multiply(cos_li ));
             final T cos_2omi_p_li  = cos_li .multiply(cos_2omi).subtract(sin_2omi.multiply(sin_li ));
             final T term1a = d2201.multiply(sin_2omi_p_li .multiply(TLEConstants.C_G22).subtract(cos_2omi_p_li .multiply(TLEConstants.S_G22))) .add(
                              d2211.multiply(sin_li        .multiply(TLEConstants.C_G22).subtract(cos_li        .multiply(TLEConstants.S_G22)))).add(
                              d3210.multiply(sin_li_p_omi  .multiply(TLEConstants.C_G32).subtract(cos_li_p_omi  .multiply(TLEConstants.S_G32)))).add(
                              d3222.multiply(sin_li_m_omi  .multiply(TLEConstants.C_G32).subtract(cos_li_m_omi  .multiply(TLEConstants.S_G32)))).add(
                              d5220.multiply(sin_li_p_omi  .multiply(TLEConstants.C_G52).subtract(cos_li_p_omi  .multiply(TLEConstants.S_G52)))).add(
                              d5232.multiply(sin_li_m_omi  .multiply(TLEConstants.C_G52).subtract(cos_li_m_omi  .multiply(TLEConstants.S_G52))));
             final T term2a = d4410.multiply(sin_2li_p_2omi.multiply(TLEConstants.C_G44).subtract(cos_2li_p_2omi.multiply(TLEConstants.S_G44))) .add(
                              d4422.multiply(sin_2li       .multiply(TLEConstants.C_G44).subtract(cos_2li       .multiply(TLEConstants.S_G44)))).add(
                              d5421.multiply(sin_2li_p_omi .multiply(TLEConstants.C_G54).subtract(cos_2li_p_omi .multiply(TLEConstants.S_G54)))).add(
                              d5433.multiply(sin_2li_m_omi .multiply(TLEConstants.C_G54).subtract(cos_2li_m_omi .multiply(TLEConstants.S_G54))));
             final T term1b = d2201.multiply(cos_2omi_p_li .multiply(TLEConstants.C_G22)     .add(sin_2omi_p_li .multiply(TLEConstants.S_G22))) .add(
                              d2211.multiply(cos_li        .multiply(TLEConstants.C_G22)     .add(sin_li        .multiply(TLEConstants.S_G22)))).add(
                              d3210.multiply(cos_li_p_omi  .multiply(TLEConstants.C_G32)     .add(sin_li_p_omi  .multiply(TLEConstants.S_G32)))).add(
                              d3222.multiply(cos_li_m_omi  .multiply(TLEConstants.C_G32)     .add(sin_li_m_omi  .multiply(TLEConstants.S_G32)))).add(
                              d5220.multiply(cos_li_p_omi  .multiply(TLEConstants.C_G52)     .add(sin_li_p_omi  .multiply(TLEConstants.S_G52)))).add(
                              d5232.multiply(cos_li_m_omi  .multiply(TLEConstants.C_G52)     .add(sin_li_m_omi  .multiply(TLEConstants.S_G52))));
             final T term2b = d4410.multiply(cos_2li_p_2omi.multiply(TLEConstants.C_G44)     .add(sin_2li_p_2omi.multiply(TLEConstants.S_G44))) .add(
                              d4422.multiply(cos_2li       .multiply(TLEConstants.C_G44)     .add(sin_2li       .multiply(TLEConstants.S_G44)))).add(
                              d5421.multiply(cos_2li_p_omi .multiply(TLEConstants.C_G54)     .add(sin_2li_p_omi .multiply(TLEConstants.S_G54)))).add(
                              d5433.multiply(cos_2li_m_omi .multiply(TLEConstants.C_G54)     .add(sin_2li_m_omi .multiply(TLEConstants.S_G54)))).multiply(2.0);
 
             derivs[0] = term1a.add(term2a);
             derivs[1] = term1b.add(term2b);
 
         }
     }
 
 }