The SGP8 Model

Source

Spacetrack Report No. 3, Hoots & Roehrich (1980), compiled by T.S. Kelso (1988) — Pages 34—49 (pages 31—46 in original numbering)

The NORAD mean element sets can be used for prediction with SGP8. 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 = 2B^*/\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'' + \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 = es\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 + 20e\eta + 5e\eta^3 + \frac{17}{2}e^2 + 34e^2\eta^2 + D_1 D_2 B_1 + C_2 \cos 2\omega + C_3 \sin \omega\right)$

$C_4 = D_1 D_7 B_2$

$C_5 = D_5 D_8 B_3$

$D_6 = 30\eta + \frac{45}{2}\eta^3$

$D_7 = 5\eta + \frac{25}{2}\eta^3$

$D_8 = 1 + \frac{27}{4}\eta^2 + \eta^4$

$\dot{e}_o = -C_o \left(4\eta + \eta^3 + 5e + 15e\eta^2 + \frac{31}{2}e^2\eta + 7e^2\eta^3 + D_1 D_6 B_1 + C_4 \cos 2\omega + C_5 \sin \omega\right)$

$\dot{\alpha}/\alpha = e\dot{e}\alpha^{-2}$

$C_6 = \frac{1}{3} \frac{\dot{n}}{n''}$

$\dot{\xi}/\xi = 2a''\xi(C_6\beta^2 + e\dot{e})$

$\dot{\eta} = (\dot{e} + e\dot{\xi}/\xi)s\xi$

$\dot{\psi}/\psi = -\eta\dot{\eta}\psi^{-2}$

$\dot{C}_o/C_o = C_6 + 4\dot{\xi}/\xi - \dot{\alpha}/\alpha - 7\dot{\psi}/\psi$

$\dot{C}_1/C_1 = \dot{n}/n'' + 4\dot{\alpha}/\alpha + \dot{C}_o/C_o$

$D_9 = 6\eta + 20e + 15e\eta^2 + 68e^2\eta$

$D_{10} = 20\eta + 5\eta^3 + 17e + 68e\eta^2$

$D_{11} = 72\eta + 18\eta^3$

$D_{12} = 30\eta + 10\eta^3$

$D_{13} = 5 + \frac{45}{4}\eta^2$

$D_{14} = \dot{\xi}/\xi - 2\dot{\psi}/\psi$

$D_{15} = 2(C_6 + e\dot{e}\beta^{-2})$

$\dot{D}_1 = D_1(D_{14} + D_{15})$

$\dot{D}_2 = \dot{\eta}D_{11}$

$\dot{D}_3 = \dot{\eta}D_{12}$

$\dot{D}_4 = \dot{\eta}D_{13}$

$\dot{D}_5 = D_5 D_{14}$

$\dot{C}_2 = B_2(\dot{D}_1 D_3 + D_1 \dot{D}_3)$

$\dot{C}_3 = B_3(\dot{D}_5 D_4 + D_5 \dot{D}_4)$

$\dot{\omega} = -\frac{3}{2} \frac{n'' k_2}{a''^2 \beta^4} (1 - 5\theta^2)$

$D_{16} = D_9 \dot{\eta} + D_{10} \dot{e} + B_1(\dot{D}_1 D_2 + D_1 \dot{D}_2) + \dot{C}_2 \cos 2\omega + \dot{C}_3 \sin \omega + \dot{\omega}(C_3 \cos \omega - 2C_2 \sin 2\omega)$

$\ddot{n}_o = \dot{n}\dot{C}_1/C_1 + C_1 D_{16}$

$\ddot{e}_o = \dot{e}\dot{C}_o/C_o - C_o \left\{\left(4 + 3\eta^2 + 30e\eta + \frac{31}{2}e^2 + 21e^2\eta^2\right)\dot{\eta} + (5 + 15\eta^2 + 31e\eta + 14e\eta^3)\dot{e}\right.$

$\quad + B_1\left[\dot{D}_1 D_6 + D_1 \dot{\eta}\left(30 + \frac{135}{2}\eta^2\right)\right] + B_2\left[\dot{D}_1 D_7 + D_1 \dot{\eta}\left(5 + \frac{75}{2}\eta^2\right)\right]\cos\omega$

$\quad + \left.B_3\left[\dot{D}_5 D_8 + D_5\eta\dot{\eta}\left(\frac{27}{2} + 4\eta^2\right)\right]\sin\omega + \dot{\omega}(C_5 \cos\omega - 2C_4 \sin 2\omega)\right\}$

$D_{17} = \ddot{n}/n'' - (\dot{n}/n'')^2$

$\ddot{\xi}/\xi = 2(\dot{\xi}/\xi - C_6)\dot{\xi}/\xi + 2a''\xi\left(\frac{1}{3}D_{17}\beta^2 - 2C_6 e\dot{e} + \dot{e}^2 + e\ddot{e}\right)$

$\ddot{\eta} = (\ddot{e} + 2\dot{e}\dot{\xi}/\xi)s\xi + \eta\ddot{\xi}/\xi$

$D_{18} = \ddot{\xi}/\xi - (\dot{\xi}/\xi)^2$

$D_{19} = -(\dot{\psi}/\psi)^2(1 + \eta^{-2}) - \eta\ddot{\eta}\psi^{-2}$

$\ddot{D}_1 = \dot{D}_1(D_{14} + D_{15}) + D_1\left(D_{18} - 2D_{19} + \frac{2}{3}D_{17} + 2\alpha^2\dot{e}^2\beta^{-4} + 2e\ddot{e}\beta^{-2}\right)$

$\dddot{n}_o = \dot{n}\left[\frac{4}{3}D_{17} + 3\dot{e}^2\alpha^{-2} + 3e\ddot{e}\alpha^{-2} - 6(\dot{\alpha}/\alpha)^2 + 4D_{18} - 7D_{19}\right]$

$\quad + \ddot{n}\dot{C}_1/C_1 + C_1\left\{D_{16}\dot{C}_1/C_1 + D_9\ddot{\eta} + D_{10}\ddot{e} + \dot{\eta}^2(6 + 30e\eta + 68e^2)\right.$

$\quad + \dot{\eta}\dot{e}(40 + 30\eta^2 + 272e\eta) + \dot{e}^2(17 + 68\eta^2)$

$\quad + B_1[\ddot{D}_1 D_2 + 2\dot{D}_1 \dot{D}_2 + D_1(\ddot{\eta}D_{11} + \dot{\eta}^2(72 + 54\eta^2))]$

$\quad + B_2[\ddot{D}_1 D_3 + 2\dot{D}_1 \dot{D}_3 + D_1(\ddot{\eta}D_{12} + \dot{\eta}^2(30 + 30\eta^2))]\cos 2\omega$

$\quad + B_3\left[(\dot{D}_5 D_{14} + D_5(D_{18} - 2D_{19}))D_4 + 2\dot{D}_4\dot{D}_5 + D_5\left(\ddot{\eta}D_{13} + \frac{45}{2}\eta\dot{\eta}^2\right)\right]\sin\omega$

$\quad + \dot{\omega}[(7C_6 + 4e\dot{e}\beta^{-2})(C_3 \cos\omega - 2C_2 \sin 2\omega) + 2C_3 \cos\omega$

$\quad\left. - 4C_2 \sin 2\omega - \dot{\omega}(C_3 \sin\omega + 4C_2 \cos 2\omega)]\right\}$

$p = \frac{2\ddot{n}_o^2 - \dot{n}_o\,\dddot{n}_o}{\ddot{n}_o^2 - \dot{n}_o\,\dddot{n}_o}$

$\gamma = -\frac{\dddot{n}_o}{\ddot{n}_o}\frac{1}{(p - 2)}$

$n_D = \frac{\dot{n}_o}{p\gamma}$

$q = 1 - \frac{\ddot{e}_o}{\dot{e}_o \gamma}$

$e_D = \frac{\dot{e}_o}{q\gamma}$

where all quantities are epoch values.

The secular effects of atmospheric drag and gravitation are included by

$n = n''_o + n_D[1 - (1 - \gamma(t - t_o))^p]$

$e = e_o + e_D[1 - (1 - \gamma(t - t_o))^q]$

$\omega = \omega_o + \dot{\omega}_1\left[(t - t_o) + \frac{7}{3}\frac{1}{n''_o}Z_1\right] + \dot{\omega}_2(t - t_o)$

$\Omega = \Omega''_o + \dot{\Omega}_1\left[(t - t_o) + \frac{7}{3}\frac{1}{n''_o}Z_1\right] + \dot{\Omega}_2(t - t_o)$

$M = M_o + n''_o(t - t_o) + Z_1 + \dot{M}_1\left[(t - t_o) + \frac{7}{3}\frac{1}{n''_o}Z_1\right] + \dot{M}_2(t - t_o)$

where

$Z_1 = \frac{\dot{n}_o}{p\gamma}\left\{(t - t_o) + \frac{1}{\gamma(p + 1)}[(1 - \gamma(t - t_o))^{p+1} - 1]\right\}.$

If drag is very small ($\frac{\dot{n}}{n''_o}$ less than $1.5 \times 10^{-6}$/min) then the secular equations for $n$, $e$, and $Z_1$ should be replaced by

$n = n''_o + \dot{n}(t - t_o)$

$e = e''_o + \dot{e}(t - t_o)$

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

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

$\dot{e} = -\frac{2}{3}\frac{\dot{n}_o}{n''_o}(1 - e_o).$

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{nae}{\beta}\sin f$

$(r\dot{f})'' = \frac{na^2\beta}{r}$

$\delta r = \frac{1}{2}\frac{k_2}{a\beta^2}[(1 - \theta^2)\cos 2u + 3(1 - 3\theta^2)] - \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 + na\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_o}{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_o}{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 = 2y_4(y_5\sin\lambda - y_4\cos\lambda) + \cos\lambda$

$U_y = -2y_4(y_5\cos\lambda + y_4\sin\lambda) + \sin\lambda$

$U_z = 2y_4\cos\frac{I}{2}$

$V_x = 2y_5(y_5\sin\lambda - y_4\cos\lambda) - \sin\lambda$

$V_y = -2y_5(y_5\cos\lambda + y_4\sin\lambda) + \cos\lambda$

$V_z = 2y_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 SGP8 is given below.

---

FORTRAN Implementation

*         SGP8                                            14 NOV 80
      SUBROUTINE SGP8(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/

      IF (IFLAG .EQ. 0) GO TO 100

*     RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
*     FROM INPUT ELEMENTS --------- CALCULATE BALLISTIC COEFFICIENT
*     (B TERM) FROM INPUT B* DRAG TERM

      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

*     INITIALIZATION

      ISIMP=0
      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+    C4*COS2G+C5*SING)
      XNDTN=XNDT/XNODP

*     IF DRAG IS VERY SMALL, THE ISIMP FLAG IS SET AND THE
*     EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN MEAN
*     MOTION AND QUADRATIC VARIATION IN MEAN ANOMALY

      IF(ABS(XNDTN*XMNPDA) .LT. 2.16E-3) GO TO 50
      D6=ETA*(30.+22.5*ETA2)
      D7=ETA*(5.+12.5*ETA2)
      D8=1.+ETA2*(6.75+ETA2)
      C8=D1*D7*B2
      C9=D5*D8*B3
      EDOT=-C0*(
     1  ETA*(4.+ETA2+EOSQ*(15.5+7.*ETA2))+EO*(5.+15.*ETA2)+
     1  D1*D6*B1 +
     1  C8*COS2G+C9*SING)
      D20=.5*TOTHRD*XNDTN
      ALDTAL=EO*EDOT/ALPHA2
      TSDTTS=2.*AODP*TSI*(D20*BETAO2+EO*EDOT)
      ETDT=(EDOT+EO*TSDTTS)*TSI*S
      PSDTPS=-ETA*ETDT*PSIM2
      SIN2G=2.*SING*COSG
      C0DTC0=D20+4.*TSDTTS-ALDTAL-7.*PSDTPS
      C1DTC1=XNDTN+4.*ALDTAL+C0DTC0
      D9=ETA*(6.+68.*EOSQ)+EO*(20.+15.*ETA2)
      D10=5.*ETA*(4.+ETA2)+EO*(17.+68.*ETA2)
      D11=ETA*(72.+18.*ETA2)
      D12=ETA*(30.+10.*ETA2)
      D13=5.+11.25*ETA2
      D14=TSDTTS-2.*PSDTPS
      D15=2.*(D20+EO*EDOT/BETAO2)
      D1DT=D1*(D14+D15)
      D2DT=ETDT*D11
      D3DT=ETDT*D12
      D4DT=ETDT*D13
      D5DT=D5*D14
      C4DT=B2*(D1DT*D3+D1*D3DT)
      C5DT=B3*(D5DT*D4+D5*D4DT)
      D16=
     1  D9*ETDT+D10*EDOT +
     1  B1*(D1DT*D2+D1*D2DT) +
     1  C4DT*COS2G+C5DT*SING+XGDT1*(C5*COSG-2.*C4*SIN2G)
      XNDDT=C1DTC1*XNDT+C1*D16
      EDDOT=C0DTC0*EDOT-C0*(
     1  (4.+3.*ETA2+30.*EETA+EOSQ*(15.5+21.*ETA2))*ETDT+(5.+15.*ETA2
     '  +EETA*(31.+14.*ETA2))*EDOT +
     1  B1*(D1DT*D6+D1*ETDT*(30.+67.5*ETA2))                    +
     1  B2*(D1DT*D7+D1*ETDT*(5.+37.5*ETA2))*COS2G+
     1  B3*(D5DT*D8+D5*ETDT*ETA*(13.5+4.*ETA2))*SING+XGDT1*(C9*
     '  COSG-2.*C8*SIN2G))
      D25=EDOT**2
      D17=XNDDT/XNODP-XNDTN**2
      TSDDTS=2.*TSDTTS*(TSDTTS-D20)+AODP*TSI*(TOTHRD*BETAO2*D17-4.*D20*
     '  EO*EDOT+2.*(D25+EO*EDDOT))
      ETDDT =(EDDOT+2.*EDOT*TSDTTS)*TSI*S+TSDDTS*ETA
      D18=TSDDTS-TSDTTS**2
      D19=-PSDTPS**2/ETA2-ETA*ETDDT*PSIM2-PSDTPS**2
      D23=ETDT*ETDT
      D1DDT=D1DT*(D14+D15)+D1*(D18-2.*D19+TOTHRD*D17+2.*(ALPHA2*D25
     '  /BETAO2+EO*EDDOT)/BETAO2)
      XNTRDT=XNDT*(2.*TOTHRD*D17+3.*
     1  (D25+EO*EDDOT)/ALPHA2-6.*ALDTAL**2 +
     1  4.*D18-7.*D19 )                                          +
     1  C1DTC1*XNDDT+C1*(C1DTC1*D16+
     1  D9*ETDDT+D10*EDDOT+D23*(6.+30.*EETA+68.*EOSQ)+
     1  ETDT*EDOT*(40.+30.*
     '  ETA2+272.*EETA)+D25*(17.+68.*ETA2) +
     1  B1*(D1DDT*D2+2.*D1DT*D2DT+D1*(ETDDT*D11+D23*(72.+54.*ETA2))) +
     1  B2*(D1DDT*D3+2.*D1DT*D3DT+D1*(ETDDT*D12+D23*(30.+30.*ETA2))) *
     1  COS2G+
     1  B3*((D5DT*D14+D5*(D18-2.*D19)) *
     1 D4+2.*D4DT*D5DT+D5*(ETDDT*D13+22.5*ETA*D23)) *SING+XGDT1*
     1  ((7.*D20+4.*EO*EDOT/BETAO2)*
     '  (C5*COSG-2.*C4*SIN2G)
     '  +((2.*C5DT*COSG-4.*C4DT*SIN2G)-XGDT1*(C5*SING+4.*
     '  C4*COS2G))))
      TMNDDT=XNDDT*1.E9
      TEMP=TMNDDT**2-XNDT*1.E18*XNTRDT
      PP=(TEMP+TMNDDT**2)/TEMP
      GAMMA=-XNTRDT/(XNDDT*(PP-2.))
      XND=XNDT/(PP*GAMMA)
      QQ=1.-EDDOT/(EDOT*GAMMA)
      ED=EDOT/(QQ*GAMMA)
      OVGPP=1./(GAMMA*(PP+1.))
      GO TO 70
   50 ISIMP=1
      EDOT=-TOTHRD*XNDTN*(1.-EO)
   70 IFLAG=0

*     UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG

  100 XMAM=FMOD2P(XMO+XLLDOT*TSINCE)
      OMGASM=OMEGAO+OMGDT*TSINCE
      XNODES=XNODEO+XNODOT*TSINCE
      IF(ISIMP .EQ. 1) GO TO 105
      TEMP=1.-GAMMA*TSINCE
      TEMP1=TEMP**PP
      XN=XNODP+XND*(1.-TEMP1)
      EM=EO+ED*(1.-TEMP**QQ)
      Z1=XND*(TSINCE+OVGPP*(TEMP*TEMP1-1.))
      GO TO 108
  105 XN=XNODP+XNDT*TSINCE
      EM=EO+EDOT*TSINCE
      Z1=.5*XNDT*TSINCE*TSINCE
  108 Z7=3.5*TOTHRD*Z1/XNODP
      XMAM=FMOD2P(XMAM+Z1+Z7*XMDT1)
      OMGASM=OMGASM+Z7*XGDT1
      XNODES=XNODES+Z7*XHDT1

*     SOLVE KEPLERS EQUATION

      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

*     SHORT PERIOD PRELIMINARY QUANTITIES

  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

*     UPDATE FOR SHORT PERIOD PERIODICS

      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=SINIO2*SNFG+CSFG*SNI2DU+.5*SNFG*COSIO2*DI
      Y5=SINIO2*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

*     ORIENTATION VECTORS

      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

*     POSITION AND VELOCITY

      X=R*UX
      Y=R*UY
      Z=R*UZ
      XDOT=RDOT*UX+RVDOT*VX
      YDOT=RDOT*UY+RVDOT*VY
      ZDOT=RDOT*UZ+RVDOT*VZ

      RETURN
      END