#include <math.h>
#include "norad.h"
#include "norad_in.h"

#define zns      1.19459E-5
#define zes      0.01675
#define znl      1.5835218E-4
#define zel      0.05490
#define thdt     4.3752691E-3

      /* Previously,  the integration step was given as two variables:      */
      /* 'stepp' (positive step = +720) and 'stepn' (negative step = -720). */
      /* Exactly why this should be made a variable,  much less _different_ */
      /* variables for positive and negative,  is entirely unclear...       */
      /* (8 Apr 2003) INTEGRATION_STEP is now a maximum integration step.   */
      /* The code in 'dpsec' splits the integration range into equally-sized */
      /* pieces of 720 minutes (half a day) or smaller.                      */
#define INTEGRATION_STEP 720.
#define STEP2 (INTEGRATION_STEP * INTEGRATION_STEP / 2.)

static double ThetaG( double jd);

/* DEEP */
void Deep_dpinit( const tle_t *tle, deep_arg_t *deep_arg)
{
   const double sinq = sin(tle->xnodeo);
   const double cosq = cos(tle->xnodeo);
   const double aqnv = 1/deep_arg->aodp;
   const double c1ss   =  2.9864797E-6;
   const double day = tle->epoch - 2415020.;
                       /* days since 1900 Jan 0.5 = JD 2415020. */
   double zcosi =  0.91744867;
   double zsini =  0.39785416;
   double zsing = -0.98088458;
   double zcosg =  0.1945905;
   double bfact, cc = c1ss, se;
   double ze = zes, zn = zns;
   double sgh, sh, si;
   double zsinh = sinq, zcosh = cosq;
   double sl;
   int iteration;

   deep_arg->thgr = ThetaG( tle->epoch);
   deep_arg->xnq = deep_arg->xnodp;
   deep_arg->xqncl = tle->xincl;
   deep_arg->omegaq = tle->omegao;

   /* If the epoch has changed,  recompute (or initialize) the lunar and */
   /* solar terms.  I (BJG) am tempted to rewrite this so that they're   */
   /* _always_ recomputed,  since it appears that no two TLEs ever have  */
   /* the same epoch anyway;  all this careful effort to preserve data   */
   /* for the next pass is wasted.                                       */

   if (day != deep_arg->preep)
      {
      const double xnodce = 4.5236020 - 9.2422029E-4 * day;
      const double stem = sin(xnodce);
      const double ctem = cos(xnodce);
      const double c_minus_gam = 0.228027132 * day - 1.1151842;
      const double gam = 5.8351514 + 0.0019443680 * day;
      double zx, zy;

      deep_arg->preep = day;
      deep_arg->zcosil = 0.91375164 - 0.03568096 * ctem;
      deep_arg->zsinil = sqrt(1. - deep_arg->zcosil * deep_arg->zcosil);
      deep_arg->zsinhl = 0.089683511 * stem / deep_arg->zsinil;
      deep_arg->zcoshl = sqrt(1. - deep_arg->zsinhl*deep_arg->zsinhl);
      deep_arg->zmol = FMod2p( c_minus_gam);
      zx = 0.39785416 * stem / deep_arg->zsinil;
      zy = deep_arg->zcoshl * ctem + 0.91744867 * deep_arg->zsinhl * stem;
      zx = atan2( zx, zy) + gam - xnodce;
      deep_arg->zcosgl = cos( zx);
      deep_arg->zsingl = sin( zx);
      deep_arg->zmos = FMod2p( 6.2565837 + 0.017201977 * day);
      } /* End if(day != deep_arg->preep) */

   /* Do solar terms */
   deep_arg->savtsn = 1E20;

   /* 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( iteration = 0; iteration < 2; iteration++)
      {
      const double c1l = 4.7968065E-7;
      const double a1 = zcosg * zcosh + zsing * zcosi * zsinh;
      const double a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
      const double a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
      const double a8 = zsing * zsini;
      const double a9 = zsing * zsinh + zcosg * zcosi * zcosh;
      const double a10 = zcosg * zsini;
      const double a2 = deep_arg->cosio * a7 + deep_arg->sinio * a8;
      const double a4 = deep_arg->cosio * a9 + deep_arg->sinio * a10;
      const double a5 = -deep_arg->sinio * a7 + deep_arg->cosio * a8;
      const double a6 = -deep_arg->sinio * a9 + deep_arg->cosio * a10;
      const double x1 = a1 * deep_arg->cosg + a2 * deep_arg->sing;
      const double x2 = a3 * deep_arg->cosg + a4 * deep_arg->sing;
      const double x3 = -a1 * deep_arg->sing + a2 * deep_arg->cosg;
      const double x4 = -a3 * deep_arg->sing + a4 * deep_arg->cosg;
      const double x5 = a5 * deep_arg->sing;
      const double x6 = a6 * deep_arg->sing;
      const double x7 = a5 * deep_arg->cosg;
      const double x8 = a6 * deep_arg->cosg;
      const double z31 = 12 * x1 * x1 - 3 * x3 * x3;
      const double z32 = 24 * x1 * x2 - 6 * x3 * x4;
      const double z33 = 12 * x2 * x2 - 3 * x4 * x4;
      const double z11 = -6 * a1 * a5 + deep_arg->eosq * (-24 * x1 * x7 - 6 * x3 * x5);
      const double z12 = -6 * (a1 * a6 + a3 * a5) +  deep_arg->eosq *
                (-24 * (x2 * x7 + x1 * x8) - 6 * (x3 * x6 + x4 * x5));
      const double z13 = -6 * a3 * a6 + deep_arg->eosq * (-24 * x2 * x8 - 6 * x4 * x6);
      const double z21 = 6 * a2 * a5 + deep_arg->eosq * (24 * x1 * x5 - 6 * x3 * x7);
      const double z22 = 6 * (a4 * a5 + a2 * a6) +  deep_arg->eosq *
                (24 * (x2 * x5 + x1 * x6) - 6 * (x4 * x7 + x3 * x8));
      const double z23 = 6 * a4 * a6 + deep_arg->eosq * (24 * x2 * x6 - 6 * x4 * x8);
      const double s3 = cc / deep_arg->xnq;
      const double s2 = -0.5 * s3 / deep_arg->betao;
      const double s4 = s3 * deep_arg->betao;
      const double s1 = -15 * tle->eo * s4;
      const double s5 = x1 * x3 + x2 * x4;
      const double s6 = x2 * x3 + x1 * x4;
      const double s7 = x2 * x4 - x1 * x3;
      double z1 = 3 * (a1 * a1 + a2 * a2) + z31 * deep_arg->eosq;
      double z2 = 6 * (a1 * a3 + a2 * a4) + z32 * deep_arg->eosq;
      double z3 = 3 * (a3 * a3 + a4 * a4) + z33 * deep_arg->eosq;

      z1 = z1 + z1 + deep_arg->betao2 * z31;
      z2 = z2 + z2 + deep_arg->betao2 * z32;
      z3 = z3 + z3 + deep_arg->betao2 * z33;
      se = s1*zn*s5;
      si = s2*zn*(z11+z13);
      sl = -zn*s3*(z1+z3-14-6*deep_arg->eosq);
      sgh = s4*zn*(z31+z33-6);
      if (deep_arg->xqncl < 5.2359877E-2)
         sh = 0;
      else
         sh = -zn*s2*(z21+z23);
      deep_arg->ee2 = 2*s1*s6;
      deep_arg->e3 = 2*s1*s7;
      deep_arg->xi2 = 2*s2*z12;
      deep_arg->xi3 = 2*s2*(z13-z11);
      deep_arg->xl2 = -2*s3*z2;
      deep_arg->xl3 = -2*s3*(z3-z1);
      deep_arg->xl4 = -2*s3*(-21-9*deep_arg->eosq)*ze;
      deep_arg->xgh2 = 2*s4*z32;
      deep_arg->xgh3 = 2*s4*(z33-z31);
      deep_arg->xgh4 = -18*s4*ze;
      deep_arg->xh2 = -2*s2*z22;
      deep_arg->xh3 = -2*s2*(z23-z21);

      if( !iteration)   /* we compute lunar terms only on the first pass: */
         {
         deep_arg->sse = se;
         deep_arg->ssi = si;
         deep_arg->ssl = sl;
         deep_arg->ssh = sh/deep_arg->sinio;
         deep_arg->ssg = sgh-deep_arg->cosio*deep_arg->ssh;
         deep_arg->se2 = deep_arg->ee2;
         deep_arg->si2 = deep_arg->xi2;
         deep_arg->sl2 = deep_arg->xl2;
         deep_arg->sgh2 = deep_arg->xgh2;
         deep_arg->sh2 = deep_arg->xh2;
         deep_arg->se3 = deep_arg->e3;
         deep_arg->si3 = deep_arg->xi3;
         deep_arg->sl3 = deep_arg->xl3;
         deep_arg->sgh3 = deep_arg->xgh3;
         deep_arg->sh3 = deep_arg->xh3;
         deep_arg->sl4 = deep_arg->xl4;
         deep_arg->sgh4 = deep_arg->xgh4;
         zcosg = deep_arg->zcosgl;
         zsing = deep_arg->zsingl;
         zcosi = deep_arg->zcosil;
         zsini = deep_arg->zsinil;
         zcosh = deep_arg->zcoshl*cosq+deep_arg->zsinhl*sinq;
         zsinh = sinq*deep_arg->zcoshl-cosq*deep_arg->zsinhl;
         zn = znl;
         cc = c1l;
         ze = zel;
         }
      }

   deep_arg->sse += se;
   deep_arg->ssi += si;
   deep_arg->ssl += sl;
   deep_arg->ssg += sgh - deep_arg->cosio / deep_arg->sinio * sh;
   deep_arg->ssh += sh / deep_arg->sinio;

   if( deep_arg->xnq >= 0.00826 && deep_arg->xnq <= 0.00924 && tle->eo >= .5)
      {           /* start of 12-hour orbit, e >.5 section */
      const double root22 = 1.7891679E-6;
      const double root32 = 3.7393792E-7;
      const double root44 = 7.3636953E-9;
      const double root52 = 1.1428639E-7;
      const double root54 = 2.1765803E-9;
      const double g201 = -0.306 - (tle->eo - 0.64) * 0.440;
      const double eoc = tle->eo * deep_arg->eosq;
      const double sini2 = deep_arg->sinio*deep_arg->sinio;
      const double f220 = 0.75*(1+2*deep_arg->cosio+deep_arg->theta2);
      const double f221 = 1.5 * sini2;
      const double f321 = 1.875 * deep_arg->sinio * (1 - 2 *\
               deep_arg->cosio - 3 * deep_arg->theta2);
      const double f322 = -1.875*deep_arg->sinio*(1+2*
               deep_arg->cosio-3*deep_arg->theta2);
      const double f441 = 35 * sini2 * f220;
      const double f442 = 39.3750 * sini2 * sini2;
      const double f522 = 9.84375*deep_arg->sinio*(sini2*(1-2*deep_arg->cosio-5*
                 deep_arg->theta2)+0.33333333*(-2+4*deep_arg->cosio+
                 6*deep_arg->theta2));
      const double f523 = deep_arg->sinio*(4.92187512*sini2*(-2-4*
                 deep_arg->cosio+10*deep_arg->theta2)+6.56250012
                 *(1+2*deep_arg->cosio-3*deep_arg->theta2));
      const double f542 = 29.53125*deep_arg->sinio*(2-8*
                 deep_arg->cosio+deep_arg->theta2*
                 (-12+8*deep_arg->cosio+10*deep_arg->theta2));
      const double f543 = 29.53125*deep_arg->sinio*(-2-8*deep_arg->cosio+
                 deep_arg->theta2*(12+8*deep_arg->cosio-10*
                 deep_arg->theta2));
      double g410, g422, g520, g521, g532, g533;
      double g211, g310, g322;
      double temp, temp1;

      deep_arg->resonance_flag = 1;       /* it _is_ resonant... */
      deep_arg->synchronous_flag = 0;     /* but it's not synchronous */
             /* Geopotential resonance initialization for 12 hour orbits: */
      if (tle->eo <= 0.65)
         {
         g211 = 3.616-13.247*tle->eo+16.290*deep_arg->eosq;
         g310 = -19.302+117.390*tle->eo-228.419*
                    deep_arg->eosq+156.591*eoc;
         g322 = -18.9068+109.7927*tle->eo-214.6334*
                    deep_arg->eosq+146.5816*eoc;
         g410 = -41.122+242.694*tle->eo-471.094*
                    deep_arg->eosq+313.953*eoc;
         g422 = -146.407+841.880*tle->eo-1629.014*
                    deep_arg->eosq+1083.435*eoc;
         g520 = -532.114+3017.977*tle->eo-5740*
                    deep_arg->eosq+3708.276*eoc;
         }
      else
         {
         g211 = -72.099+331.819*tle->eo-508.738*
                    deep_arg->eosq+266.724*eoc;
         g310 = -346.844+1582.851*tle->eo-2415.925*
                    deep_arg->eosq+1246.113*eoc;
         g322 = -342.585+1554.908*tle->eo-2366.899*
                    deep_arg->eosq+1215.972*eoc;
         g410 = -1052.797+4758.686*tle->eo-7193.992*
                    deep_arg->eosq+3651.957*eoc;
         g422 = -3581.69+16178.11*tle->eo-24462.77*
                    deep_arg->eosq+ 12422.52*eoc;
         if (tle->eo <= 0.715)
            g520 = 1464.74-4664.75*tle->eo+3763.64*deep_arg->eosq;
         else
            g520 = -5149.66+29936.92*tle->eo-54087.36*
                         deep_arg->eosq+31324.56*eoc;
         } /* End if (tle->eo <= 0.65) */

      if (tle->eo < 0.7)
         {
         g533 = -919.2277+4988.61*tle->eo-9064.77*
                    deep_arg->eosq+5542.21*eoc;
         g521 = -822.71072+4568.6173*tle->eo-8491.4146*
                    deep_arg->eosq+5337.524*eoc;
         g532 = -853.666+4690.25*tle->eo-8624.77*
                    deep_arg->eosq+ 5341.4*eoc;
         }
      else
         {
         g533 = -37995.78+161616.52*tle->eo-229838.2*
                    deep_arg->eosq+109377.94*eoc;
         g521 = -51752.104+218913.95*tle->eo-309468.16*
                    deep_arg->eosq+146349.42*eoc;
         g532 = -40023.88+170470.89*tle->eo-242699.48*
                    deep_arg->eosq+115605.82*eoc;
         } /* End if (tle->eo <= 0.7) */

      temp1 = 3 * deep_arg->xnq * deep_arg->xnq * aqnv * aqnv;
      temp = temp1*root22;
      deep_arg->d2201 = temp * f220 * g201;
      deep_arg->d2211 = temp * f221 * g211;
      temp1 *= aqnv;
      temp = temp1*root32;
      deep_arg->d3210 = temp * f321 * g310;
      deep_arg->d3222 = temp * f322 * g322;
      temp1 *= aqnv;
      temp = 2*temp1*root44;
      deep_arg->d4410 = temp * f441 * g410;
      deep_arg->d4422 = temp * f442 * g422;
      temp1 *= aqnv;
      temp = temp1*root52;
      deep_arg->d5220 = temp * f522 * g520;
      deep_arg->d5232 = temp * f523 * g532;
      temp = 2*temp1*root54;
      deep_arg->d5421 = temp * f542 * g521;
      deep_arg->d5433 = temp * f543 * g533;
      deep_arg->xlamo = tle->xmo+tle->xnodeo+tle->xnodeo-deep_arg->thgr-deep_arg->thgr;
      bfact = deep_arg->xmdot + deep_arg->xnodot+
                   deep_arg->xnodot - thdt - thdt;
      bfact += deep_arg->ssl + deep_arg->ssh + deep_arg->ssh;
      }        /* end of 12-hour orbit, e >.5 section */
   else if( (deep_arg->xnq < 0.0052359877) && (deep_arg->xnq > 0.0034906585))
      {
      const double q22    =  1.7891679E-6;
      const double q31    =  2.1460748E-6;
      const double q33    =  2.2123015E-7;
      const double cosio_plus_1 = 1. + deep_arg->cosio;
      const double g200 = 1+deep_arg->eosq*(-2.5+0.8125*deep_arg->eosq);
      const double g300 = 1+deep_arg->eosq*(-6+6.60937*deep_arg->eosq);
      const double f311 = 0.9375*deep_arg->sinio*deep_arg->sinio*
             (1+3*deep_arg->cosio)-0.75*cosio_plus_1;
      const double g310 = 1+2*deep_arg->eosq;
      const double f220 = 0.75 * cosio_plus_1 * cosio_plus_1;
      const double f330 = 2.5 * f220 * cosio_plus_1;

      deep_arg->resonance_flag = deep_arg->synchronous_flag = 1;
      /* Synchronous resonance terms initialization */
      deep_arg->del1 = 3*deep_arg->xnq*deep_arg->xnq*aqnv*aqnv;
      deep_arg->del2 = 2*deep_arg->del1*f220*g200*q22;
      deep_arg->del3 = 3*deep_arg->del1*f330*g300*q33*aqnv;
      deep_arg->del1 = deep_arg->del1*f311*g310*q31*aqnv;
      deep_arg->xlamo = tle->xmo+tle->xnodeo+tle->omegao-deep_arg->thgr;
      bfact = deep_arg->xmdot + deep_arg->omgdot + deep_arg->xnodot - thdt;
      bfact = bfact+deep_arg->ssl+deep_arg->ssg+deep_arg->ssh;
      } /* End of geosych case */
   else              /* it's neither a high-e 12-hr orbit nor a geosynch: */
      deep_arg->resonance_flag = deep_arg->synchronous_flag = 0;

   if( deep_arg->resonance_flag)
      {
      deep_arg->xfact = bfact-deep_arg->xnq;

      /* Initialize integrator */
      deep_arg->xli = deep_arg->xlamo;
      deep_arg->xni = deep_arg->xnq;
      deep_arg->atime = 0;
      }
   /* End case dpinit: */
}

void Deep_dpsec( const tle_t *tle, deep_arg_t *deep_arg)
{
   double delt, temp;
   int n_steps;

   deep_arg->xll += deep_arg->ssl*deep_arg->t;
   deep_arg->omgadf += deep_arg->ssg*deep_arg->t;
   deep_arg->xnode += deep_arg->ssh*deep_arg->t;
   deep_arg->em = tle->eo+deep_arg->sse*deep_arg->t;
   deep_arg->xinc = tle->xincl+deep_arg->ssi*deep_arg->t;
   if (deep_arg->xinc < 0)       /* Begin April 1983 errata correction: */
      {
      deep_arg->xinc = -deep_arg->xinc;
      deep_arg->xnode = deep_arg->xnode + pi;
      deep_arg->omgadf = deep_arg->omgadf-pi;
      }                          /* End April 1983 errata correction. */
   if( !deep_arg->resonance_flag ) return;

            /* 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 */
   if( fabs( deep_arg->t) < fabs( deep_arg->t - deep_arg->atime))
      {                                    /* Epoch restart */
      deep_arg->atime = 0;
      deep_arg->xni = deep_arg->xnq;
      deep_arg->xli = deep_arg->xlamo;
      }

               /* How many integration steps does it take to get   */
               /* get from our starting time, deep_arg->atime,     */
               /* to the desired time,  deep_arg->t?:              */
   n_steps = (int)ceil( fabs( deep_arg->t - deep_arg->atime) / INTEGRATION_STEP);
   if( n_steps)
      delt = (deep_arg->t - deep_arg->atime) / (double)n_steps;
   else
      delt = 0.;

   while( n_steps--)
      {
      const double sin_li = sin( deep_arg->xli);
      const double cos_li = cos( deep_arg->xli);
      const double sin_2li = 2. * sin_li * cos_li;
      const double cos_2li = 2. * cos_li * cos_li - 1.;
      double xldot, xnddt, xndot;

                /* Dot terms calculated,  using a lot of trig add/subtract */
                /* identities to reduce the computational load... at the   */
                /* cost of making the code somewhat hard to follow:        */
      if( deep_arg->synchronous_flag )
         {
/*       const double fasx2 = 0.13130908;       */
/*       const double fasx4 = 2.8843198;        */
/*       const double fasx6 = 0.37448087;       */
         const double c_fasx2 =  0.99139134268488593;
         const double s_fasx2 =  0.13093206501640101;
         const double c_2fasx4 =  0.87051638752972937;
         const double s_2fasx4 = -0.49213943048915526;
         const double c_3fasx6 =  0.43258117585763334;
         const double s_3fasx6 =  0.90159499016666422;
         const double sin_3li = sin_2li * cos_li + cos_2li * sin_li;
         const double cos_3li = cos_2li * cos_li - sin_2li * sin_li;

         xndot = deep_arg->del1 * (sin_li  * c_fasx2  - cos_li  * s_fasx2)
               + deep_arg->del2 * (sin_2li * c_2fasx4 - cos_2li * s_2fasx4)
               + deep_arg->del3 * (sin_3li * c_3fasx6 - cos_3li * s_3fasx6);
         xnddt = deep_arg->del1 * (cos_li  * c_fasx2  + sin_li  * s_fasx2)
           + 2 * deep_arg->del2 * (cos_2li * c_2fasx4 + sin_2li * s_2fasx4)
           + 3 * deep_arg->del3 * (cos_3li * c_3fasx6 + sin_3li * s_3fasx6);
         }        /* end of geosynch case */
      else
         {        /* orbit is a 12-hour resonant one: */
/*       const double g22    =  5.7686396;      */
/*       const double g32    =  0.95240898;     */
/*       const double g44    =  1.8014998;      */
/*       const double g52    =  1.0508330;      */
/*       const double g54    =  4.4108898;      */
         const double c_g22 =  0.87051638752972937;
         const double s_g22 = -0.49213943048915526;
         const double c_g32 =  0.57972190187001149;
         const double s_g32 =  0.81481440616389245;
         const double c_g44 = -0.22866241528815548;
         const double s_g44 =  0.97350577801807991;
         const double c_g52 =  0.49684831179884198;
         const double s_g52 =  0.86783740128127729;
         const double c_g54 = -0.29695209575316894;
         const double s_g54 = -0.95489237761529999;
         const double xomi =
                   deep_arg->omegaq + deep_arg->omgdot * deep_arg->atime;
         const double sin_omi = sin( xomi), cos_omi = cos( xomi);
         const double sin_li_m_omi = sin_li * cos_omi - sin_omi * cos_li;
         const double sin_li_p_omi = sin_li * cos_omi + sin_omi * cos_li;
         const double cos_li_m_omi = cos_li * cos_omi + sin_omi * sin_li;
         const double cos_li_p_omi = cos_li * cos_omi - sin_omi * sin_li;
         const double sin_2omi = 2. * sin_omi * cos_omi;
         const double cos_2omi = 2. * cos_omi * cos_omi - 1.;
         const double sin_2li_m_omi = sin_2li * cos_omi - sin_omi * cos_2li;
         const double sin_2li_p_omi = sin_2li * cos_omi + sin_omi * cos_2li;
         const double cos_2li_m_omi = cos_2li * cos_omi + sin_omi * sin_2li;
         const double cos_2li_p_omi = cos_2li * cos_omi - sin_omi * sin_2li;
         const double sin_2li_p_2omi = sin_2li * cos_2omi + sin_2omi * cos_2li;
         const double cos_2li_p_2omi = cos_2li * cos_2omi - sin_2omi * sin_2li;
         const double sin_2omi_p_li = sin_li * cos_2omi + sin_2omi * cos_li;
         const double cos_2omi_p_li = cos_li * cos_2omi - sin_2omi * sin_li;

         xndot = deep_arg->d2201 * (sin_2omi_p_li*c_g22 - cos_2omi_p_li*s_g22)
               + deep_arg->d2211 * (sin_li * c_g22 - cos_li * s_g22)
               + deep_arg->d3210 * (sin_li_p_omi*c_g32 - cos_li_p_omi*s_g32)
               + deep_arg->d3222 * (sin_li_m_omi*c_g32 - cos_li_m_omi*s_g32)
               + deep_arg->d4410 * (sin_2li_p_2omi*c_g44 - cos_2li_p_2omi*s_g44)
               + deep_arg->d4422 * (sin_2li * c_g44 - cos_2li * s_g44)
               + deep_arg->d5220 * (sin_li_p_omi*c_g52 - cos_li_p_omi*s_g52)
               + deep_arg->d5232 * (sin_li_m_omi*c_g52 - cos_li_m_omi*s_g52)
               + deep_arg->d5421 * (sin_2li_p_omi*c_g54 - cos_2li_p_omi*s_g54)
               + deep_arg->d5433 * (sin_2li_m_omi*c_g54 - cos_2li_m_omi*s_g54);
         xnddt = deep_arg->d2201 * (cos_2omi_p_li*c_g22 + sin_2omi_p_li*s_g22)
               + deep_arg->d2211 * (cos_li * c_g22 + sin_li * s_g22)
               + deep_arg->d3210 * (cos_li_p_omi*c_g32 + sin_li_p_omi*s_g32)
               + deep_arg->d3222 * (cos_li_m_omi*c_g32 + sin_li_m_omi*s_g32)
               + deep_arg->d5220 * (cos_li_p_omi*c_g52 + sin_li_p_omi*s_g52)
               + deep_arg->d5232 * (cos_li_m_omi*c_g52 + sin_li_m_omi*s_g52)
            + 2*(deep_arg->d4410 * (cos_2li_p_2omi*c_g44 + sin_2li_p_2omi*s_g44)
               + deep_arg->d4422 * (cos_2li * c_g44 + sin_2li * s_g44)
               + deep_arg->d5421 * (cos_2li_p_omi*c_g54 + sin_2li_p_omi*s_g54)
               + deep_arg->d5433 * (cos_2li_m_omi*c_g54 + sin_2li_m_omi*s_g54));
         } /* End of 12-hr resonant case */

      xldot = deep_arg->xni+deep_arg->xfact;
      xnddt *= xldot;

      deep_arg->xli += delt * (xldot + xndot * delt / 2.);
      deep_arg->xni += delt * (xndot + xnddt * delt / 2.);
      deep_arg->atime += delt;
      }

   deep_arg->xn = deep_arg->xni;

   temp = -deep_arg->xnode + deep_arg->thgr + deep_arg->t * thdt;

   deep_arg->xll = deep_arg->xli + temp
         + (deep_arg->synchronous_flag ? -deep_arg->omgadf : temp);
   /*End case dpsec: */
}

void Deep_dpper( deep_arg_t *deep_arg)
{
               /* 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.    */
      if (fabs(deep_arg->savtsn-deep_arg->t) >= 30.)
         {
         double zf, zm, sinzf, ses, sis, sil, sel, sll, sls;
         double f2, f3, sghl, sghs, shs, sh1;

         deep_arg->savtsn = deep_arg->t;

               /* Update solar perturbations for time T: */
         zm = deep_arg->zmos+zns*deep_arg->t;
         zf = zm+2*zes*sin(zm);
         sinzf = sin(zf);
         f2 = 0.5*sinzf*sinzf-0.25;
         f3 = -0.5*sinzf*cos(zf);
         ses = deep_arg->se2*f2+deep_arg->se3*f3;
         sis = deep_arg->si2*f2+deep_arg->si3*f3;
         sls = deep_arg->sl2*f2+deep_arg->sl3*f3+deep_arg->sl4*sinzf;
         sghs = deep_arg->sgh2*f2+deep_arg->sgh3*f3+deep_arg->sgh4*sinzf;
         shs = deep_arg->sh2*f2+deep_arg->sh3*f3;

               /* Update lunar perturbations for time T: */
         zm = deep_arg->zmol+znl*deep_arg->t;
         zf = zm+2*zel*sin(zm);
         sinzf = sin(zf);
         f2 = 0.5*sinzf*sinzf-0.25;
         f3 = -0.5*sinzf*cos(zf);
         sel = deep_arg->ee2*f2+deep_arg->e3*f3;
         sil = deep_arg->xi2*f2+deep_arg->xi3*f3;
         sll = deep_arg->xl2*f2+deep_arg->xl3*f3+deep_arg->xl4*sinzf;
         sghl = deep_arg->xgh2*f2+deep_arg->xgh3*f3+deep_arg->xgh4*sinzf;
         sh1 = deep_arg->xh2*f2+deep_arg->xh3*f3;

               /* Sum the solar and lunar contributions: */
         deep_arg->pe = ses+sel;
         deep_arg->pinc = sis+sil;
         deep_arg->pl = sls+sll;
         deep_arg->pgh = sghs+sghl;
         deep_arg->ph = shs+sh1;
#ifdef RETAIN_PERTURBATION_VALUES_AT_EPOCH
         if( deep_arg->solar_lunar_init_flag)
            {
            deep_arg->pe0 = deep_arg->pe;
            deep_arg->pinc0 = deep_arg->pinc;
            deep_arg->pl0 = deep_arg->pl;
            deep_arg->pgh0 = deep_arg->pgh;
            deep_arg->ph0 = deep_arg->ph;
            }
         deep_arg->pe  -= deep_arg->pe0;
         deep_arg->pinc -= deep_arg->pinc0;
         deep_arg->pl  -= deep_arg->pl0;
         deep_arg->pgh -= deep_arg->pgh0;
         deep_arg->ph  -= deep_arg->ph0;
         if( deep_arg->solar_lunar_init_flag)
            return;        /* done all we really need to do here... */
#endif
         }

               /* In Spacetrack 3, sinis & cosis were initialized       */
               /* _before_ perturbations were added to xinc.  In        */
               /* Spacetrack 6,  it's the other way around (see below). */
#ifdef SPACETRACK_3
      const double sinis = sin( deep_arg->xinc);
      const double cosis = cos( deep_arg->xinc);
#endif
         /* Add solar/lunar perturbation correction to inclination: */
      deep_arg->xinc += deep_arg->pinc;

         /* Add solar/lunar perturbation correction to eccentricity: */
      deep_arg->em += deep_arg->pe;

      if (deep_arg->xqncl >= 0.2)
         {
                /* Apply periodics directly */
         const double temp_val = deep_arg->ph / deep_arg->sinio;

         deep_arg->omgadf += deep_arg->pgh - deep_arg->cosio * temp_val;
         deep_arg->xnode += temp_val;
         deep_arg->xll += deep_arg->pl;
         }
      else
         {
           /* Apply periodics with Lyddane modification */
         const double sinok = sin(deep_arg->xnode);
         const double cosok = cos(deep_arg->xnode);
               /* Correction from Spacetrack Report #3 to #6:      */
               /* used to be sinis & cosis were computed _before_  */
               /* adding perturbations to XINC.  Now it's _after_: */
#ifndef SPACETRACK_3
         const double sinis = sin(deep_arg->xinc);
         const double cosis = cos(deep_arg->xinc);
#endif
         const double alfdp = deep_arg->ph * cosok
                    + (deep_arg->pinc * cosis + sinis) * sinok;
         const double betdp = - deep_arg->ph * sinok
                    + (deep_arg->pinc * cosis + sinis) * cosok;
         double xls, xnoh;

         deep_arg->xnode = FMod2p(deep_arg->xnode);
         xls = deep_arg->xll + deep_arg->omgadf + cosis * deep_arg->xnode;
         xls += deep_arg->pl + deep_arg->pgh
                            - deep_arg->pinc * deep_arg->xnode * sinis;
         xnoh = deep_arg->xnode;
         deep_arg->xnode = atan2(alfdp,betdp);

          /* This is a patch to Lyddane modification suggested */
          /* by Rob Matson, streamlined very slightly by BJG, to */
          /* keep 'xnode' & 'xnoh' within 180 degrees of each other. */

         if( deep_arg->xnode < xnoh - pi)
            deep_arg->xnode += twopi;
         else if( deep_arg->xnode > xnoh + pi)
            deep_arg->xnode -= twopi;

         deep_arg->xll += deep_arg->pl;
         deep_arg->omgadf = xls - deep_arg->xll - cos(deep_arg->xinc)*
                             deep_arg->xnode;
         } /* End case dpper: */
}

/*------------------------------------------------------------------*/

/* THETAG */
static double ThetaG( double jd)
{
  /* Reference:  The 1992 Astronomical Almanac, page B6. */
  const double omega_E = 1.00273790934;
                   /* Earth rotations per sidereal day (non-constant) */
  const double UT = fmod( jd + .5, 1.);
  double t_cen, GMST, rval;

  t_cen = (jd - UT - 2451545.0) / 36525.;
  GMST = 24110.54841 + t_cen * (8640184.812866 + t_cen *
                           (0.093104 - t_cen * 6.2E-6));
  GMST = fmod( GMST + secday * omega_E * UT, secday);
  if( GMST < 0.)
     GMST += secday;
  rval = twopi * GMST / secday;

  return( rval);
} /*Function thetag*/