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

to read:

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