Changes to the 'dpsec' part of the SDP4/SDP8 code

Last updated 22 April 2003

This describes some changes made in the satellite ephemeris C/C++ source code posted on this site. Specifically, I've made changes in two different areas of the Deep_dpsec() function in deep.cpp. 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.

However, this was even more likely to happen with the older implementation, which would only reset initial conditions when t changed sign. So the situation is at least better than it was, and the actual integration error does not usually appear to be a concern (though I still want to quantify that statement.) If it bothers you, you can revise the code very slightly to restore the initial conditions whenever the "new" time has a smaller absolute value than the "old" time. That is, change

`   if( fabs( deep_arg->t) < fabs( deep_arg->t - deep_arg->atime))`

`   if( fabs( deep_arg->t) < fabs( deep_arg->atime))`

But personally, I wouldn't bother; I think you'd be making the code slower for a very minor improvement in precision, and (probably) no real improvement at all in accuracy.

Also, instead of using a fixed 720-minute step size, the integrator now divides the integration range into equal steps of 720 minutes or less. As a result, we integrate exactly to the time t, instead of to a time that may be up to 720 minutes away from t. The original implementation got around this "close, but not exact" problem by doing a final partial step of size ft (lines 119 and 120). The new implementation doesn't have to do this.

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.)

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;
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;
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: }
```