The SGP Model

Source

Spacetrack Report No. 3, Hoots & Roehrich (1980), compiled by T.S. Kelso (1988) — Pages 6—12 (pages 3—9 in original numbering)

The NORAD mean element sets can be used for prediction with SGP. All symbols not defined below are defined in the list of symbols in Section Twelve. Predictions are made by first calculating the constants

$$a_1 = \left( \frac{k_e}{n_o} \right)^{\frac{2}{3}}$$ $$\delta_1 = \frac{3}{4} J_2 \frac{{a_E}^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]$$ $$p_o = a_o (1 - {e_o}^2)$$ $$q_o = a_o (1 - e_o)$$ $$L_o = M_o + \omega_o + \Omega_o$$ $$\frac{d\Omega}{dt} = -\frac{3}{2} J_2 \frac{{a_E}^2}{{p_o}^2} n_o \cos i_o$$ $$\frac{d\omega}{dt} = \frac{3}{4} J_2 \frac{{a_E}^2}{{p_o}^2} n_o (5 \cos^2 i_o - 1).$$

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

$$a = a_o \left\{ \frac{n_o}{n_o + 2 \left( \frac{\dot{n}_o}{2} \right)(t - t_o) + 3 \left( \frac{\ddot{n}_o}{6} \right)(t - t_o)^2} \right\}^{\frac{2}{3}}$$ $$e = \begin{cases} 1 - \dfrac{q_o}{a}, & \text{for } a > q_o \\[6pt] 10^{-6}, & \text{for } a \le q_o \end{cases}$$ $$p = a(1 - e^2)$$ $$\Omega_{s_o} = \Omega_o + \frac{d\Omega}{dt}(t - t_o)$$ $$\omega_{s_o} = \omega_o + \frac{d\omega}{dt}(t - t_o)$$ $$L_s = L_o + \left( n_o + \frac{d\omega}{dt} + \frac{d\Omega}{dt} \right)(t - t_o) + \frac{\dot{n}_o}{2}(t - t_o)^2 + \frac{\ddot{n}_o}{6}(t - t_o)^3$$

where $(t - t_o)$ is time since epoch.

Long-period periodics are included through the equations

$$a_{yNSL} = e \sin \omega_{s_o} - \frac{1}{2} \frac{J_3}{J_2} \frac{a_E}{p} \sin i_o$$ $$L = L_s - \frac{1}{4} \frac{J_3}{J_2} \frac{a_E}{p} a_{xNSL} \sin i_o \left[ \frac{3 + 5 \cos i_o}{1 + \cos i_o} \right]$$

where

$$a_{xNSL} = e \cos \omega_{s_o}.$$

Solve Kepler’s equation for $E + \omega$ (by iteration to the desired accuracy), where

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

with

$$\Delta(E + \omega)_i = \frac{U - a_{yNSL} \cos(E + \omega)_i + a_{xNSL} \sin(E + \omega)_i - (E + \omega)_i}{-a_{yNSL} \sin(E + \omega)_i - a_{xNSL} \cos(E + \omega)_i + 1}$$ $$U = L - \Omega_{s_o}$$

and

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

Then calculate the intermediate (partially osculating) quantities

$$e \cos E = a_{xNSL} \cos(E + \omega) + a_{yNSL} \sin(E + \omega)$$ $$e \sin E = a_{xNSL} \sin(E + \omega) - a_{yNSL} \cos(E + \omega)$$ $${e_L}^2 = (a_{xNSL})^2 + (a_{yNSL})^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{v} = k_e \frac{\sqrt{p_L}}{r}$$ $$\sin u = \frac{a}{r} \left[ \sin(E + \omega) - a_{yNSL} - a_{xNSL} \frac{e \sin E}{1 + \sqrt{1 - {e_L}^2}} \right]$$ $$\cos u = \frac{a}{r} \left[ \cos(E + \omega) - a_{xNSL} + a_{yNSL} \frac{e \sin E}{1 + \sqrt{1 - {e_L}^2}} \right]$$ $$u = \tan^{-1} \left( \frac{\sin u}{\cos u} \right).$$

Short-period perturbations are now included by

$$r_k = r + \frac{1}{4} J_2 \frac{{a_E}^2}{p_L} \sin^2 i_o \cos 2u$$ $$u_k = u - \frac{1}{8} J_2 \frac{{a_E}^2}{{p_L}^2} (7 \cos^2 i_o - 1) \sin 2u$$ $$\Omega_k = \Omega_{s_o} + \frac{3}{4} J_2 \frac{{a_E}^2}{{p_L}^2} \cos i_o \sin 2u$$ $$i_k = i_o + \frac{3}{4} J_2 \frac{{a_E}^2}{{p_L}^2} \sin i_o \cos i_o \cos 2u.$$

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} \mathbf{U} + (r\dot{v}) \mathbf{V}.$$

A FORTRAN IV computer code listing of the subroutine SGP is given below.

---

FORTRAN Implementation

*         SGP                                         31 OCT 80
      SUBROUTINE SGP(IFLAG,TSINCE)
      COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
     1          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 19
*         INITIALIZATION
      C1= CK2*1.5
      C2= CK2/4.0
      C3= CK2/2.0
      C4= XJ3*AE**3/(4.0*CK2)
      COSIO=COS(XINCL)
      SINIO=SIN(XINCL)
      A1=(XKE/XNO)**TOTHRD
      D1=     C1/A1/A1*(3.*COSIO*COSIO-1.)/(1.-EO*EO)**1.5
      AO=A1*(1.-1./3.*D1-D1*D1-134./81.*D1*D1*D1)
      PO=AO*(1.-EO*EO)
      QO=AO*(1.-EO)
      XLO=XMO+OMEGAO+XNODEO
      D1O= C3 *SINIO*SINIO
      D2O= C2 *(7.*COSIO*COSIO-1.)
      D3O=C1*COSIO
      D4O=D3O*SINIO
      PO2NO=XNO/(PO*PO)
      OMGDT=C1*PO2NO*(5.*COSIO*COSIO-1.)
      XNODOT=-2.*D3O*PO2NO
      C5=.5*C4*SINIO*(3.+5.*COSIO)/(1.+COSIO)
      C6=C4*SINIO
      IFLAG=0
*         UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
   19 A=XNO+(2.*XNDT2O+3.*XNDD6O*TSINCE)*TSINCE
      A=AO*(XNO/A)**TOTHRD
      E=E6A
      IF(A.GT.QO) E=1.-QO/A
      P=A*(1.-E*E)
      XNODES= XNODEO+XNODOT*TSINCE
      OMGAS= OMEGAO+OMGDT*TSINCE
      XLS=FMOD2P(XLO+(XNO+OMGDT+XNODOT+(XNDT2O+XNDD6O*TSINCE)*
     1 TSINCE)*TSINCE)
*         LONG PERIOD PERIODICS
      AXNSL=E*COS(OMGAS)
      AYNSL=E*SIN(OMGAS)-C6/P
      XL=FMOD2P(XLS-C5/P*AXNSL)
*         SOLVE KEPLERS EQUATION
      U=FMOD2P(XL-XNODES)
      ITEM3=0
      EO1=U
      TEM5=1.
   20 SINEO1=SIN(EO1)
      COSEO1=COS(EO1)
      IF(ABS(TEM5).LT.E6A) GO TO 30
      IF(ITEM3.GE.10) GO TO 30
      ITEM3=ITEM3+1
      TEM5=1.-COSEO1*AXNSL-SINEO1*AYNSL
      TEM5=(U-AYNSL*COSEO1+AXNSL*SINEO1-EO1)/TEM5
      TEM2=ABS(TEM5)
      IF(TEM2.GT.1.) TEM5=TEM2/TEM5
      EO1=EO1+TEM5
      GO TO 20
*         SHORT PERIOD PRELIMINARY QUANTITIES
   30 ECOSE=AXNSL*COSEO1+AYNSL*SINEO1
      ESINE=AXNSL*SINEO1-AYNSL*COSEO1
      EL2=AXNSL*AXNSL+AYNSL*AYNSL
      PL=A*(1.-EL2)
      PL2=PL*PL
      R=A*(1.-ECOSE)
      RDOT=XKE*SQRT(A)/R*ESINE
      RVDOT=XKE*SQRT(PL)/R
      TEMP=ESINE/(1.+SQRT(1.-EL2))
      SINU=A/R*(SINEO1-AYNSL-AXNSL*TEMP)
      COSU=A/R*(COSEO1-AXNSL+AYNSL*TEMP)
      SU=ACTAN(SINU,COSU)
*         UPDATE FOR SHORT PERIODICS
      SIN2U=(COSU+COSU)*SINU
      COS2U=1.-2.*SINU*SINU
      RK=R+D1O/PL*COS2U
      UK=SU-D2O/PL2*SIN2U
      XNODEK=XNODES+D3O*SIN2U/PL2
      XINCK =XINCL+D4O/PL2*COS2U
*         ORIENTATION VECTORS
      SINUK=SIN(UK)
      COSUK=COS(UK)
      SINNOK=SIN(XNODEK)
      COSNOK=COS(XNODEK)
      SINIK=SIN(XINCK)
      COSIK=COS(XINCK)
      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=RDOT*UX
      YDOT=RDOT*UY
      ZDOT=RDOT*UZ
      XDOT=RVDOT*VX+XDOT
      YDOT=RVDOT*VY+YDOT
      ZDOT=RVDOT*VZ+ZDOT
      RETURN
      END