DEEP Subroutine

Source

Spacetrack Report No. 3, Hoots & Roehrich (1980) — Pages 67–80 of the source listing. 451 lines, the largest single compilation unit in the codebase.

The DEEP subroutine is the most complex component of the STR#3 codebase. It implements lunar-solar gravitational perturbations for deep-space orbits (period $\geq$ 225 minutes) and handles both synchronous (24-hour, GEO) and 12-hour (Molniya) resonance terms.

Architecture: Three Entry Points

DEEP is actually three routines sharing a single set of local variables via the FORTRAN IV ENTRY statement:

Entry Point	Called	Purpose
DPINIT(...)	Once per TLE	Initialize lunar-solar and resonance coefficients
DPSEC(...)	Each timestep	Apply secular perturbations + numerical integrator
DPPER(...)	Each timestep	Apply periodic perturbations (Lyddane modification)

Cross-entry argument sharing

In FORTRAN IV, all entry points shared a single fixed argument area. DPPER accesses T (a DPSEC parameter) and SINIQ/COSIQ (both DPINIT parameters). Modern gfortran gives each ENTRY its own calling convention, so accessing another entry’s arguments dereferences a stale pointer. The fix saves cross-entry arguments into local variables: TSAVE, SIQSAV, CIQSAV, OGDSAV.

Source Code

*     DEEP SPACE                                        31 OCT 80
      SUBROUTINE DEEP
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
     1          XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
     1          XJ3,XKE,XKMPER,XMNPDA,AE
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      DOUBLE PRECISION EPOCH, DS50
      DOUBLE PRECISION
     *    DAY,PREEP,XNODCE,ATIME,DELT,SAVTSN,STEP2,STEPN,STEPP,T,
     *    TSAVE,SIQSAV,CIQSAV,OGDSAV
      DATA            ZNS,              C1SS,           ZES/
     A                1.19459E-5,       2.9864797E-6,   .01675/
      DATA            ZNL,              C1L,            ZEL/
     A                1.5835218E-4,     4.7968065E-7,   .05490/
      DATA            ZCOSIS,           ZSINIS,         ZSINGS/
     A                .91744867,        .39785416,      -.98088458/
      DATA            ZCOSGS,           ZCOSHS,         ZSINHS/
     A                .1945905,         1.0,            0.0/
      DATA Q22,Q31,Q33/1.7891679E-6,2.1460748E-6,2.2123015E-7/
      DATA G22,G32/5.7686396,0.95240898/
      DATA G44,G52/1.8014998,1.0508330/
      DATA G54/4.4108898/
      DATA ROOT22,ROOT32/1.7891679E-6,3.7393792E-7/
      DATA ROOT44,ROOT52/7.3636953E-9,1.1428639E-7/
      DATA ROOT54/2.1765803E-9/
      DATA THDT/4.3752691E-3/

DPINIT — Initialization

*
*     ENTRANCE FOR DEEP SPACE INITIALIZATION
*
      ENTRY DPINIT(EQSQ,SINIQ,COSIQ,RTEQSQ,AO,COSQ2,SINOMO,COSOMO,
     1             BSQ,XLLDOT,OMGDT,XNODOT,XNODP)
      SIQSAV=SINIQ
      CIQSAV=COSIQ
      OGDSAV=OMGDT
      THGR=THETAG(EPOCH)
      EQ = EO
      XNQ = XNODP
      AQNV = 1./AO
      XQNCL = XINCL
      XMAO=XMO
      XPIDOT=OMGDT+XNODOT
      SINQ = SIN(XNODEO)
      COSQ = COS(XNODEO)
      OMEGAQ = OMEGAO

The initialization computes lunar and solar perturbation coefficients using an ASSIGN/GO TO dispatch pattern — first solar terms (label 20→30), then lunar terms (label 20→40) reusing the same coefficient calculation code. This is essentially a “function pointer to a label,” standard FORTRAN 77 but removed in Fortran 95.

Resonance Detection

The orbit is classified as resonant by checking the mean motion:

*     GEOPOTENTIAL RESONANCE INITIALIZATION FOR 12 HOUR ORBITS
      IRESFL=0
      ISYNFL=0
      IF(XNQ.LT.(.0052359877).AND.XNQ.GT.(.0034906585)) GO TO 70
      IF (XNQ.LT.(8.26E-3) .OR. XNQ.GT.(9.24E-3))       RETURN
      IF (EQ.LT.0.5)        RETURN
      IRESFL =1

• Synchronous (24h): Mean motion near $2\pi/1440$ rad/min (GEO orbits)
• 12-hour: Mean motion near $4\pi/1440$, eccentricity $> 0.5$ (Molniya orbits)

DPSEC — Secular Integrator

*     ENTRANCE FOR DEEP SPACE SECULAR EFFECTS
      ENTRY DPSEC(XLL,OMGASM,XNODES,EM,XINC,XN,T)
      TSAVE=T
      XLL=XLL+SSL*T
      OMGASM=OMGASM+SSG*T
      XNODES=XNODES+SSH*T
      EM=EO+SSE*T
      XINC=XINCL+SSI*T

For resonant orbits, the integrator steps in 720-minute intervals using a midpoint method. The integrator always restarts from epoch (label 170), which avoids the direction-dependent bug documented in the Rev-1 paper.

DPPER — Periodic Perturbations with Lyddane Modification

      ENTRY DPPER(EM,XINC,OMGASM,XNODES,XLL)
      SINIS = SIN(XINC)
      COSIS = COS(XINC)
      IF (DABS(SAVTSN-TSAVE).LT.(30.D0))    GO TO 210

The Lyddane modification (label 220) avoids singularities near zero inclination by using auxiliary variables ALFDP, BETDP instead of the node angle directly. The threshold for switching between direct application and Lyddane modification is inclination $< 0.2$ radians ($\approx 11.5°$).

Notable Design Patterns

ASSIGN/GO TO dispatch — Used to execute the same coefficient calculation code for both solar and lunar terms, with different input constants. The ASSIGN 30 TO LS and ASSIGN 40 TO LS statements set the return target, and GO TO LS dispatches to it.

720-minute step integrator — The resonance integrator uses fixed 720-minute steps (STEPP=720.D0), saving the intermediate state (ATIME, XLI, XNI). The STEP2 = 259200.D0 value is $720^2/2$, the half-step factor for the midpoint integration formula.

Static state — All local variables persist between calls via -fno-automatic. The 24 DATA constants, the resonance coefficients, and the integrator state all survive across the DPINIT → DPSEC → DPPER call sequence within a single propagation run, and across multiple timestep calls.