The SGP4 Model

The NORAD mean element sets can be used for prediction with SGP4. All symbols not defined below are defined in the list of symbols in Section Twelve. The original mean motion ( $n''_o$ ) and semimajor axis ( $a''_o$ ) are first recovered from the input elements by the equations

$a_1 = \left(\frac{k_e}{n_o}\right)^{\frac{2}{3}}$

$\delta_1 = \frac{3}{2} \frac{k_2}{a_1^{\ 2}} \frac{(3\cos^2 i_o - 1)}{(1 - e_o^{\ 2})^{\frac{3}{2}}}$

$a_o = a_1 \left(1 - \frac{1}{3}\delta_1 - \delta_1^{\ 2} - \frac{134}{81}\delta_1^{\ 3}\right)$

$\delta_o = \frac{3}{2} \frac{k_2}{a_o^{\ 2}} \frac{(3\cos^2 i_o - 1)}{(1 - e_o^{\ 2})^{\frac{3}{2}}}$

$n''_o = \frac{n_o}{1 + \delta_o}$

$a''_o = \frac{a_o}{1 - \delta_o}.$

For perigee between 98 kilometers and 156 kilometers, the value of the constant $s$ used in SGP4 is changed to

$s^* = a''_o(1 - e_o) - s + a_E$

For perigee below 98 kilometers, the value of $s$ is changed to

$s^* = 20/\text{XKMPER} + a_E.$

If the value of $s$ is changed, then the value of $(q_o - s)^4$ must be replaced by

$(q_o - s^*)^4 = \left[[(q_o - s)^4]^{\frac{1}{4}} + s - s^*\right]^4.$

Then calculate the constants (using the appropriate values of $s$ and $(q_o - s)^4$ )

$\theta = \cos i_o$

$\xi = \frac{1}{a''_o - s}$

$\beta_o = (1 - e_o^{\ 2})^{\frac{1}{2}}$

$\eta = a''_o e_o \xi$

$C_2 = (q_o - s)^4 \xi^4 n''_o (1 - \eta^2)^{-\frac{7}{2}} \left[a''_o \left(1 + \frac{3}{2}\eta^2 + 4e_o\eta + e_o\eta^3\right) + \frac{3}{2} \frac{k_2 \xi}{(1 - \eta^2)} \left(-\frac{1}{2} + \frac{3}{2}\theta^2\right)(8 + 24\eta^2 + 3\eta^4)\right]$

$C_1 = B^* C_2$

$C_3 = \frac{(q_o - s)^4 \xi^5 A_{3,0} n''_o a_E \sin i_o}{k_2 e_o}$

$C_4 = 2n''_o (q_o - s)^4 \xi^4 a''_o \beta_o^{\ 2} (1 - \eta^2)^{-\frac{7}{2}} \left(\left[2\eta(1 + e_o\eta) + \frac{1}{2}e_o + \frac{1}{2}\eta^3\right] - \frac{2k_2\xi}{a''_o(1 - \eta^2)} \times \left[3(1 - 3\theta^2)\left(1 + \frac{3}{2}\eta^2 - 2e_o\eta - \frac{1}{2}e_o\eta^3\right) + \frac{3}{4}(1 - \theta^2)(2\eta^2 - e_o\eta - e_o\eta^3)\cos 2\omega_o\right]\right)$

$C_5 = 2(q_o - s)^4 \xi^4 a''_o \beta_o^{\ 2} (1 - \eta^2)^{-\frac{7}{2}} \left[1 + \frac{11}{4}\eta(\eta + e_o) + e_o\eta^3\right]$

$D_2 = 4a''_o \xi C_1^{\ 2}$

$D_3 = \frac{4}{3}a''_o \xi^2 (17a''_o + s) C_1^{\ 3}$

$D_4 = \frac{2}{3}a''_o \xi^3 (221a''_o + 31s) C_1^{\ 4}.$

The secular effects of atmospheric drag and gravitation are included through the equations

$M_{DF} = M_o + \left[1 + \frac{3k_2(-1 + 3\theta^2)}{2a''^{\ 2}_o \beta_o^{\ 3}} + \frac{3k_2^{\ 2}(13 - 78\theta^2 + 137\theta^4)}{16a''^{\ 4}_o \beta_o^{\ 7}}\right] n''_o(t - t_o)$

$\omega_{DF} = \omega_o + \left[\frac{-3k_2(1 - 5\theta^2)}{2a''^{\ 2}_o \beta_o^{\ 4}} + \frac{3k_2^{\ 2}(7 - 114\theta^2 + 395\theta^4)}{16a''^{\ 4}_o \beta_o^{\ 8}} + \frac{5k_4(3 - 36\theta^2 + 49\theta^4)}{4a''^{\ 4}_o \beta_o^{\ 8}}\right] n''_o(t - t_o)$

$\Omega_{DF} = \Omega_o + \left[\frac{-3k_2\theta}{a''^{\ 2}_o \beta_o^{\ 4}} + \frac{3k_2^{\ 2}(4\theta - 19\theta^3)}{2a''^{\ 4}_o \beta_o^{\ 8}} + \frac{5k_4\theta(3 - 7\theta^2)}{2a''^{\ 4}_o \beta_o^{\ 8}}\right] n''_o(t - t_o)$

$\delta\omega = B^* C_3 (\cos \omega_o)(t - t_o)$

$\delta M = -\frac{2}{3}(q_o - s)^4 B^* \xi^4 \frac{a_E}{e_o \eta}\left[(1 + \eta \cos M_{DF})^3 - (1 + \eta \cos M_o)^3\right]$

$M_p = M_{DF} + \delta\omega + \delta M$

$\omega = \omega_{DF} - \delta\omega - \delta M$

$\Omega = \Omega_{DF} - \frac{21}{2} \frac{n''_o k_2 \theta}{a''^{\ 2}_o \beta_o^{\ 2}} C_1 (t - t_o)^2$

$e = e_o - B^* C_4 (t - t_o) - B^* C_5 (\sin M_p - \sin M_o)$

$a = a''_o \left[1 - C_1(t - t_o) - D_2(t - t_o)^2 - D_3(t - t_o)^3 - D_4(t - t_o)^4\right]^2$

$\mathbb{L} = M_p + \omega + \Omega + n''_o \left[\frac{3}{2}C_1(t - t_o)^2 + (D_2 + 2C_1^{\ 2})(t - t_o)^3 + \frac{1}{4}(3D_3 + 12C_1 D_2 + 10C_1^{\ 3})(t - t_o)^4 + \frac{1}{5}(3D_4 + 12C_1 D_3 + 6D_2^{\ 2} + 30C_1^{\ 2} D_2 + 15C_1^{\ 4})(t - t_o)^5\right]$

$\beta = \sqrt{(1 - e^2)}$

$n = k_e \Big/ a^{\frac{3}{2}}$

where $(t - t_o)$ is time since epoch. It should be noted that when epoch perigee height is less than 220 kilometers, the equations for $a$ and $\mathbb{L}$ are truncated after the $C_1$ term, and the terms involving $C_5$ , $\delta\omega$ , and $\delta M$ are dropped.

Add the long-period periodic terms

$a_{xN} = e \cos \omega$

$\mathbb{L}_L = \frac{A_{3,0} \sin i_o}{8 k_2 a \beta^2} (e \cos \omega) \left(\frac{3 + 5\theta}{1 + \theta}\right)$

$a_{yNL} = \frac{A_{3,0} \sin i_o}{4 k_2 a \beta^2}$

$\mathbb{L}_T = \mathbb{L} + \mathbb{L}_L$

$a_{yN} = e \sin \omega + a_{yNL}.$

Solve Kepler’s equation for $(E + \omega)$ by defining

$U = \mathbb{L}_T - \Omega$

and using the iteration equation

$(E + \omega)_{i+1} = (E + \omega)_i + \Delta(E + \omega)_i$

with

$\Delta(E + \omega)_i = \frac{U - a_{yN}\cos(E + \omega)_i + a_{xN}\sin(E + \omega)_i - (E + \omega)_i}{-a_{yN}\sin(E + \omega)_i - a_{xN}\cos(E + \omega)_i + 1}$

and

$(E + \omega)_1 = U.$

The following equations are used to calculate preliminary quantities needed for short-period periodics.

$e \cos E = a_{xN} \cos(E + \omega) + a_{yN} \sin(E + \omega)$

$e \sin E = a_{xN} \sin(E + \omega) - a_{yN} \cos(E + \omega)$

$e_L = (a_{xN}^{\ 2} + a_{yN}^{\ 2})^{\frac{1}{2}}$

$p_L = a(1 - e_L^{\ 2})$

$r = a(1 - e \cos E)$

$\dot{r} = k_e \frac{\sqrt{a}}{r} e \sin E$

$r\dot{f} = k_e \frac{\sqrt{p_L}}{r}$

$\cos u = \frac{a}{r} \left[\cos(E + \omega) - a_{xN} + \frac{a_{yN}(e \sin E)}{1 + \sqrt{1 - e_L^{\ 2}}}\right]$

$\sin u = \frac{a}{r} \left[\sin(E + \omega) - a_{yN} - \frac{a_{xN}(e \sin E)}{1 + \sqrt{1 - e_L^{\ 2}}}\right]$

$u = \tan^{-1}\left(\frac{\sin u}{\cos u}\right)$

$\Delta r = \frac{k_2}{2 p_L}(1 - \theta^2) \cos 2u$

$\Delta u = -\frac{k_2}{4 p_L^{\ 2}}(7\theta^2 - 1) \sin 2u$

$\Delta\Omega = \frac{3 k_2 \theta}{2 p_L^{\ 2}} \sin 2u$

$\Delta i = \frac{3 k_2 \theta}{2 p_L^{\ 2}} \sin i_o \cos 2u$

$\Delta\dot{r} = -\frac{k_2 n}{p_L}(1 - \theta^2) \sin 2u$

$\Delta r\dot{f} = \frac{k_2 n}{p_L}\left[(1 - \theta^2) \cos 2u - \frac{3}{2}(1 - 3\theta^2)\right]$

The short-period periodics are added to give the osculating quantities

$r_k = r \left[1 - \frac{3}{2} k_2 \frac{\sqrt{1 - e_L^{\ 2}}}{p_L^{\ 2}}(3\theta^2 - 1)\right] + \Delta r$

$u_k = u + \Delta u$

$\Omega_k = \Omega + \Delta\Omega$

$i_k = i_o + \Delta i$

$\dot{r}_k = \dot{r} + \Delta\dot{r}$

$r\dot{f}_k = r\dot{f} + \Delta r\dot{f}.$

Then unit orientation vectors are calculated by

$\mathbf{U} = \mathbf{M} \sin u_k + \mathbf{N} \cos u_k$

$\mathbf{V} = \mathbf{M} \cos u_k - \mathbf{N} \sin u_k$

where

$\mathbf{M} = \begin{cases} M_x = -\sin \Omega_k \cos i_k \\ M_y = \cos \Omega_k \cos i_k \\ M_z = \sin i_k \end{cases}$

$\mathbf{N} = \begin{cases} N_x = \cos \Omega_k \\ N_y = \sin \Omega_k \\ N_z = 0 \end{cases}.$

Then position and velocity are given by

$\mathbf{r} = r_k \mathbf{U}$

and

$\dot{\mathbf{r}} = \dot{r}_k \mathbf{U} + (r\dot{f})_k \mathbf{V}.$

A FORTRAN IV computer code listing of the subroutine SGP4 is given below. These equations contain all currently anticipated changes to the SCC operational program. These changes are scheduled for implementation in March, 1981.

FORTRAN Implementation

*     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