DeepSDP4.java

  1. /* Copyright 2002-2013 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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;

  19. import org.apache.commons.math3.util.FastMath;
  20. import org.apache.commons.math3.util.MathUtils;
  21. import org.orekit.attitudes.AttitudeProvider;
  22. import org.orekit.errors.OrekitException;
  23. import org.orekit.time.AbsoluteDate;
  24. import org.orekit.time.TimeScalesFactory;
  25. import org.orekit.utils.Constants;


  28. /** This class contains the methods that compute deep space perturbation terms.
  29.  * <p>
  30.  * The user should not bother in this class since it is handled internaly by the
  31.  * {@link TLEPropagator}.
  32.  * </p>
  33.  * <p>This implementation is largely inspired from the paper and source code <a
  34.  * href="http://www.celestrak.com/publications/AIAA/2006-6753/">Revisiting Spacetrack
  35.  * Report #3</a> and is fully compliant with its results and tests cases.</p>
  36.  * @author Felix R. Hoots, Ronald L. Roehrich, December 1980 (original fortran)
  37.  * @author David A. Vallado, Paul Crawford, Richard Hujsak, T.S. Kelso (C++ translation and improvements)
  38.  * @author Fabien Maussion (java translation)
  39.  */
  40. class DeepSDP4 extends SDP4 {

  42.     // CHECKSTYLE: stop JavadocVariable check

  44.     // Internal constants
  45.     private static final double ZNS      = 1.19459E-5;
  46.     private static final double ZES      = 0.01675;
  47.     private static final double ZNL      = 1.5835218E-4;
  48.     private static final double ZEL      = 0.05490;
  49.     private static final double THDT     = 4.3752691E-3;
  50.     private static final double C1SS     =  2.9864797E-6;
  51.     private static final double C1L      = 4.7968065E-7;

  53.     private static final double ROOT22   = 1.7891679E-6;
  54.     private static final double ROOT32   = 3.7393792E-7;
  55.     private static final double ROOT44   = 7.3636953E-9;
  56.     private static final double ROOT52   = 1.1428639E-7;
  57.     private static final double ROOT54   = 2.1765803E-9;

  59.     private static final double Q22      =  1.7891679E-6;
  60.     private static final double Q31      =  2.1460748E-6;
  61.     private static final double Q33      =  2.2123015E-7;

  63.     private static final double C_FASX2  =  0.99139134268488593;
  64.     private static final double S_FASX2  =  0.13093206501640101;
  65.     private static final double C_2FASX4 =  0.87051638752972937;
  66.     private static final double S_2FASX4 = -0.49213943048915526;
  67.     private static final double C_3FASX6 =  0.43258117585763334;
  68.     private static final double S_3FASX6 =  0.90159499016666422;

  70.     private static final double C_G22    =  0.87051638752972937;
  71.     private static final double S_G22    = -0.49213943048915526;
  72.     private static final double C_G32    =  0.57972190187001149;
  73.     private static final double S_G32    =  0.81481440616389245;
  74.     private static final double C_G44    = -0.22866241528815548;
  75.     private static final double S_G44    =  0.97350577801807991;
  76.     private static final double C_G52    =  0.49684831179884198;
  77.     private static final double S_G52    =  0.86783740128127729;
  78.     private static final double C_G54    = -0.29695209575316894;
  79.     private static final double S_G54    = -0.95489237761529999;

  81.     /** Integration step (seconds). */
  82.     private static final double SECULAR_INTEGRATION_STEP  = 720.0;

  84.     /** Integration order. */
  85.     private static final int    SECULAR_INTEGRATION_ORDER = 2;

  87.     /** Intermediate values. */
  88.     private double thgr;
  89.     private double xnq;
  90.     private double omegaq;
  91.     private double zcosil;
  92.     private double zsinil;
  93.     private double zsinhl;
  94.     private double zcoshl;
  95.     private double zmol;
  96.     private double zcosgl;
  97.     private double zsingl;
  98.     private double zmos;
  99.     private double savtsn;

  101.     private double ee2;
  102.     private double e3;
  103.     private double xi2;
  104.     private double xi3;
  105.     private double xl2;
  106.     private double xl3;
  107.     private double xl4;
  108.     private double xgh2;
  109.     private double xgh3;
  110.     private double xgh4;
  111.     private double xh2;
  112.     private double xh3;

  114.     private double d2201;
  115.     private double d2211;
  116.     private double d3210;
  117.     private double d3222;
  118.     private double d4410;
  119.     private double d4422;
  120.     private double d5220;
  121.     private double d5232;
  122.     private double d5421;
  123.     private double d5433;
  124.     private double xlamo;

  126.     private double sse;
  127.     private double ssi;
  128.     private double ssl;
  129.     private double ssh;
  130.     private double ssg;
  131.     private double se2;
  132.     private double si2;
  133.     private double sl2;
  134.     private double sgh2;
  135.     private double sh2;
  136.     private double se3;
  137.     private double si3;
  138.     private double sl3;
  139.     private double sgh3;
  140.     private double sh3;
  141.     private double sl4;
  142.     private double sgh4;

  144.     private double del1;
  145.     private double del2;
  146.     private double del3;
  147.     private double xfact;
  148.     private double xli;
  149.     private double xni;
  150.     private double atime;

  152.     private double pe;
  153.     private double pinc;
  154.     private double pl;
  155.     private double pgh;
  156.     private double ph;

  158.     private double[] derivs;

  160.     // CHECKSTYLE: resume JavadocVariable check

  162.     /** Flag for resonant orbits. */
  163.     private boolean resonant;

  165.     /** Flag for synchronous orbits. */
  166.     private boolean synchronous;

  168.     /** Flag for compliance with Dundee modifications. */
  169.     private boolean isDundeeCompliant = true;

  171.     /** Constructor for a unique initial TLE.
  172.      * @param initialTLE the TLE to propagate.
  173.      * @param attitudeProvider provider for attitude computation
  174.      * @param mass spacecraft mass (kg)
  175.      * @exception OrekitException if some specific error occurs
  176.      */
  177.     protected DeepSDP4(final TLE initialTLE, final AttitudeProvider attitudeProvider,
  178.                        final double mass) throws OrekitException {
  179.         super(initialTLE, attitudeProvider, mass);
  180.     }

  182.     /** Computes luni - solar terms from initial coordinates and epoch.
  183.      * @exception OrekitException when UTC time steps can't be read
  184.      */
  185.     protected void luniSolarTermsComputation() throws OrekitException {

  187.         final double sing = FastMath.sin(tle.getPerigeeArgument());
  188.         final double cosg = FastMath.cos(tle.getPerigeeArgument());

  190.         final double sinq = FastMath.sin(tle.getRaan());
  191.         final double cosq = FastMath.cos(tle.getRaan());
  192.         final double aqnv = 1.0 / a0dp;

  194.         // Compute julian days since 1900
  195.         final double daysSince1900 =
  196.             (tle.getDate().durationFrom(AbsoluteDate.JULIAN_EPOCH) +
  197.              tle.getDate().timeScalesOffset(TimeScalesFactory.getUTC(), TimeScalesFactory.getTT())) / Constants.JULIAN_DAY - 2415020;


  200.         double cc = C1SS;
  201.         double ze = ZES;
  202.         double zn = ZNS;
  203.         double zsinh = sinq;
  204.         double zcosh = cosq;

  206.         thgr = thetaG(tle.getDate());
  207.         xnq = xn0dp;
  208.         omegaq = tle.getPerigeeArgument();

  210.         final double xnodce = 4.5236020 - 9.2422029e-4 * daysSince1900;
  211.         final double stem = FastMath.sin(xnodce);
  212.         final double ctem = FastMath.cos(xnodce);
  213.         final double c_minus_gam = 0.228027132 * daysSince1900 - 1.1151842;
  214.         final double gam = 5.8351514 + 0.0019443680 * daysSince1900;

  216.         zcosil = 0.91375164 - 0.03568096 * ctem;
  217.         zsinil = FastMath.sqrt(1.0 - zcosil * zcosil);
  218.         zsinhl = 0.089683511 * stem / zsinil;
  219.         zcoshl = FastMath.sqrt(1.0 - zsinhl * zsinhl);
  220.         zmol = MathUtils.normalizeAngle(c_minus_gam, FastMath.PI);

  222.         double zx = 0.39785416 * stem / zsinil;
  223.         final double zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
  224.         zx = FastMath.atan2( zx, zy) + gam - xnodce;
  225.         zcosgl = FastMath.cos( zx);
  226.         zsingl = FastMath.sin( zx);
  227.         zmos = MathUtils.normalizeAngle(6.2565837 + 0.017201977 * daysSince1900, FastMath.PI);

  229.         // Do solar terms
  230.         savtsn = 1e20;

  232.         double zcosi =  0.91744867;
  233.         double zsini =  0.39785416;
  234.         double zsing = -0.98088458;
  235.         double zcosg =  0.1945905;

  237.         double se = 0;
  238.         double sgh = 0;
  239.         double sh = 0;
  240.         double si = 0;
  241.         double sl = 0;

  243.         // There was previously some convoluted logic here, but it boils
  244.         // down to this:  we compute the solar terms,  then the lunar terms.
  245.         // On a second pass,  we recompute the solar terms, taking advantage
  246.         // of the improved data that resulted from computing lunar terms.
  247.         for (int iteration = 0; iteration < 2; ++iteration) {
  248.             final double a1 = zcosg * zcosh + zsing * zcosi * zsinh;
  249.             final double a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
  250.             final double a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
  251.             final double a8 = zsing * zsini;
  252.             final double a9 = zsing * zsinh + zcosg * zcosi * zcosh;
  253.             final double a10 = zcosg * zsini;
  254.             final double a2 = cosi0 * a7 + sini0 * a8;
  255.             final double a4 = cosi0 * a9 + sini0 * a10;
  256.             final double a5 = -sini0 * a7 + cosi0 * a8;
  257.             final double a6 = -sini0 * a9 + cosi0 * a10;
  258.             final double x1 = a1 * cosg + a2 * sing;
  259.             final double x2 = a3 * cosg + a4 * sing;
  260.             final double x3 = -a1 * sing + a2 * cosg;
  261.             final double x4 = -a3 * sing + a4 * cosg;
  262.             final double x5 = a5 * sing;
  263.             final double x6 = a6 * sing;
  264.             final double x7 = a5 * cosg;
  265.             final double x8 = a6 * cosg;
  266.             final double z31 = 12 * x1 * x1 - 3 * x3 * x3;
  267.             final double z32 = 24 * x1 * x2 - 6 * x3 * x4;
  268.             final double z33 = 12 * x2 * x2 - 3 * x4 * x4;
  269.             final double z11 = -6 * a1 * a5 + e0sq * (-24 * x1 * x7 - 6 * x3 * x5);
  270.             final double z12 = -6 * (a1 * a6 + a3 * a5) +
  271.                                e0sq * (-24 * (x2 * x7 + x1 * x8) - 6 * (x3 * x6 + x4 * x5));
  272.             final double z13 = -6 * a3 * a6 + e0sq * (-24 * x2 * x8 - 6 * x4 * x6);
  273.             final double z21 = 6 * a2 * a5 + e0sq * (24 * x1 * x5 - 6 * x3 * x7);
  274.             final double z22 = 6 * (a4 * a5 + a2 * a6) +
  275.                                e0sq * (24 * (x2 * x5 + x1 * x6) - 6 * (x4 * x7 + x3 * x8));
  276.             final double z23 = 6 * a4 * a6 + e0sq * (24 * x2 * x6 - 6 * x4 * x8);
  277.             final double s3 = cc / xnq;
  278.             final double s2 = -0.5 * s3 / beta0;
  279.             final double s4 = s3 * beta0;
  280.             final double s1 = -15 * tle.getE() * s4;
  281.             final double s5 = x1 * x3 + x2 * x4;
  282.             final double s6 = x2 * x3 + x1 * x4;
  283.             final double s7 = x2 * x4 - x1 * x3;
  284.             double z1 = 3 * (a1 * a1 + a2 * a2) + z31 * e0sq;
  285.             double z2 = 6 * (a1 * a3 + a2 * a4) + z32 * e0sq;
  286.             double z3 = 3 * (a3 * a3 + a4 * a4) + z33 * e0sq;

  288.             z1 = z1 + z1 + beta02 * z31;
  289.             z2 = z2 + z2 + beta02 * z32;
  290.             z3 = z3 + z3 + beta02 * z33;
  291.             se = s1 * zn * s5;
  292.             si = s2 * zn * (z11 + z13);
  293.             sl = -zn * s3 * (z1 + z3 - 14 - 6 * e0sq);
  294.             sgh = s4 * zn * (z31 + z33 - 6);
  295.             if (tle.getI() < (FastMath.PI / 60.0)) {
  296.                 // inclination smaller than 3 degrees
  297.                 sh = 0;
  298.             } else {
  299.                 sh = -zn * s2 * (z21 + z23);
  300.             }
  301.             ee2  =  2 * s1 * s6;
  302.             e3   =  2 * s1 * s7;
  303.             xi2  =  2 * s2 * z12;
  304.             xi3  =  2 * s2 * (z13 - z11);
  305.             xl2  = -2 * s3 * z2;
  306.             xl3  = -2 * s3 * (z3 - z1);
  307.             xl4  = -2 * s3 * (-21 - 9 * e0sq) * ze;
  308.             xgh2 =  2 * s4 * z32;
  309.             xgh3 =  2 * s4 * (z33 - z31);
  310.             xgh4 = -18 * s4 * ze;
  311.             xh2  = -2 * s2 * z22;
  312.             xh3  = -2 * s2 * (z23 - z21);

  314.             if (iteration == 0) { // we compute lunar terms only on the first pass:
  315.                 sse = se;
  316.                 ssi = si;
  317.                 ssl = sl;
  318.                 ssh = (tle.getI() < (FastMath.PI / 60.0)) ? 0 : sh / sini0;
  319.                 ssg = sgh - cosi0 * ssh;
  320.                 se2 = ee2;
  321.                 si2 = xi2;
  322.                 sl2 = xl2;
  323.                 sgh2 = xgh2;
  324.                 sh2 = xh2;
  325.                 se3 = e3;
  326.                 si3 = xi3;
  327.                 sl3 = xl3;
  328.                 sgh3 = xgh3;
  329.                 sh3 = xh3;
  330.                 sl4 = xl4;
  331.                 sgh4 = xgh4;
  332.                 zcosg = zcosgl;
  333.                 zsing = zsingl;
  334.                 zcosi = zcosil;
  335.                 zsini = zsinil;
  336.                 zcosh = zcoshl * cosq + zsinhl * sinq;
  337.                 zsinh = sinq * zcoshl - cosq * zsinhl;
  338.                 zn = ZNL;
  339.                 cc = C1L;
  340.                 ze = ZEL;
  341.             }
  342.         } // end of solar - lunar - solar terms computation

  344.         sse += se;
  345.         ssi += si;
  346.         ssl += sl;
  347.         ssg += sgh - ((tle.getI() < (FastMath.PI / 60.0)) ? 0 : (cosi0 / sini0 * sh));
  348.         ssh += (tle.getI() < (FastMath.PI / 60.0)) ? 0 : sh / sini0;



  352.         //        Start the resonant-synchronous tests and initialization

  354.         double bfact = 0;

  356.         // if mean motion is 1.893053 to 2.117652 revs/day, and eccentricity >= 0.5,
  357.         // start of the 12-hour orbit, e > 0.5 section
  358.         if ((xnq >= 0.00826) && (xnq <= 0.00924) && (tle.getE() >= 0.5)) {

  360.             final double g201 = -0.306 - (tle.getE() - 0.64) * 0.440;
  361.             final double eoc = tle.getE() * e0sq;
  362.             final double sini2 = sini0 * sini0;
  363.             final double f220 = 0.75 * (1 + 2 * cosi0 + theta2);
  364.             final double f221 = 1.5 * sini2;
  365.             final double f321 =  1.875 * sini0 * (1 - 2 * cosi0 - 3 * theta2);
  366.             final double f322 = -1.875 * sini0 * (1 + 2 * cosi0 - 3 * theta2);
  367.             final double f441 = 35 * sini2 * f220;
  368.             final double f442 = 39.3750 * sini2 * sini2;
  369.             final double f522 = 9.84375 * sini0 * (sini2 * (1 - 2 * cosi0 - 5 * theta2) +
  370.                                                    0.33333333 * (-2 + 4 * cosi0 + 6 * theta2));
  371.             final double f523 = sini0 * (4.92187512 * sini2 * (-2 - 4 * cosi0 + 10 * theta2) +
  372.                                          6.56250012 * (1 + 2 * cosi0 - 3 * theta2));
  373.             final double f542 = 29.53125 * sini0 * (2 - 8 * cosi0 + theta2 * (-12 + 8 * cosi0 + 10 * theta2));
  374.             final double f543 = 29.53125 * sini0 * (-2 - 8 * cosi0 + theta2 * (12 + 8 * cosi0 - 10 * theta2));
  375.             double g211;
  376.             double g310;
  377.             double g322;
  378.             double g410;
  379.             double g422;
  380.             double g520;

  382.             resonant = true;       // it is resonant...
  383.             synchronous = false;     // but it's not synchronous

  385.             // Geopotential resonance initialization for 12 hour orbits :
  386.             if (tle.getE() <= 0.65) {
  387.                 g211 =    3.616  -   13.247  * tle.getE() +   16.290  * e0sq;
  388.                 g310 =  -19.302  +  117.390  * tle.getE() -  228.419  * e0sq +  156.591  * eoc;
  389.                 g322 =  -18.9068 +  109.7927 * tle.getE() -  214.6334 * e0sq +  146.5816 * eoc;
  390.                 g410 =  -41.122  +  242.694  * tle.getE() -  471.094  * e0sq +  313.953  * eoc;
  391.                 g422 = -146.407  +  841.880  * tle.getE() - 1629.014  * e0sq + 1083.435  * eoc;
  392.                 g520 = -532.114  + 3017.977  * tle.getE() - 5740.032  * e0sq + 3708.276  * eoc;
  393.             } else  {
  394.                 g211 =   -72.099 +   331.819 * tle.getE() -   508.738 * e0sq +   266.724 * eoc;
  395.                 g310 =  -346.844 +  1582.851 * tle.getE() -  2415.925 * e0sq +  1246.113 * eoc;
  396.                 g322 =  -342.585 +  1554.908 * tle.getE() -  2366.899 * e0sq +  1215.972 * eoc;
  397.                 g410 = -1052.797 +  4758.686 * tle.getE() -  7193.992 * e0sq +  3651.957 * eoc;
  398.                 g422 = -3581.69  + 16178.11  * tle.getE() - 24462.77  * e0sq + 12422.52  * eoc;
  399.                 if (tle.getE() <= 0.715) {
  400.                     g520 = 1464.74 - 4664.75 * tle.getE() + 3763.64 * e0sq;
  401.                 } else {
  402.                     g520 = -5149.66 + 29936.92 * tle.getE() - 54087.36 * e0sq + 31324.56 * eoc;
  403.                 }
  404.             }

  406.             double g533;
  407.             double g521;
  408.             double g532;
  409.             if (tle.getE() < 0.7) {
  410.                 g533 = -919.2277  + 4988.61   * tle.getE() - 9064.77   * e0sq + 5542.21  * eoc;
  411.                 g521 = -822.71072 + 4568.6173 * tle.getE() - 8491.4146 * e0sq + 5337.524 * eoc;
  412.                 g532 = -853.666   + 4690.25   * tle.getE() - 8624.77   * e0sq + 5341.4   * eoc;
  413.             } else {
  414.                 g533 = -37995.78  + 161616.52 * tle.getE() - 229838.2  * e0sq + 109377.94 * eoc;
  415.                 g521 = -51752.104 + 218913.95 * tle.getE() - 309468.16 * e0sq + 146349.42 * eoc;
  416.                 g532 = -40023.88  + 170470.89 * tle.getE() - 242699.48 * e0sq + 115605.82 * eoc;
  417.             }

  419.             double temp1 = 3 * xnq * xnq * aqnv * aqnv;
  420.             double temp = temp1 * ROOT22;
  421.             d2201 = temp * f220 * g201;
  422.             d2211 = temp * f221 * g211;
  423.             temp1 *= aqnv;
  424.             temp = temp1 * ROOT32;
  425.             d3210 = temp * f321 * g310;
  426.             d3222 = temp * f322 * g322;
  427.             temp1 *= aqnv;
  428.             temp = 2 * temp1 * ROOT44;
  429.             d4410 = temp * f441 * g410;
  430.             d4422 = temp * f442 * g422;
  431.             temp1 *= aqnv;
  432.             temp = temp1 * ROOT52;
  433.             d5220 = temp * f522 * g520;
  434.             d5232 = temp * f523 * g532;
  435.             temp = 2 * temp1 * ROOT54;
  436.             d5421 = temp * f542 * g521;
  437.             d5433 = temp * f543 * g533;
  438.             xlamo = tle.getMeanAnomaly() + tle.getRaan() + tle.getRaan() - thgr - thgr;
  439.             bfact = xmdot + xnodot + xnodot - THDT - THDT;
  440.             bfact += ssl + ssh + ssh;
  441.         } else if ((xnq < 0.0052359877) && (xnq > 0.0034906585)) {
  442.             // if mean motion is .8 to 1.2 revs/day : (geosynch)

  444.             final double cosio_plus_1 = 1.0 + cosi0;
  445.             final double g200 = 1 + e0sq * (-2.5 + 0.8125  * e0sq);
  446.             final double g300 = 1 + e0sq * (-6   + 6.60937 * e0sq);
  447.             final double f311 = 0.9375 * sini0 * sini0 * (1 + 3 * cosi0) - 0.75 * cosio_plus_1;
  448.             final double g310 = 1 + 2 * e0sq;
  449.             final double f220 = 0.75 * cosio_plus_1 * cosio_plus_1;
  450.             final double f330 = 2.5 * f220 * cosio_plus_1;

  452.             resonant = true;
  453.             synchronous = true;

  455.             // Synchronous resonance terms initialization
  456.             del1 = 3 * xnq * xnq * aqnv * aqnv;
  457.             del2 = 2 * del1 * f220 * g200 * Q22;
  458.             del3 = 3 * del1 * f330 * g300 * Q33 * aqnv;
  459.             del1 = del1 * f311 * g310 * Q31 * aqnv;
  460.             xlamo = tle.getMeanAnomaly() + tle.getRaan() + tle.getPerigeeArgument() - thgr;
  461.             bfact = xmdot + omgdot + xnodot - THDT;
  462.             bfact = bfact + ssl + ssg + ssh;
  463.         } else {
  464.             // it's neither a high-e 12-hours orbit nor a geosynchronous:
  465.             resonant = false;
  466.             synchronous = false;
  467.         }

  469.         if (resonant) {
  470.             xfact = bfact - xnq;

  472.             // Initialize integrator
  473.             xli   = xlamo;
  474.             xni   = xnq;
  475.             atime = 0;
  476.         }
  477.         derivs = new double[SECULAR_INTEGRATION_ORDER];
  478.     }

  480.     /** Computes secular terms from current coordinates and epoch.
  481.      * @param t offset from initial epoch (minutes)
  482.      */
  483.     protected void deepSecularEffects(final double t)  {

  485.         xll    += ssl * t;
  486.         omgadf += ssg * t;
  487.         xnode  += ssh * t;
  488.         em      = tle.getE() + sse * t;
  489.         xinc    = tle.getI() + ssi * t;

  491.         if (resonant) {
  492.             // If we're closer to t = 0 than to the currently-stored data
  493.             // from the previous call to this function,  then we're
  494.             // better off "restarting",  going back to the initial data.
  495.             // The Dundee code rigs things up to _always_ take 720-minute
  496.             // steps from epoch to end time,  except for the final step.
  497.             // Easiest way to arrange similar behavior in this code is
  498.             // just to always do a restart,  if we're in Dundee-compliant
  499.             // mode.
  500.             if (FastMath.abs(t) < FastMath.abs(t - atime) || isDundeeCompliant)  {
  501.                 // Epoch restart
  502.                 atime = 0;
  503.                 xni = xnq;
  504.                 xli = xlamo;
  505.             }
  506.             boolean lastIntegrationStep = false;
  507.             // if |step|>|step max| then do one step at step max
  508.             while (!lastIntegrationStep) {
  509.                 double delt = t - atime;
  510.                 if (delt > SECULAR_INTEGRATION_STEP) {
  511.                     delt = SECULAR_INTEGRATION_STEP;
  512.                 } else if (delt < -SECULAR_INTEGRATION_STEP) {
  513.                     delt = -SECULAR_INTEGRATION_STEP;
  514.                 } else {
  515.                     lastIntegrationStep = true;
  516.                 }

  518.                 computeSecularDerivs();

  520.                 final double xldot = xni + xfact;

  522.                 double xlpow = 1.;
  523.                 xli += delt * xldot;
  524.                 xni += delt * derivs[0];
  525.                 double delt_factor = delt;
  526.                 for (int j = 2; j <= SECULAR_INTEGRATION_ORDER; ++j) {
  527.                     xlpow *= xldot;
  528.                     derivs[j - 1] *= xlpow;
  529.                     delt_factor *= delt / (double) j;
  530.                     xli += delt_factor * derivs[j - 2];
  531.                     xni += delt_factor * derivs[j - 1];
  532.                 }
  533.                 atime += delt;
  534.             }
  535.             xn = xni;
  536.             final double temp = -xnode + thgr + t * THDT;
  537.             xll = xli + temp + (synchronous ? -omgadf : temp);
  538.         }
  539.     }

  541.     /** Computes periodic terms from current coordinates and epoch.
  542.      * @param t offset from initial epoch (min)
  543.      */
  544.     protected void deepPeriodicEffects(final double t)  {

  546.         // If the time didn't change by more than 30 minutes,
  547.         // there's no good reason to recompute the perturbations;
  548.         // they don't change enough over so short a time span.
  549.         // However,  the Dundee code _always_ recomputes,  so if
  550.         // we're attempting to replicate its results,  we've gotta
  551.         // recompute everything,  too.
  552.         if ((FastMath.abs(savtsn - t) >= 30.0) || isDundeeCompliant)  {

  554.             savtsn = t;

  556.             // Update solar perturbations for time T
  557.             double zm = zmos + ZNS * t;
  558.             double zf = zm + 2 * ZES * FastMath.sin(zm);
  559.             double sinzf = FastMath.sin(zf);
  560.             double f2 = 0.5 * sinzf * sinzf - 0.25;
  561.             double f3 = -0.5 * sinzf * FastMath.cos(zf);
  562.             final double ses = se2 * f2 + se3 * f3;
  563.             final double sis = si2 * f2 + si3 * f3;
  564.             final double sls = sl2 * f2 + sl3 * f3 + sl4 * sinzf;
  565.             final double sghs = sgh2 * f2 + sgh3 * f3 + sgh4 * sinzf;
  566.             final double shs = sh2 * f2 + sh3 * f3;

  568.             // Update lunar perturbations for time T
  569.             zm = zmol + ZNL * t;
  570.             zf = zm + 2 * ZEL * FastMath.sin(zm);
  571.             sinzf = FastMath.sin(zf);
  572.             f2 =  0.5 * sinzf * sinzf - 0.25;
  573.             f3 = -0.5 * sinzf * FastMath.cos(zf);
  574.             final double sel = ee2 * f2 + e3 * f3;
  575.             final double sil = xi2 * f2 + xi3 * f3;
  576.             final double sll = xl2 * f2 + xl3 * f3 + xl4 * sinzf;
  577.             final double sghl = xgh2 * f2 + xgh3 * f3 + xgh4 * sinzf;
  578.             final double sh1 = xh2 * f2 + xh3 * f3;

  580.             // Sum the solar and lunar contributions
  581.             pe   = ses  + sel;
  582.             pinc = sis  + sil;
  583.             pl   = sls  + sll;
  584.             pgh  = sghs + sghl;
  585.             ph   = shs  + sh1;
  586.         }

  588.         xinc += pinc;

  590.         final double sinis = FastMath.sin( xinc);
  591.         final double cosis = FastMath.cos( xinc);

  593.         /* Add solar/lunar perturbation correction to eccentricity: */
  594.         em     += pe;
  595.         xll    += pl;
  596.         omgadf += pgh;
  597.         xinc    = MathUtils.normalizeAngle(xinc, 0);

  599.         if (FastMath.abs(xinc) >= 0.2) {
  600.             // Apply periodics directly
  601.             final double temp_val = ph / sinis;
  602.             omgadf -= cosis * temp_val;
  603.             xnode += temp_val;
  604.         } else {
  605.             // Apply periodics with Lyddane modification
  606.             final double sinok = FastMath.sin(xnode);
  607.             final double cosok = FastMath.cos(xnode);
  608.             final double alfdp =  ph * cosok + (pinc * cosis + sinis) * sinok;
  609.             final double betdp = -ph * sinok + (pinc * cosis + sinis) * cosok;
  610.             final double delta_xnode = MathUtils.normalizeAngle(FastMath.atan2(alfdp, betdp) - xnode, 0);
  611.             final double dls = -xnode * sinis * pinc;
  612.             omgadf += dls - cosis * delta_xnode;
  613.             xnode  += delta_xnode;
  614.         }
  615.     }

  617.     /** Computes internal secular derivs. */
  618.     private void computeSecularDerivs() {

  620.         final double sin_li = FastMath.sin(xli);
  621.         final double cos_li = FastMath.cos(xli);
  622.         final double sin_2li = 2. * sin_li * cos_li;
  623.         final double cos_2li = 2. * cos_li * cos_li - 1.;

  625.         // Dot terms calculated :
  626.         if (synchronous)  {
  627.             final double sin_3li = sin_2li * cos_li + cos_2li * sin_li;
  628.             final double cos_3li = cos_2li * cos_li - sin_2li * sin_li;
  629.             double term1a = del1 * (sin_li  * C_FASX2  - cos_li  * S_FASX2);
  630.             double term2a = del2 * (sin_2li * C_2FASX4 - cos_2li * S_2FASX4);
  631.             double term3a = del3 * (sin_3li * C_3FASX6 - cos_3li * S_3FASX6);
  632.             double term1b = del1 * (cos_li  * C_FASX2  + sin_li  * S_FASX2);
  633.             double term2b = 2.0 * del2 * (cos_2li * C_2FASX4 + sin_2li * S_2FASX4);
  634.             double term3b = 3.0 * del3 * (cos_3li * C_3FASX6 + sin_3li * S_3FASX6);

  636.             for (int j = 0; j < SECULAR_INTEGRATION_ORDER; j += 2)  {
  637.                 derivs[j]     = term1a + term2a + term3a;
  638.                 derivs[j + 1] = term1b + term2b + term3b;
  639.                 if ((i + 2) < SECULAR_INTEGRATION_ORDER) {
  640.                     term1a  = -term1a;
  641.                     term2a *= -4.0;
  642.                     term3a *= -9.0;
  643.                     term1b = -term1b;
  644.                     term2b *= -4.0;
  645.                     term3b *= -9.0;
  646.                 }
  647.             }
  648.         } else {
  649.             // orbit is a 12-hour resonant one
  650.             final double xomi = omegaq + omgdot * atime;
  651.             final double sin_omi = FastMath.sin(xomi);
  652.             final double cos_omi = FastMath.cos(xomi);
  653.             final double sin_li_m_omi = sin_li * cos_omi - sin_omi * cos_li;
  654.             final double sin_li_p_omi = sin_li * cos_omi + sin_omi * cos_li;
  655.             final double cos_li_m_omi = cos_li * cos_omi + sin_omi * sin_li;
  656.             final double cos_li_p_omi = cos_li * cos_omi - sin_omi * sin_li;
  657.             final double sin_2omi = 2. * sin_omi * cos_omi;
  658.             final double cos_2omi = 2. * cos_omi * cos_omi - 1.;
  659.             final double sin_2li_m_omi = sin_2li * cos_omi - sin_omi * cos_2li;
  660.             final double sin_2li_p_omi = sin_2li * cos_omi + sin_omi * cos_2li;
  661.             final double cos_2li_m_omi = cos_2li * cos_omi + sin_omi * sin_2li;
  662.             final double cos_2li_p_omi = cos_2li * cos_omi - sin_omi * sin_2li;
  663.             final double sin_2li_p_2omi = sin_2li * cos_2omi + sin_2omi * cos_2li;
  664.             final double cos_2li_p_2omi = cos_2li * cos_2omi - sin_2omi * sin_2li;
  665.             final double sin_2omi_p_li = sin_li * cos_2omi + sin_2omi * cos_li;
  666.             final double cos_2omi_p_li = cos_li * cos_2omi - sin_2omi * sin_li;
  667.             double term1a = d2201 * (sin_2omi_p_li * C_G22 - cos_2omi_p_li * S_G22) +
  668.                             d2211 * (sin_li * C_G22 - cos_li * S_G22) +
  669.                             d3210 * (sin_li_p_omi * C_G32 - cos_li_p_omi * S_G32) +
  670.                             d3222 * (sin_li_m_omi * C_G32 - cos_li_m_omi * S_G32) +
  671.                             d5220 * (sin_li_p_omi * C_G52 - cos_li_p_omi * S_G52) +
  672.                             d5232 * (sin_li_m_omi * C_G52 - cos_li_m_omi * S_G52);
  673.             double term2a = d4410 * (sin_2li_p_2omi * C_G44 - cos_2li_p_2omi * S_G44) +
  674.                             d4422 * (sin_2li * C_G44 - cos_2li * S_G44) +
  675.                             d5421 * (sin_2li_p_omi * C_G54 - cos_2li_p_omi * S_G54) +
  676.                             d5433 * (sin_2li_m_omi * C_G54 - cos_2li_m_omi * S_G54);
  677.             double term1b = d2201 * (cos_2omi_p_li * C_G22 + sin_2omi_p_li * S_G22) +
  678.                             d2211 * (cos_li * C_G22 + sin_li * S_G22) +
  679.                             d3210 * (cos_li_p_omi * C_G32 + sin_li_p_omi * S_G32) +
  680.                             d3222 * (cos_li_m_omi * C_G32 + sin_li_m_omi * S_G32) +
  681.                             d5220 * (cos_li_p_omi * C_G52 + sin_li_p_omi * S_G52) +
  682.                             d5232 * (cos_li_m_omi * C_G52 + sin_li_m_omi * S_G52);
  683.             double term2b = 2.0 * (d4410 * (cos_2li_p_2omi * C_G44 + sin_2li_p_2omi * S_G44) +
  684.                                    d4422 * (cos_2li * C_G44 + sin_2li * S_G44) +
  685.                                    d5421 * (cos_2li_p_omi * C_G54 + sin_2li_p_omi * S_G54) +
  686.                                    d5433 * (cos_2li_m_omi * C_G54 + sin_2li_m_omi * S_G54));

  688.             for (int j = 0; j < SECULAR_INTEGRATION_ORDER; j += 2) {
  689.                 derivs[j]     = term1a + term2a;
  690.                 derivs[j + 1] = term1b + term2b;
  691.                 if ((j + 2) < SECULAR_INTEGRATION_ORDER)  {
  692.                     term1a  = -term1a;
  693.                     term2a *= -4.0;
  694.                     term1b  = -term1b;
  695.                     term2b *= -4.0;
  696.                 }
  697.             }
  698.         }
  699.     }

  701. }
Created with JaCoCo 0.6.3.201306030806