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