The SDP8 Model

Source

Spacetrack Report No. 3, Hoots & Roehrich (1980), compiled by T.S. Kelso (1988) — Pages 50—56 (pages 47—53 in original numbering); FORTRAN listing pages 57—60 (54—57)

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

The ballistic coefficient ($B$ term) is then calculated from the $B^*$ drag term by

$$B = 2 B^* / \rho_o$$

where

$$\rho_o = (2.461 \times 10^{-5}) \text{ XKMPER kg/m}^2\text{/Earth radii}$$

is a reference value of atmospheric density.

Then calculate the constants

$$\beta^2 = 1 - e^2$$ $$\theta = \cos i$$ $$\dot{M}_1 = -\frac{3}{2} \frac{n'' k_2}{{a''}^2 \beta^3} (1 - 3\theta^2)$$ $$\dot{\omega}_1 = -\frac{3}{2} \frac{n'' k_2}{{a''}^2 \beta^4} (1 - 5\theta^2)$$ $$\dot{\Omega}_1 = -3 \frac{n'' k_2}{{a''}^2 \beta^4} \theta$$ $$\dot{M}_2 = \frac{3}{16} \frac{n'' {k_2}^2}{{a''}^4 \beta^7} (13 - 78\theta^2 + 137\theta^4)$$ $$\dot{\omega}_2 = \frac{3}{16} \frac{n'' {k_2}^2}{{a''}^4 \beta^8} (7 - 114\theta^2 + 395\theta^4) + \frac{5}{4} \frac{n'' k_4}{{a''}^4 \beta^8} (3 - 36\theta^2 + 49\theta^4)$$ $$\dot{\Omega}_2 = \frac{3}{2} \frac{n'' {k_2}^2}{{a''}^4 \beta^8} \theta (4 - 19\theta^2) + \frac{5}{2} \frac{n'' k_4}{{a''}^4 \beta^8} \theta (3 - 7\theta^2)$$ $$\dot{\ell} = n''_o + \dot{M}_1 + \dot{M}_2$$ $$\dot{\omega} = \dot{\omega}_1 + \dot{\omega}_2$$ $$\dot{\Omega} = \dot{\Omega}_1 + \dot{\Omega}_2$$ $$\xi = \frac{1}{a'' \beta^2 - s}$$ $$\eta = e s \xi$$ $$\psi = \sqrt{1 - \eta^2}$$ $$\alpha^2 = 1 + e^2$$ $$C_o = \frac{1}{2} B \rho_o (q_o - s)^4 n'' a'' \xi^4 \alpha^{-1} \psi^{-7}$$ $$C_1 = \frac{3}{2} n'' \alpha^4 C_o$$ $$D_1 = \xi \psi^{-2} / a'' \beta^2$$ $$D_2 = 12 + 36\eta^2 + \frac{9}{2} \eta^4$$ $$D_3 = 15\eta^2 + \frac{5}{2} \eta^4$$ $$D_4 = 5\eta + \frac{15}{4} \eta^3$$ $$D_5 = \xi \psi^{-2}$$ $$B_1 = -k_2 (1 - 3\theta^2)$$ $$B_2 = -k_2 (1 - \theta^2)$$ $$B_3 = \frac{A_{3,0}}{k_2} \sin i$$ $$C_2 = D_1 D_3 B_2$$ $$C_3 = D_4 D_5 B_3$$ $$\dot{n}_o = C_1 \left( 2 + 3\eta^2 + 20 e\eta + 5 e\eta^3 + \frac{17}{2} e^2 + 34 e^2 \eta^2 + D_1 D_2 B_1 + C_2 \cos 2\omega + C_3 \sin \omega \right)$$ $$\dot{e}_o = -\frac{2}{3} \frac{\dot{n}}{n''} (1 - e)$$

where all quantities are epoch values.

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

The secular effect of gravity is included in mean anomaly by

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

and the secular effects of gravity and atmospheric drag are included in argument of perigee and longitude of ascending node by

$$\omega = \omega_o + \dot{\omega}(t - t_o) + \dot{\omega}_1 Z_7$$ $$\Omega = \Omega_o + \dot{\Omega}(t - t_o) + \dot{\Omega}_1 Z_7$$

where

$$Z_7 = \frac{7}{3} Z_1 / n''_o$$

with

$$Z_1 = \frac{1}{2} \dot{n}_o (t - t_o)^2.$$

Next, SDP8 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).

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

$$n = n_{DS} + \dot{n}_o (t - t_o)$$ $$e = e_{DS} + \dot{e}_o (t - t_o)$$ $$M = M_{DS} + Z_1 + \dot{M}_1 Z_7$$

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

Here, SDP8 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.

Solve Kepler’s equation for $E$ by using the iteration equation

$$E_{i+1} = E_i + \Delta E_i$$

with

$$\Delta E_i = \frac{M + e \sin E_i - E_i}{1 - e \cos E_i}$$

and

$$E_1 = M + e \sin M + \frac{1}{2} e^2 \sin 2M.$$

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

$$a = \left( \frac{k_e}{n} \right)^{\frac{2}{3}}$$ $$\beta = (1 - e^2)^{\frac{1}{2}}$$ $$\sin f = \frac{\beta \sin E}{1 - e \cos E}$$ $$\cos f = \frac{\cos E - e}{1 - e \cos E}$$ $$u = f + \omega$$ $$r'' = \frac{a \beta^2}{1 + e \cos f}$$ $$\dot{r}'' = \frac{n a e}{\beta} \sin f$$ $$(r\dot{f})'' = \frac{n a^2 \beta}{r}$$ $$\delta r = \frac{1}{2} \frac{k_2}{a \beta^2} \left[ (1 - \theta^2) \cos 2u + 3(1 - 3\theta^2) \right] - \frac{1}{4} \frac{A_{3,0}}{k_2} \sin i_o \sin u$$ $$\delta \dot{r} = -n \left( \frac{a}{r} \right)^2 \left[ \frac{k_2}{a \beta^2} (1 - \theta^2) \sin 2u + \frac{1}{4} \frac{A_{3,0}}{k_2} \sin i_o \cos u \right]$$ $$\delta I = \theta \left[ \frac{3}{2} \frac{k_2}{a^2 \beta^4} \sin i_o \cos 2u - \frac{1}{4} \frac{A_{3,0}}{k_2 a \beta^2} e \sin \omega \right]$$ $$\delta(r\dot{f}) = -n \left( \frac{a}{r} \right)^2 \delta r + n a \left( \frac{a}{r} \right) \frac{\sin i_o}{\theta} \delta I$$ $$\delta u = \frac{1}{2} \frac{k_2}{a^2 \beta^4} \left[ \frac{1}{2} (1 - 7\theta^2) \sin 2u - 3(1 - 5\theta^2)(f - M + e \sin f) \right]$$ $$\quad - \frac{1}{4} \frac{A_{3,0}}{k_2 a \beta^2} \left[ \sin i_o \cos u (2 + e \cos f) + \frac{1}{2} \frac{\theta^2}{\sin i_o/2 \; \cos i_o/2} e \cos \omega \right]$$ $$\delta\lambda = \frac{1}{2} \frac{k_2}{a^2 \beta^4} \left[ \frac{1}{2} (1 + 6\theta - 7\theta^2) \sin 2u - 3(1 + 2\theta - 5\theta^2)(f - M + e \sin f) \right]$$ $$\quad + \frac{1}{4} \frac{A_{3,0}}{k_2 a \beta^2} \sin i_o \left[ \frac{e\theta}{1 + \theta} \cos \omega - (2 + e \cos f) \cos u \right]$$

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

$$r = r'' + \delta r$$ $$\dot{r} = \dot{r}'' + \delta \dot{r}$$ $$r\dot{f} = (r\dot{f})'' + \delta(r\dot{f})$$ $$y_4 = \sin \frac{I}{2} \sin u + \cos u \sin \frac{i_o}{2} \, \delta u + \frac{1}{2} \sin u \cos \frac{i_o}{2} \, \delta I$$ $$y_5 = \sin \frac{I}{2} \cos u - \sin u \sin \frac{i_o}{2} \, \delta u + \frac{1}{2} \cos u \cos \frac{i_o}{2} \, \delta I$$ $$\lambda = u + \Omega + \delta\lambda.$$

Unit orientation vectors are calculated by

$$U_x = 2 y_4 (y_5 \sin \lambda - y_4 \cos \lambda) + \cos \lambda$$ $$U_y = -2 y_4 (y_5 \cos \lambda + y_4 \sin \lambda) + \sin \lambda$$ $$U_z = 2 y_4 \cos \frac{I}{2}$$ $$V_x = 2 y_5 (y_5 \sin \lambda - y_4 \cos \lambda) - \sin \lambda$$ $$V_y = -2 y_5 (y_5 \cos \lambda + y_4 \sin \lambda) + \cos \lambda$$ $$V_z = 2 y_5 \cos \frac{I}{2}$$

where

$$\cos \frac{I}{2} = \sqrt{1 - {y_4}^2 - {y_5}^2}.$$

Position and velocity are given by

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

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

---

FORTRAN Implementation

*     SDP8                                            14 NOV 80
      SUBROUTINE SDP8(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
      DATA RHO/.15696615/
C
      IF (IFLAG .EQ. 0) GO TO 100
C
*     RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
*     FROM INPUT ELEMENTS --------- CALCULATE BALLISTIC COEFFICIENT
*     (B TERM) FROM INPUT B* DRAG TERM
C
      A1=(XKE/XNO)**TOTHRD
      COSI=COS(XINCL)
      THETA2=COSI*COSI
      TTHMUN=3.*THETA2-1.
      EOSQ=EO*EO
      BETAO2=1.-EOSQ
      BETAO=SQRT(BETAO2)
      DEL1=1.5*CK2*TTHMUN/(A1*A1*BETAO*BETAO2)
      AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
      DELO=1.5*CK2*TTHMUN/(AO*AO*BETAO*BETAO2)
      AODP=AO/(1.-DELO)
      XNODP=XNO/(1.+DELO)
      B=2.*BSTAR/RHO
C
*     INITIALIZATION
C
      PO=AODP*BETAO2
      POM2=1./(PO*PO)
      SINI=SIN(XINCL)
      SING=SIN(OMEGAO)
      COSG=COS(OMEGAO)
      TEMP=.5*XINCL
      SINIO2=SIN(TEMP)
      COSIO2=COS(TEMP)
      THETA4=THETA2**2
      UNM5TH=1.-5.*THETA2
      UNMTH2=1.-THETA2
      A3COF=-XJ3/CK2*AE**3
      PARDT1=3.*CK2*POM2*XNODP
      PARDT2=PARDT1*CK2*POM2
      PARDT4=1.25*CK4*POM2*POM2*XNODP
      XMDT1=.5*PARDT1*BETAO*TTHMUN
      XGDT1=-.5*PARDT1*UNM5TH
      XHDT1=-PARDT1*COSI
      XLLDOT=XNODP+XMDT1+
     2          .0625*PARDT2*BETAO*(13.-78.*THETA2+137.*THETA4)
      OMGDT=XGDT1+
     1    .0625*PARDT2*(7.-114.*THETA2+395.*THETA4)+PARDT4*(3.-36.*
     2        THETA2+49.*THETA4)
      XNODOT=XHDT1+
     1          (.5*PARDT2*(4.-19.*THETA2)+2.*PARDT4*(3.-7.*THETA2))*COSI
      TSI=1./(PO-S)
      ETA=EO*S*TSI
      ETA2=ETA**2
      PSIM2=ABS(1./(1.-ETA2))
      ALPHA2=1.+EOSQ
      EETA=EO*ETA
      COS2G=2.*COSG**2-1.
      D5=TSI*PSIM2
      D1=D5/PO
      D2=12.+ETA2*(36.+4.5*ETA2)
      D3=ETA2*(15.+2.5*ETA2)
      D4=ETA*(5.+3.75*ETA2)
      B1=CK2*TTHMUN
      B2=-CK2*UNMTH2
      B3=A3COF*SINI
      C0=.5*B*RHO*QOMS2T*XNODP*AODP*TSI**4*PSIM2**3.5/SQRT(ALPHA2)
      C1=1.5*XNODP*ALPHA2**2*C0
      C4=D1*D3*B2
      C5=D5*D4*B3
      XNDT=C1*(
     1  (2.+ETA2*(3.+34.*EOSQ)+5.*EETA*(4.+ETA2)+8.5*EOSQ)+
     1    D1*D2*B1+
     1    C4*COS2G+C5*SING)
      XNDTN=XNDT/XNODP
      EDOT=-TOTHRD*XNDTN*(1.-EO)
      IFLAG=0
      CALL DPINIT(EOSQ,SINI,COSI,BETAO,AODP,THETA2,SING,COSG,
     1          BETAO2,XLLDOT,OMGDT,XNODOT,XNODP)
C
*     UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
C
  100 Z1=.5*XNDT*TSINCE*TSINCE
      Z7=3.5*TOTHRD*Z1/XNODP
      XMAMDF=XMO+XLLDOT*TSINCE
      OMGASM=OMEGAO+OMGDT*TSINCE+Z7*XGDT1
      XNODES=XNODEO+XNODOT*TSINCE+Z7*XHDT1
      XN=XNODP
      CALL DPSEC(XMAMDF,OMGASM,XNODES,EM,XINC,XN,TSINCE)
      XN=XN+XNDT*TSINCE
      EM=EM+EDOT*TSINCE
      XMAM=XMAMDF+Z1+Z7*XMDT1
      CALL DPPER(EM,XINC,OMGASM,XNODES,XMAM)
      XMAM=FMOD2P(XMAM)
C
*     SOLVE KEPLERS EQUATION
C
      ZC2=XMAM+EM*SIN(XMAM)*(1.+EM*COS(XMAM))
      DO 130 I=1,10
      SINE=SIN(ZC2)
      COSE=COS(ZC2)
      ZC5=1./(1.-EM*COSE)
      CAPE=(XMAM+EM*SINE-ZC2)*
     1    ZC5+ZC2
      IF(ABS(CAPE-ZC2) .LE. E6A) GO TO 140
  130 ZC2=CAPE
C
*     SHORT PERIOD PRELIMINARY QUANTITIES
C
  140 AM=(XKE/XN)**TOTHRD
      BETA2M=1.-EM*EM
      SINOS=SIN(OMGASM)
      COSOS=COS(OMGASM)
      AXNM=EM*COSOS
      AYNM=EM*SINOS
      PM=AM*BETA2M
      G1=1./PM
      G2=.5*CK2*G1
      G3=G2*G1
      BETA=SQRT(BETA2M)
      G4=.25*A3COF*SINI
      G5=.25*A3COF*G1
      SNF=BETA*SINE*ZC5
      CSF=(COSE-EM)*ZC5
      FM=ACTAN(SNF,CSF)
      SNFG=SNF*COSOS+CSF*SINOS
      CSFG=CSF*COSOS-SNF*SINOS
      SN2F2G=2.*SNFG*CSFG
      CS2F2G=2.*CSFG**2-1.
      ECOSF=EM*CSF
      G10=FM-XMAM+EM*SNF
      RM=PM/(1.+ECOSF)
      AOVR=AM/RM
      G13=XN*AOVR
      G14=-G13*AOVR
      DR=G2*(UNMTH2*CS2F2G-3.*TTHMUN)-G4*SNFG
      DIWC=3.*G3*SINI*CS2F2G-G5*AYNM
      DI=DIWC*COSI
      SINI2=SIN(.5*XINC)
C
*     UPDATE FOR SHORT PERIOD PERIODICS
C
      SNI2DU=SINIO2*(
     1  G3*(.5*(1.-7.*THETA2)*SN2F2G-3.*UNM5TH*G10)-G5*SINI*CSFG*(2.+
     2        ECOSF))-.5*G5*THETA2*AXNM/COSIO2
      XLAMB=FM+OMGASM+XNODES+G3*(.5*(1.+6.*COSI-7.*THETA2)*SN2F2G-3.*
     1        (UNM5TH+2.*COSI)*G10)+G5*SINI*(COSI*AXNM/(1.+COSI)-(2.
     2        +ECOSF)*CSFG)
      Y4=SINI2*SNFG+CSFG*SNI2DU+.5*SNFG*COSIO2*DI
      Y5=SINI2*CSFG-SNFG*SNI2DU+.5*CSFG*COSIO2*DI
      R=RM+DR
      RDOT=XN*AM*EM*SNF/BETA+G14*(2.*G2*UNMTH2*SN2F2G+G4*CSFG)
      RVDOT=XN*AM**2*BETA/RM+
     1        G14*DR+AM*G13*SINI*DIWC
C
*     ORIENTATION VECTORS
C
      SNLAMB=SIN(XLAMB)
      CSLAMB=COS(XLAMB)
      TEMP=2.*(Y5*SNLAMB-Y4*CSLAMB)
      UX=Y4*TEMP+CSLAMB
      VX=Y5*TEMP-SNLAMB
      TEMP=2.*(Y5*CSLAMB+Y4*SNLAMB)
      UY=-Y4*TEMP+SNLAMB
      VY=-Y5*TEMP+CSLAMB
      TEMP=2.*SQRT(1.-Y4*Y4-Y5*Y5)
      UZ=Y4*TEMP
      VZ=Y5*TEMP
C
*     POSITION AND VELOCITY
C
      X=R*UX
      Y=R*UY
      Z=R*UZ
      XDOT=RDOT*UX+RVDOT*VX
      YDOT=RDOT*UY+RVDOT*VY
      ZDOT=RDOT*UZ+RVDOT*VZ
C
      RETURN
      END