SGP4 & SGP8 Models

SGP4 — Brouwer Mean Elements (215 lines)

SGP4 is the workhorse model, used for all near-earth objects (period < 225 minutes) in the NORAD catalog. It implements an analytic theory based on Brouwer’s solution to the Vinti problem, with Lane’s atmospheric drag model.

The code structure follows an initialization/propagation pattern controlled by IFLAG:

IFLAG=1: Compute all constants from the mean elements (runs once per TLE)
IFLAG=0: Propagate to TSINCE minutes from epoch using saved constants

The perigee height determines three drag regimes:

Above 220 km: Full equations with higher-order secular terms ( $D_2$ , $D_3$ , $D_4$ )
98–156 km: Modified atmospheric parameter $s^*$ replaces the standard $s$
Below 98 km: Simplified $s = 20/\text{XKMPER} + a_E$

*     SGP4                                           3 NOV 80
      SUBROUTINE SGP4(IFLAG,TSINCE)
      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
      DOUBLE PRECISION EPOCH, DS50
      IF (IFLAG .EQ. 0) GO TO 100
*
*     RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
*     FROM INPUT ELEMENTS
*
      A1=(XKE/XNO)**TOTHRD
      COSIO=COS(XINCL)
      THETA2=COSIO*COSIO
      X3THM1=3.*THETA2-1.
      EOSQ=EO*EO
      BETAO2=1.-EOSQ
      BETAO=SQRT(BETAO2)
      DEL1=1.5*CK2*X3THM1/(A1*A1*BETAO*BETAO2)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=1.5*CK2*X3THM1/(AO*AO*BETAO*BETAO2)
      XNODP=XNO/(1.+DELO)
      AODP=AO/(1.-DELO)
*
*     INITIALIZATION
*
*     FOR PERIGEE LESS THAN 220 KILOMETERS, THE ISIMP FLAG IS SET AND
*     THE EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN SQRT A AND
*     QUADRATIC VARIATION IN MEAN ANOMALY.  ALSO, THE C3 TERM, THE
*     DELTA OMEGA TERM, AND THE DELTA M TERM ARE DROPPED.
*
      ISIMP=0
      IF((AODP*(1.-EO)/AE) .LT. (220./XKMPER+AE)) ISIMP=1
*
*     FOR PERIGEE BELOW 156 KM, THE VALUES OF
*     S AND QOMS2T ARE ALTERED
*
      S4=S
      QOMS24=QOMS2T
      PERIGE=(AODP*(1.-EO)-AE)*XKMPER
      IF(PERIGE .GE. 156.) GO TO 10
      S4=PERIGE-78.
      IF(PERIGE .GT. 98.) GO TO 9
      S4=20.
    9 QOMS24=((120.-S4)*AE/XKMPER)**4
      S4=S4/XKMPER+AE
   10 PINVSQ=1./(AODP*AODP*BETAO2*BETAO2)
      TSI=1./(AODP-S4)
      ETA=AODP*EO*TSI
      ETASQ=ETA*ETA
      EETA=EO*ETA
      PSISQ=ABS(1.-ETASQ)
      COEF=QOMS24*TSI**4
      COEF1=COEF/PSISQ**3.5
      C2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*
     1   CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)))
      C1=BSTAR*C2
      SINIO=SIN(XINCL)
      A3OVK2=-XJ3/CK2*AE**3
      C3=COEF*TSI*A3OVK2*XNODP*AE*SINIO/EO
      X1MTH2=1.-THETA2
      C4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*
     1   (2.+.5*ETASQ)+EO*(.5+2.*ETASQ)-2.*CK2*TSI/
     2   (AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*
     3   (1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*
     4   (1.+ETASQ))*COS(2.*OMEGAO)))
      C5=2.*COEF1*AODP*BETAO2*(1.+2.75*(ETASQ+EETA)+EETA*ETASQ)
      THETA4=THETA2*THETA2
      TEMP1=3.*CK2*PINVSQ*XNODP
      TEMP2=TEMP1*CK2*PINVSQ
      TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP
      XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*
     1   (13.-78.*THETA2+137.*THETA4)
      X1M5TH=1.-5.*THETA2
      OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+
     1   395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4)
      XHDOT1=-TEMP1*COSIO
      XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-
     1   7.*THETA2))*COSIO
      OMGCOF=BSTAR*C3*COS(OMEGAO)
      XMCOF=-TOTHRD*COEF*BSTAR*AE/EETA
      XNODCF=3.5*BETAO2*XHDOT1*C1
      T2COF=1.5*C1
      XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO)
      AYCOF=.25*A3OVK2*SINIO
      DELMO=(1.+ETA*COS(XMO))**3
      SINMO=SIN(XMO)
      X7THM1=7.*THETA2-1.
      IF(ISIMP .EQ. 1) GO TO 90
      C1SQ=C1*C1
      D2=4.*AODP*TSI*C1SQ
      TEMP=D2*TSI*C1/3.
      D3=(17.*AODP+S4)*TEMP
      D4=.5*TEMP*AODP*TSI*(221.*AODP+31.*S4)*C1
      T3COF=D2+2.*C1SQ
      T4COF=.25*(3.*D3+C1*(12.*D2+10.*C1SQ))
      T5COF=.2*(3.*D4+12.*C1*D3+6.*D2*D2+15.*C1SQ*(
     1           2.*D2+C1SQ))
   90 IFLAG=0
*
*     UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
*
  100 XMDF=XMO+XMDOT*TSINCE
      OMGADF=OMEGAO+OMGDOT*TSINCE
      XNODDF=XNODEO+XNODOT*TSINCE
      OMEGA=OMGADF
      XMP=XMDF
      TSQ=TSINCE*TSINCE
      XNODE=XNODDF+XNODCF*TSQ
      TEMPA=1.-C1*TSINCE
      TEMPE=BSTAR*C4*TSINCE
      TEMPL=T2COF*TSQ
      IF(ISIMP .EQ. 1) GO TO 110
      DELOMG=OMGCOF*TSINCE
      DELM=XMCOF*((1.+ETA*COS(XMDF))**3-DELMO)
      TEMP=DELOMG+DELM
      XMP=XMDF+TEMP
      OMEGA=OMGADF-TEMP
      TCUBE=TSQ*TSINCE
      TFOUR=TSINCE*TCUBE
      TEMPA=TEMPA-D2*TSQ-D3*TCUBE-D4*TFOUR
      TEMPE=TEMPE+BSTAR*C5*(SIN(XMP)-SINMO)
      TEMPL=TEMPL+T3COF*TCUBE+
     1           TFOUR*(T4COF+TSINCE*T5COF)
  110 A=AODP*TEMPA**2
      E=EO-TEMPE
      XL=XMP+OMEGA+XNODE+XNODP*TEMPL
      BETA=SQRT(1.-E*E)
      XN=XKE/A**1.5
*
*     LONG PERIOD PERIODICS
*
      AXN=E*COS(OMEGA)
      TEMP=1./(A*BETA*BETA)
      XLL=TEMP*XLCOF*AXN
      AYNL=TEMP*AYCOF
      XLT=XL+XLL
      AYN=E*SIN(OMEGA)+AYNL
*
*     SOLVE KEPLERS EQUATION
*
      CAPU=FMOD2P(XLT-XNODE)
      TEMP2=CAPU
      DO 130 I=1,10
      SINEPW=SIN(TEMP2)
      COSEPW=COS(TEMP2)
      TEMP3=AXN*SINEPW
      TEMP4=AYN*COSEPW
      TEMP5=AXN*COSEPW
      TEMP6=AYN*SINEPW
      EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2
      IF(ABS(EPW-TEMP2) .LE. E6A) GO TO 140
  130 TEMP2=EPW
*
*     SHORT PERIOD PRELIMINARY QUANTITIES
*
  140 ECOSE=TEMP5+TEMP6
      ESINE=TEMP3-TEMP4
      ELSQ=AXN*AXN+AYN*AYN
      TEMP=1.-ELSQ
      PL=A*TEMP
      R=A*(1.-ECOSE)
      TEMP1=1./R
      RDOT=XKE*SQRT(A)*ESINE*TEMP1
      RFDOT=XKE*SQRT(PL)*TEMP1
      TEMP2=A*TEMP1
      BETAL=SQRT(TEMP)
      TEMP3=1./(1.+BETAL)
      COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3)
      SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3)
      U=ACTAN(SINU,COSU)
      SIN2U=2.*SINU*COSU
      COS2U=2.*COSU*COSU-1.
      TEMP=1./PL
      TEMP1=CK2*TEMP
      TEMP2=TEMP1*TEMP
*
*     UPDATE FOR SHORT PERIODICS
*
      RK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U
      UK=U-.25*TEMP2*X7THM1*SIN2U
      XNODEK=XNODE+1.5*TEMP2*COSIO*SIN2U
      XINCK=XINCL+1.5*TEMP2*COSIO*SINIO*COS2U
      RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U
      RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5*X3THM1)
*
*     ORIENTATION VECTORS
*
      SINUK=SIN(UK)
      COSUK=COS(UK)
      SINIK=SIN(XINCK)
      COSIK=COS(XINCK)
      SINNOK=SIN(XNODEK)
      COSNOK=COS(XNODEK)
      XMX=-SINNOK*COSIK
      XMY=COSNOK*COSIK
      UX=XMX*SINUK+COSNOK*COSUK
      UY=XMY*SINUK+SINNOK*COSUK
      UZ=SINIK*SINUK
      VX=XMX*COSUK-COSNOK*SINUK
      VY=XMY*COSUK-SINNOK*SINUK
      VZ=SINIK*COSUK
*
*     POSITION AND VELOCITY
*
      X=RK*UX
      Y=RK*UY
      Z=RK*UZ
      XDOT=RDOTK*UX+RFDOTK*VX
      YDOT=RDOTK*UY+RFDOTK*VY
      ZDOT=RDOTK*UZ+RFDOTK*VZ
      RETURN
      END

Key Implementation Notes

COMMON blocks — All models share two COMMON blocks:

/E1/ — Orbital elements and output position/velocity (input/output interface)
/C1/ — Physical and geopotential constants (WGS-72)

Static locals — The IFLAG pattern relies on FORTRAN IV’s static allocation: variables from the initialization pass (lines 8–100) survive to subsequent propagation calls (line 100+). This requires -fno-automatic on modern compilers.

Kepler’s equation — Solved via Newton-Raphson iteration (lines 143–154) with tolerance E6A = 1.0E-6. Limited to 10 iterations, which can diverge for $e > 0.9$ . The Rev-1 paper addresses this with a clamped first correction.

Perigee thresholds — Lines 33–47 implement three atmospheric density regimes. The threshold constants (220 km, 156 km, 98 km) are specific to the Lane drag model and must not be changed independently of the atmospheric density coefficients.