Utility Functions

Source

Spacetrack Report No. 3, Hoots & Roehrich (1980) — Pages 85–86 of the source listing

Three small utility functions support the propagation models. Despite their brevity, each solves a specific numerical problem that the standard FORTRAN IV library didn’t address.

ACTAN — Four-Quadrant Arctangent (22 lines)

The standard FORTRAN IV ATAN function returns values in $[-\pi/2, \pi/2]$. Orbital mechanics needs the full $[0, 2\pi)$ range. ACTAN takes separate sine and cosine arguments to determine the correct quadrant — the same concept as C’s atan2, which wasn’t available in FORTRAN IV.

      FUNCTION  ACTAN(SINX,COSX)
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      ACTAN=0.
      IF (COSX.EQ.0. ) GO TO 5
      IF (COSX.GT.0. ) GO TO 1
      ACTAN=PI
      GO TO 7
    1 IF (SINX.EQ.0. ) GO TO 8
      IF (SINX.GT.0. ) GO TO 7
      ACTAN=TWOPI
      GO TO 7
    5 IF (SINX.EQ.0. ) GO TO 8
      IF (SINX.GT.0. ) GO TO 6
      ACTAN=X3PIO2
      GO TO 8
    6 ACTAN=PIO2
      GO TO 8
    7 TEMP=SINX/COSX
      ACTAN=ACTAN+ATAN(TEMP)
    8 RETURN
      END

The logic manually determines the base angle ($0$, $\pi/2$, $\pi$, $3\pi/2$, or $2\pi$) based on the signs of SINX and COSX, then adds the standard ATAN correction. This guarantees a return value in $[0, 2\pi)$.

FMOD2P — Modulo $2\pi$ (11 lines)

Reduces any angle to the range $[0, 2\pi)$. Used after every angle accumulation to prevent floating-point overflow during long propagation runs.

      FUNCTION FMOD2P(X)
      COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
      FMOD2P=X
      I=FMOD2P/TWOPI
      FMOD2P=FMOD2P-I*TWOPI
      IF(FMOD2P.LT.0) FMOD2P=FMOD2P+TWOPI
      RETURN
      END

The integer truncation I=FMOD2P/TWOPI exploits FORTRAN’s implicit REAL→INTEGER conversion (truncation toward zero). The final IF handles negative inputs.

THETAG — Greenwich Sidereal Angle (22 lines)

Computes the Greenwich Sidereal Time (GST) at a given TLE epoch, expressed as an angle in radians. This is the rotation angle of the Earth relative to the vernal equinox — essential for converting the TEME output frame to Earth-fixed coordinates.

      FUNCTION THETAG(EP)
      COMMON /E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
     1 X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
      DOUBLE PRECISION EPOCH,D,THETA,TWOPI,YR,TEMP,EP,DS50
      TWOPI=6.28318530717959D0
      YR=(EP+2.D-7)*1.D-3
      JY=YR
      YR=JY
      D=EP-YR*1.D3
      IF(JY.LT.10) JY=JY+80
      N=(JY-69)/4
      IF(JY.LT.70) N=(JY-72)/4
      DS50=7305.D0 + 365.D0*(JY-70) +N + D
      THETA=1.72944494D0 + 6.3003880987D0*DS50
      TEMP=THETA/TWOPI
      I=TEMP
      TEMP=I
      THETAG=THETA-TEMP*TWOPI
      IF(THETAG.LT.0.D0) THETAG=THETAG+TWOPI
      RETURN
      END

The epoch format is YYDDD.DDDDDDDD (two-digit year + fractional day of year). The +2.D-7 offset on line 6 prevents truncation errors when the epoch falls exactly on a year boundary. The Y2K handling on line 10 (IF(JY.LT.10) JY=JY+80) maps years 00–09 to 2000–2009 and years 10–99 to 1910–1999.

The constant 1.72944494 is the Greenwich sidereal angle at the reference epoch (January 0.0, 1950), and 6.3003880987 is the Earth’s rotation rate in radians per day ($2\pi \times 1.00273790934$, accounting for the sidereal/solar day difference).

Side effect: THETAG also computes DS50 (days since 1950 January 0.0) as a side effect, storing it in the /E1/ COMMON block. The DEEP subroutine depends on this value being set before DPINIT is called.