Changes to the dpsec part of the SDP4/SDP8 code

Last updated 2019 August 19

This describes some changes made to my satellite ephemeris C/C++ source code. Specifically, I've made changes to the Deep_dpsec() function in deep.cpp that may be not particularly clear (though still much clearer than the original code was). This code does a numerical integration to determine secular variations in the orbits of geosynchronous and high-eccentricity 12-hour-period objects. I've improved the implementation of the numerical integration, and (of somewhat lesser importance) improved its speed through the use of trigonometric identities.

How the numerical integration in Deep_dpsec() used to work: The "normal" implementation of this function has some bizarre aspects to it, resulting in more complex code and, in some cases, unnecessary slowdowns and increased integration error. Click here for the original implementation.

Figuring out how this code works took a while. Basically, secular perturbations xli and xni have to be numerically integrated. Their initial values for tsince=0 are computed as part of Deep_dpinit(), and are stored as xlamo and xnq, respectively.

The integration has to be done from a given point (initially t=0) to the "present". However, the results of that integration are saved between calls to Deep_dpsec(), so the integrator can pick up where it left off. The time of the previous call is stored in atime.

If you call the function first for a positive time, then a negative time, or vice versa, it would make sense to go back to the initial t=0 values. That saves some integrating and you get rid of the accumulated integration error. Indeed, lines 22-30 of the original code do just that.

The actual integration was done with a fixed step size of 720 minutes. Unless t changed sign, resulting in restoring the initial integrator state, integration was always done from the previous time atime to the current time. Thus, if you evaluated a satellite for times t=1000000 and t=50, the integrator would first take 1000000/720 steps forward , then take the same number of steps backward to reach t=50. Because of the fixed 720-minute step size, however, it would actually go back to t=0... so only at the end of the integration would it "wake up" and realize that it ought to use the initial values!

The revised code now looks to see if t is "closer" to t=0 or t=atime. If the former, it restores the initial conditions. If the latter, it uses the results of the previous integration. This can still result in some troubles if, for example, you invoked the function for t=100000, t=51000, t=100000, t=51000, etc. The code would integrate forward and backward, gradually accumulating integration error. You can then be puzzled as to why the code gave one result for t=100000 the first time, and different results the next time.

As a result, the Dundee SGP4/SDP4 code and the Vallado et. al. code restarts the integrator. Invoke the function for t=100000 and t=99000 seconds, and the second integration will restart instead of "going backwards" 1000 seconds. For a particular value of t, you always get repeatable results. My code doesn't default to doing this, but one can set Dundee-compliant mode to get that behavior.

Speedups through use of trigonometric identities: The actual quantities being integrated are six trig terms for geosynchs (lines 60 to 69 of the original code) or 20 trig terms for 12-hour orbits (lines 73 to 102). Doing this numerical integration is often the slowest part of computing satellite positions, so reducing that "per-step" workload is important to me.

The nature of the terms is such that they are very amenable to handling through trigonometric addition formulae:

sin(a+b) = sin(a) cos(b) + sin(b) cos(a)
cos(a+b) = cos(a) cos(b) + sin(b) sin(a) 

In the revised implementation, sines and cosines are computed twice per step for synchronous orbits, four times for resonant orbits. These are three- and fivefold reductions, respectively. (Unfortunately, this is one of the few cases where the "improvement" made the code less easy to understand.)

Derivatives computed in a separate function : The numerical integration works by computing the derivatives of xli and xni. I've broken that out into a separate compute_dpsec_derivs() function. In the process, I realized that computing higher-order derivatives for those variables would be just about trivial, and would allow for a larger step size. I did this, and it works, and speeds things up nicely. The only drawback is that it's not compliant with the way JSpOC and DoD do things, and you won't get exactly the same results they do. So by default, we just compute the first derivatives and stop there.

If backward compatibility wasn't an issue, we'd almost certainly compute at least the second derivative and do a full-fledged Taylor series integration. But backward compatibility is an issue.

Older (reference) Deep_dpsec():

001: void Deep_dpsec( const tle_t *tle, deep_arg_t *deep_arg)
002: {
003:    double xl, xldot, xnddt, xndot, temp;
004:    double delt = 0., ft = 0.;
005:    int do_loop_flag = 0, epoch_restart_flag = 0;
006:
007:    deep_arg->xll += deep_arg->ssl*deep_arg->t;
008:    deep_arg->omgadf += deep_arg->ssg*deep_arg->t;
009:    deep_arg->xnode += deep_arg->ssh*deep_arg->t;
010:    deep_arg->em = tle->eo+deep_arg->sse*deep_arg->t;
011:    deep_arg->xinc = tle->xincl+deep_arg->ssi*deep_arg->t;
012:    if (deep_arg->xinc < 0)       /* Begin April 1983 errata correction: */
013:       {
014:       deep_arg->xinc = -deep_arg->xinc;
015:       deep_arg->xnode = deep_arg->xnode + pi;
016:       deep_arg->omgadf = deep_arg->omgadf-pi;
017:       }                          /* End April 1983 errata correction. */
018:    if( !deep_arg->resonance_flag ) return;
019:
020:    do
021:       {
022:       if( (deep_arg->atime == 0) ||
023:             ((deep_arg->t >= 0) && (deep_arg->atime < 0)) ||
024:             ((deep_arg->t < 0) && (deep_arg->atime >= 0)) )
025:          {
026:             /* Epoch restart */
027:          delt = ((deep_arg->t >= 0) ? INTEGRATION_STEP : -INTEGRATION_STEP);
028:          deep_arg->atime = 0;
029:          deep_arg->xni = deep_arg->xnq;
030:          deep_arg->xli = deep_arg->xlamo;
031:          }
032:       else
033:          {
034:          if( fabs(deep_arg->t) >= fabs(deep_arg->atime) )
035:             delt = ((deep_arg->t > 0) ? INTEGRATION_STEP : -INTEGRATION_STEP);
036:          }
037:
038:       do
039:          {
040:          if ( fabs(deep_arg->t-deep_arg->atime) >= INTEGRATION_STEP )
041:             {
042:             do_loop_flag = 1;
043:             epoch_restart_flag = 0;
044:             }
045:          else
046:             {
047:             ft = deep_arg->t-deep_arg->atime;
048:             do_loop_flag = 0;
049:             }
050:
051:          if( fabs(deep_arg->t) < fabs(deep_arg->atime) )
052:             {
053:             delt = ((deep_arg->t >= 0) ? -INTEGRATION_STEP : INTEGRATION_STEP);
054:             do_loop_flag = epoch_restart_flag = 1;
055:             }
056:
057:          /* Dot terms calculated */
058:          if( deep_arg->synchronous_flag )
059:             {
060:             const double fasx2 = 0.13130908;
061:             const double fasx4 = 2.8843198;
062:             const double fasx6 = 0.37448087;
063:
064:             xndot = deep_arg->del1 * sin(     deep_arg->xli - fasx2)
065:                   + deep_arg->del2 * sin(2 * (deep_arg->xli - fasx4))
066:                   + deep_arg->del3 * sin(3 * (deep_arg->xli - fasx6));
067:             xnddt = deep_arg->del1 * cos(     deep_arg->xli - fasx2)
068:               + 2 * deep_arg->del2 * cos(2 * (deep_arg->xli - fasx4))
069:               + 3 * deep_arg->del3 * cos(3 * (deep_arg->xli - fasx6));
070:             }        /* end of geosynch case */
071:          else
072:             {        /* orbit is a 12-hour resonant one: */
073:             const double g22    =  5.7686396;
074:             const double g32    =  0.95240898;
075:             const double g44    =  1.8014998;
076:             const double g52    =  1.0508330;
077:             const double g54    =  4.4108898;
078:             const double xomi =
079:                       deep_arg->omegaq + deep_arg->omgdot * deep_arg->atime;
080:             const double x2omi = xomi + xomi;
081:             const double x2li = deep_arg->xli + deep_arg->xli;
082:
083:             xndot = deep_arg->d2201 * sin( x2omi+deep_arg->xli-g22)
084:                   + deep_arg->d2211 * sin( deep_arg->xli-g22)
085:                   + deep_arg->d3210 * sin(  xomi+deep_arg->xli-g32)
086:                   + deep_arg->d3222 * sin( -xomi+deep_arg->xli-g32)
087:                   + deep_arg->d4410 * sin( x2omi+x2li-g44)
088:                   + deep_arg->d4422 * sin( x2li-g44)
089:                   + deep_arg->d5220 * sin(  xomi+deep_arg->xli-g52)
090:                   + deep_arg->d5232 * sin( -xomi+deep_arg->xli-g52)
091:                   + deep_arg->d5421 * sin(  xomi+x2li-g54)
092:                   + deep_arg->d5433 * sin( -xomi+x2li-g54);
093:             xnddt = deep_arg->d2201 * cos( x2omi+deep_arg->xli-g22)
094:                   + deep_arg->d2211 * cos( deep_arg->xli-g22)
095:                   + deep_arg->d3210 * cos(  xomi+deep_arg->xli-g32)
096:                   + deep_arg->d3222 * cos( -xomi+deep_arg->xli-g32)
097:                   + deep_arg->d5220 * cos(  xomi+deep_arg->xli-g52)
098:                   + deep_arg->d5232 * cos( -xomi+deep_arg->xli-g52)
099:                   + 2*(deep_arg->d4410 * cos( x2omi+x2li-g44)
100:                   + deep_arg->d4422 * cos( x2li-g44)
101:                   + deep_arg->d5421 * cos(  xomi+x2li-g54)
102:                   + deep_arg->d5433 * cos( -xomi+x2li-g54));
103:             } /* End of 12-hr resonant case */
104:
105:          xldot = deep_arg->xni+deep_arg->xfact;
106:          xnddt *= xldot;
107:
108:          if( do_loop_flag)
109:             {
110:             deep_arg->xli += xldot * delt + xndot * STEP2;
111:             deep_arg->xni += xndot * delt + xnddt * STEP2;
112:             deep_arg->atime += delt;
113:             }
114:          }
115:       while( do_loop_flag && !epoch_restart_flag);
116:       }
117:    while( do_loop_flag && epoch_restart_flag);
118:
119:    deep_arg->xn = deep_arg->xni + ft * (xndot + xnddt * ft * 0.5);
120:    xl = deep_arg->xli + ft * (xldot + xndot * ft * 0.5);
121:    temp = -deep_arg->xnode + deep_arg->thgr + deep_arg->t * thdt;
122:
123:    deep_arg->xll = xl + temp
124:          + (deep_arg->synchronous_flag ? -deep_arg->omgadf : temp);
125:    /*End case dpsec: */
126: }