The SDP4 Model

Source

Spacetrack Report No. 3, Hoots & Roehrich (1980), compiled by T.S. Kelso (1988) — Pages 24—33 (pages 21—30 in original numbering)

Recovery of Original Mean Motion and Semimajor Axis

(Page 24 / original p. 21)

The NORAD mean element sets can be used for prediction with SDP4. 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}.$

Perigee-Dependent Constants

For perigee between 98 kilometers and 156 kilometers, the value of the constant $s$ used in SDP4 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.$

Initialization Constants

(Page 25 / original p. 22)

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_4 = 2 n_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)$

$\dot{M} = \left[ 1 + \frac{3k_2(-1 + 3\theta^2)}{2{a_o''}^2 {\beta_o}^3} + \frac{3{k_2}^2(13 - 78\theta^2 + 137\theta^4)}{16{a_o''}^4 {\beta_o}^7} \right] n_o''$

$\dot{\omega} = \left[ \frac{-3k_2(1 - 5\theta^2)}{2{a_o''}^2 {\beta_o}^4} + \frac{3{k_2}^2(7 - 114\theta^2 + 395\theta^4)}{16{a_o''}^4 {\beta_o}^8} + \frac{5k_4(3 - 36\theta^2 + 49\theta^4)}{4{a_o''}^4 {\beta_o}^8} \right] n_o''$

$\dot{\Omega}_1 = -\frac{3k_2\theta}{{a_o''}^2 {\beta_o}^4} n_o''$

$\dot{\Omega} = \dot{\Omega}_1 + \left[ \frac{3{k_2}^2(4\theta - 19\theta^3)}{2{a_o''}^4 {\beta_o}^8} + \frac{5k_4\theta(3 - 7\theta^2)}{2{a_o''}^4 {\beta_o}^8} \right] n_o''.$

Deep-Space Initialization

At this point SDP4 calls the initialization section of DEEP which calculates all initialized quantities needed for the deep-space perturbations (see Section Ten).

Secular Effects of Gravity

The secular effects of gravity are included by

$M_{DF} = M_o + \dot{M}(t - t_o)$

$\omega_{DF} = \omega_o + \dot{\omega}(t - t_o)$

(Page 26 / original p. 23)

$\Omega_{DF} = \Omega_o + \dot{\Omega}(t - t_o)$

where $(t - t_o)$ is time since epoch. The secular effect of drag on longitude of ascending node is included by

$\Omega = \Omega_{DF} - \frac{21}{2} \frac{n_o'' k_2 \theta}{{a_o''}^2 {\beta_o}^2} C_1 (t - t_o)^2.$

Deep-Space Secular and Resonance Effects

Next, SDP4 calls the secular section of DEEP which adds the deep-space secular effects and long-period resonance effects to the six classical orbital elements (see Section Ten).

Secular Effects of Drag

The secular effects of drag are included in the remaining elements by

$a = a_{DS} [1 - C_1(t - t_o)]^2$

$e = e_{DS} - B^* C_4 (t - t_o)$

$\mathcal{L} = M_{DS} + \omega_{DS} + \Omega_{DS} + n_o'' \left[ \frac{3}{2} C_1 (t - t_o)^2 \right]$

where $a_{DS}$, $e_{DS}$, $M_{DS}$, $\omega_{DS}$, and $\Omega_{DS}$, are the values of $n_o$, $e_o$, $M_{DF}$, $\omega_{DF}$, and $\Omega$ after deep-space secular and resonance perturbations have been applied.

Deep-Space Lunar-Solar Periodics

Here SDP4 calls the periodics section of DEEP which adds the deep-space lunar and solar periodics to the orbital elements (see Section Ten). From this point on, it will be assumed that $n$, $e$, $I$, $\omega$, $\Omega$, and $M$ are the mean motion, eccentricity, inclination, argument of perigee, longitude of ascending node, and mean anomaly after lunar-solar periodics have been added.

Long-Period Periodic Terms

Add the long-period periodic terms

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

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

$\mathcal{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}$

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

(Page 27 / original p. 24)

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

Kepler’s Equation

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

$U = \mathcal{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.$

Short-Period Preliminary Quantities

(Page 27—28 / original p. 24—25)

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]$

(Page 28 / original p. 25)

$\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)$

Short-Period Perturbations

$\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]$

Osculating Quantities

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 + \Delta i$

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

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

Orientation Vectors and Position/Velocity

(Page 29 / original p. 26)

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 SDP4 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

(Pages 30—33 / original p. 27—30)

*         SDP4                                         3 NOV 80
      SUBROUTINE SDP4(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 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)
      SING=SIN(OMEGAO)
      COSG=COS(OMEGAO)
      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
      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)))
      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
      XNODCF=3.5*BETAO2*XHDOT1*C1
      T2COF=1.5*C1
      XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO)
      AYCOF=.25*A3OVK2*SINIO
      X7THM1=7.*THETA2-1.
   90 IFLAG=0
      CALL DPINIT(EOSQ,SINIO,COSIO,BETAO,AODP,THETA2,
     1   SING,COSG,BETAO2,XMDOT,OMGDOT,XNODOT,XNODP)

*     UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG

  100 XMDF=XMO+XMDOT*TSINCE
      OMGADF=OMEGAO+OMGDOT*TSINCE
      XNODDF=XNODEO+XNODOT*TSINCE
      TSQ=TSINCE*TSINCE
      XNODE=XNODDF+XNODCF*TSQ
      TEMPA=1.-C1*TSINCE
      TEMPE=BSTAR*C4*TSINCE
      TEMPL=T2COF*TSQ
      XN=XNODP
      CALL DPSEC(XMDF,OMGADF,XNODE,EM,XINC,XN,TSINCE)
      A=(XKE/XN)**TOTHRD*TEMPA**2
      E=EM-TEMPE
      XMAM=XMDF+XNODP*TEMPL
      CALL DPPER(E,XINC,OMGADF,XNODE,XMAM)
      XL=XMAM+OMGADF+XNODE
      BETA=SQRT(1.-E*E)
      XN=XKE/A**1.5

*     LONG PERIOD PERIODICS

      AXN=E*COS(OMGADF)
      TEMP=1./(A*BETA*BETA)
      XLL=TEMP*XLCOF*AXN
      AYNL=TEMP*AYCOF
      XLT=XL+XLL
      AYN=E*SIN(OMGADF)+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=XINC+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

---

Notes on SDP4 vs SGP4

SDP4 extends SGP4 for deep-space satellites (orbital period $\geq$ 225 minutes). The key differences are:

1. Shared initialization — The recovery of original mean motion ($n_o''$, $a_o''$), perigee-dependent constant adjustments, and the computation of $C_1$, $C_2$, $C_4$, $\dot{M}$, $\dot{\omega}$, $\dot{\Omega}$ are identical to SGP4.

2. Deep-space calls — SDP4 inserts three calls to the DEEP subroutine (Section Ten) that SGP4 does not make:

   • DPINIT — initialization of deep-space quantities (called once)
   • DPSEC — secular effects from lunar/solar gravity and resonance perturbations (called each timestep)
   • DPPER — periodic effects from lunar/solar gravity (called each timestep)

3. Simplified drag model — SDP4 uses a simpler atmospheric drag model than SGP4. There are no $C_5$, $D_2$, $D_3$, $D_4$ terms and no cubic/quartic/quintic secular drag corrections. The drag terms are truncated at second order in $C_1$ because deep-space satellites experience negligible atmospheric drag.

4. Identical short-period corrections — The long-period periodics, Kepler’s equation solution, short-period preliminary quantities, short-period perturbations, osculating elements, orientation vectors, and final position/velocity computation are identical between SGP4 and SDP4.