\documentstyle[11pt,fleqn]{article}
\oddsidemargin 0in
\evensidemargin 0in
\textwidth 6.5in
\textheight 9in
\topmargin -36pt
\headheight 12pt
\headsep 24pt
\parindent 0pt
\def\dfrac#1#2{{\displaystyle{#1\over#2}}}
\def\dt{(t-t_o)}
\def\Epo{(E+\omega)}
\def\dddot#1{#1}
\def\dn{\dot{n}}
\def\dr{\dot{r}}
\def\dv{\dot{v}}
\def\df{\dot{f}}
\def\dM{\dot{M}}
\def\dC{\dot{C}}
\def\dD{\dot{D}}
\def\dom{\dot{\omega}}
\def\dOm{\dot{\Omega}}
\def\dal{\dot{\alpha}}
\def\dxi{\dot{\xi}}
\def\de{\dot{e}}
\def\deta{\dot{\eta}}
\def\dpsi{\dot{\psi}}
\def\ddn{\ddot{n}}
\def\dde{\ddot{e}}
\def\ddxi{\ddot{\xi}}
\def\ddeta{\ddot{\eta}}
\def\ddD{\ddot{D}}
\def\dddn{\stackrel{{\ldotp\ldotp\ldotp}}{n}}
\def\ppa{a''}
\def\ppn{n''}
\def\ppao{a''_o}
\def\ppno{n''_o}
\def\Mvec{{\bf M}}
\def\Nvec{{\bf N}}
\def\Uvec{{\bf U}}
\def\Vvec{{\bf V}}
\def\rvec{{\bf r}}
\def\drvec{\dot{\rvec}}
\def\L{{I\hspace{-3pt}L}}
\title{SPACETRACK REPORT NO. 3\\[12pt]
Models for Propagation of\\
NORAD Element Sets}
\author{ {\sc Felix R.\@ Hoots}\\
{\sc Ronald L.\@ Roehrich}\\[12pt]
{\sc December 1980}\\[12pt]
Package Compiled by\\
TS Kelso}
\date{31 December 1988}
\begin{document}
\pagenumbering{roman}
\addcontentsline{toc}{section}{ABSTRACT}
\maketitle

General perturbations element sets generated by NORAD can be used to predict
position and velocity of Earth-orbiting objects.  To do this one must be
careful to use a prediction method which is compatible with the way in which
the elements were generated.  Equations for five compatible models are given
here along with corresponding FORTRAN IV computer code.  With this information
a user will be able to make satellite predictions which are completely
compatible with NORAD predictions.

\newpage
\addcontentsline{toc}{section}{CONTENTS}
\tableofcontents
\newpage
APPROVED FOR PUBLIC RELEASE; DISTRIBUTION UNLIMITED

\vspace{1in}
Requests for additional copies by agencies of the Department of Defense, their
contractors, and other government agencies should be directed to
the:\footnote[1]{Editor's note: These are the addresses from the original document
and may no longer be valid.}

\begin{quote}
Defense Documentation Center\\
Cameron Station\\
Alexandria VA  22314
\end{quote}

All other persons and organizations should apply to the:

\begin{quote}
Department of Commerce\\
National Technical Information Service\\
5285 Port Royal Road\\
Springfield VA  22161
\end{quote}

Address special inquiries to:

\begin{quote}
Project Spacetrack Reports\\
Office of Astrodynamics\\
Aerospace Defense Center\\
ADC/DO6\\
Peterson AFB CO  80914
\end{quote}
\newpage
\parindent 2em
\parskip 12pt
\pagenumbering{arabic}
\section{INTRODUCTION}
NORAD maintains general perturbation element sets on all resident space
objects.  These element sets are periodically refined so as to maintain a
reasonable prediction capability on all space objects.  In turn, these element
sets are provided to users.  The purpose of this report is to provide the user
with a means of propagating these element sets in time to obtain a position
and velocity of the space object.

The \underline{most important} point to be noted is that not just any
prediction model will suffice.  The NORAD element sets are ``mean'' values
obtained by removing periodic variations in a particular way.  In order to
obtain good predictions, these periodic variations must be reconstructed (by
the prediction model) in exactly the same way they were removed by NORAD.
Hence, inputting NORAD element sets into a different model (even though the
model may be more accurate or even a numerical integrator) will result in
degraded predictions.  The NORAD element sets \underline{must} be used with
one of the models described in this report in order to retain maximum
prediction accuracy.

All space objects are classified by NORAD as near-Earth (period less than 225
minutes) or deep-space (period greater than or equal 225 minutes).  Depending
on the period, the NORAD element sets are automatically generated with the
near-Earth or deep-space model.  The user can then calculate the satellite
period and know which prediction model to use.

\section[The Propagation Models]{THE PROPAGATION MODELS}
Five mathematical models for prediction of satellite position and velocity are
available.  The first of these, SGP, was developed by Hilton \& Kuhlman (1966)
and is used for near-Earth satellites.  This model uses a simplification of
the work of Kozai (1959) for its gravitational model and it takes the drag
effect on mean motion as linear in time.  This assumption dictates a quadratic
variation of mean anomaly with time.  The drag effect on eccentricity is
modeled in such a way that perigee height remains constant.

The second model, SGP4, was developed by Ken Cranford in 1970 (see Lane and
Hoots 1979) and is used for near-Earth satellites.  This model was obtained by
simplification of the more extensive analytical theory of Lane and Cranford
(1969) which uses the solution of Brouwer (1959) for its gravitational model
and a power density function for its atmospheric model (see Lane, et al.\
1962).

The next model, SDP4, is an extension of SGP4 to be used for deep-space
satellites.  The deep-space equations were developed by Hujsak (1979) and
model the gravitational effects of the moon and sun as well as certain
sectoral and tesseral Earth harmonics which are of particular importance for
half-day and one-day period orbits.

The SGP8 model (see Hoots 1980) is used for near-Earth satellites and is
obtained by simplification of an extensive analytical theory of Hoots (to
appear) which uses the same gravitational and atmospheric models as Lane and
Cranford did but integrates the differential equations in a much different
manner.

Finally, the SDP8 model is an extension of SGP8 to be used for deep-space
satellites.  The deep-space effects are modeled in SDP8 with the same
equations used in SDP4.

\section[Compatibility with NORAD Element Sets]{COMPATIBILITY WITH NORAD ELEMENT SETS}
The NORAD element sets are currently generated with either SGP4 or SDP4
depending on whether the satellite is near-Earth or deep-space.  For element
sets sent to external users, the value of mean motion is altered slightly and
a pseudo-drag term ($\dn/2$) is generated.  These changes allow an SGP user
to make compatible predictions in the following manner.  If the satellite is
near-Earth, then the pseudo-drag term used in SGP simulates the drag effect of
the SGP4 model.  If the satellite is deep-space, then the pseudo-drag term
used in SGP simulates the deep-space secular effects of SDP4.

For SGP4 and SDP4 users, the mean motion is first recovered from its altered
form and the drag effect is obtained from the SGP4 drag term ($B^*$) with the
pseudo-drag term being ignored.  The value of the mean motion can be used to
determine whether the satellite is near-Earth or deep-space (and hence whether
SGP4 or SDP4 was used to generate the element set).  From this information the
user can decide whether to use SGP4 or SDP4 for propagation and hence be
assured of agreement with NORAD predictions.

The SGP8 and SDP8 models have the same gravitational and atmospheric models as
SGP4 and SDP4, although the form of the solution equations is quite different.
Additionally, SGP8 and SDP8 use a ballistic coefficient ($B$ term) in the drag
equations rather than the $B^*$ drag term.  However, compatible predictions
can be made with NORAD element sets by first calculating a $B$ term from the
SGP4 $B^*$ drag term.

At the present time consideration is being given to replacing SGP4 and SDP4 by
SGP8 and SDP8 as the NORAD satellite models.  In such a case the new NORAD
element sets would still give compatible predictions for SGP, SGP4, and SDP4
users and, for SGP8 and SDP8 users, would give agreement with NORAD
predictions.

\section[General Program Description]{GENERAL PROGRAM DESCRIPTION}
The five ephemeris packages cited in Section Two have each been programmed in
FORTRAN IV as stand-alone subroutines.  They each access the two function
subroutines ACTAN and FMOD2P and the deep-space equations access the function
subroutine THETAG.  The function subroutine ACTAN is a two argument (quadrant
preserving) arctangent subroutine which has been specifically designed to
return the angle within the range of 0 to $2\pi$.  The function subroutine
FMOD2P takes an angle and returns the modulo by $2\pi$ of that angle.  The
function subroutine THETAG calculates the epoch time in days since 1950 Jan
0.0 UTC, stores this in COMMON, and returns the right ascension of Greenwich
at epoch.

One additional subroutine DEEP is accessed by SDP4 and SDP8 to obtain the
deep-space perturbations to be added to the main equations of motion.

The main program DRIVER reads the input NORAD 2-line element set in either
G-card internal format or T-card transmission format and calls the appropriate
ephemeris package as specified by the user.  The DRIVER converts the elements
to the units of radians and minutes before calling the appropriate subroutine.
The ephemeris package returns position and velocity in units of Earth radii
and minutes.  These are converted by the DRIVER to kilometers and seconds for
printout.

All physical constants are contained in the constants COMMON C1 and can be
changed through the data statements in the DRIVER.  The one exception is the
physical constants used only in DEEP which are set in the data statements in
DEEP.

In the following sections the equations and program listing are given for each
ephemeris model.  Every effort has been made to maintain a strict parallel
structure between the equations and the computer code.

\section[The SGP Model]{THE SGP MODEL}
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(\dfrac{k_e}{n_o}\right)^{\frac23}\]
\[\delta_1=\dfrac34J_2\dfrac{a_E{}^2}{a_1{}^2}\dfrac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^{\frac32}}\]
\[a_o=a_1\left[1-\dfrac13\delta_1-\delta_1{}^2
-\dfrac{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\]
\[\dfrac{d\Omega}{dt}=-\dfrac32J_2\dfrac{a_E{}^2}{p_o{}^2}n_o\cos i_o\]
\[\dfrac{d\omega}{dt}=\dfrac34J_2\dfrac{a_E{}^2}{p_o{}^2}n_o(5\cos^2i_o-1).\]

The secular effects of atmospheric drag and gravitation are included through
the equations
\[a=a_o\left\{\dfrac{n_o}{n_o+2\left(\dfrac{\dn_o}2\right)
\dt+3\left(\dfrac{\ddn_o}6\right)\dt^2}\right\}^{\frac23}\]
\[e=\left\{\begin{array}{ll}
              1-\dfrac{q_o}a, & \mbox{for } a > q_o\\[12pt]
              10^{-6}, & \mbox{for } a \leq q_o
              \end{array}\right\}\]
\[p=a(1-e^2)\]
\[\Omega_{s_o}=\Omega_o+\dfrac{d\Omega}{dt}\dt\]
\[\omega_{s_o}=\omega_o+\dfrac{d\omega}{dt}\dt\]
\[L_s=L_o+\left(n_o+\dfrac{d\omega}{dt}+\dfrac{d\Omega}{dt}\right)
\dt+\dfrac{\dn_o}2\dt^2+\dfrac{\ddn_o}6\dt^3\]
where $\dt$ is time since epoch.

Long-period periodics are included through the equations
\[a_{yNSL}=e\sin\omega_{s_o}-\dfrac12\dfrac{J_3}{J_2}\dfrac{a_E}p\sin i_o\]
\[L=L_s-\dfrac14\dfrac{J_3}{J_2}\dfrac{a_E}pa_{xNSL}\sin i_o
\left[\dfrac{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
\[\Epo_{i+1}=\Epo_i+\Delta\Epo_i\]
with
\[\Delta\Epo_i=\dfrac{U-a_{yNSL}\cos\Epo_i+a_{xNSL}\sin\Epo_i-\Epo_i}
{-a_{yNSL}\sin\Epo_i-a_{xNSL}\cos\Epo_i+1}\]
\[U=L-\Omega_{s_o}\]
and
\[\Epo_1=U.\]
Then calculate the intermediate (partially osculating) quantities
\[e\cos E=a_{xNSL}\cos\Epo+a_{yNSL}\sin\Epo\]
\[e\sin E=a_{xNSL}\sin\Epo-a_{yNSL}\cos\Epo\]
\[e_L{}^2=(a_{xNSL})^2+(a_{yNSL})^2\]
\[p_L=a(1-e_L{}^2)\]
\[r=a(1-e\cos E)\]
\[\dr=k_e\dfrac{\sqrt{a}}re\sin E\]
\[r\dv=k_e\dfrac{\sqrt{p_L}}r\]
\[\sin u=\dfrac{a}r\left[\sin\Epo-a_{yNSL}-a_{xNSL}
\dfrac{e\sin E}{1+\sqrt{1-e_L{}^2}}\right]\]
\[\cos u=\dfrac{a}r\left[\cos\Epo-a_{xNSL}+a_{yNSL}
\dfrac{e\sin E}{1+\sqrt{1-e_L{}^2}}\right]\]
\[u=\tan^{-1}\left(\dfrac{\sin u}{\cos u}\right).\]

Short-period perturbations are now included by
\[r_k=r+\dfrac14J_2\dfrac{a_E{}^2}{p_L}\sin^2i_o\cos 2u\]
\[u_k=u-\dfrac18J_2\dfrac{a_E{}^2}{p_L{}^2}(7\cos^2i_o-1)\sin 2u\]
\[\Omega_k=\Omega_{s_o}+\dfrac34J_2\dfrac{a_E{}^2}{p_L{}^2}\cos i_o\sin 2u\]
\[i_k=i_o+\dfrac34J_2\dfrac{a_E{}^2}{p_L{}^2}\sin i_o\cos i_o\cos 2u.\]
Then unit orientation vectors are calculated by
\[\Uvec=\Mvec\sin u_k+\Nvec\cos u_k\]
\[\Vvec=\Mvec\cos u_k-\Nvec\sin u_k\]
where
\[\Mvec=\left\{\begin{array}{l}
                          M_x=-\sin\Omega_k\cos i_k\\
                          M_y=\cos\Omega_k\cos i_k\\
                          M_z=\sin i_k
                         \end{array}\right\}\]
\[\Nvec=\left\{\begin{array}{l}
                          N_x=\cos\Omega_k\\
                          N_y=\sin\Omega_k\\
                          N_z=0
                         \end{array}\right\}.\]
Then position and velocity are given by
\[\rvec=r_k\Uvec\]
and
\[\drvec=\dr\Uvec+(r\dv)\Vvec.\]

A FORTRAN IV computer code listing of the subroutine SGP is given below.
\newpage
\input{SGP.FOR}
\newpage
\section[The SGP4 Model]{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 ($\ppno$) and semimajor axis ($\ppao$) are
first recovered from the input elements by the equations
\[a_1=\left(\dfrac{k_e}{n_o}\right)^{\frac23}\]
\[\delta_1=\dfrac32\dfrac{k_2}{a_1{}^2}\dfrac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^{\frac32}}\]
\[a_o=a_1\left(1-\dfrac13\delta_1-\delta_1{}^2
-\dfrac{134}{81}\delta_1{}^3\right)\]
\[\delta_o=\dfrac32\dfrac{k_2}{a_o{}^2}\dfrac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^{\frac32}}\]
\[\ppno=\dfrac{n_o}{1+\delta_o}\]
\[\ppao=\dfrac{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^*=\ppao(1-e_o)-s+a_E\]
For perigee below 98 kilometers, the value of $s$ is changed to
\[s^*=20/\mbox{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]^{\frac14}+s-s^*\right]^4.\]
Then calculate the constants (using the appropriate values of $s$ and
$(q_o-s)^4$)
\[\theta=\cos i_o\]
\[\xi=\dfrac{1}{\ppao-s}\]
\[\beta_o=(1-e_o{}^2)^{\frac12}\]
\[\eta=\ppao e_o\xi\]
\[C_2=\begin{array}[t]{l}(q_o-s)^4\xi^4\ppno(1-\eta^2)^{-\frac72}
\left[\ppao\left(1+\dfrac32\eta^2+4e_o\eta+e_o\eta^3\right)\right.\\[12pt]
\left.+\dfrac32\dfrac{k_2\xi}{(1-\eta^2)}
\left(-\dfrac12+\dfrac32\theta^2\right)(8+24\eta^2+3\eta^4)\right]
\end{array}\]
\[C_1=B^*C_2\]
\[C_3=\dfrac{(q_o-s)^4\xi^5A_{3,0}\ppno a_E\sin i_o}{k_2e_o}\]
\[C_4=\begin{array}[t]{l}2\ppno(q_o-s)^4\xi^4\ppao\beta_o{}^2(1-\eta^2)^{-\frac72}
\biggl(\left[2\eta(1+e_o\eta)
+\dfrac12e_o+\dfrac12\eta^3\right]-\dfrac{2k_2\xi}{\ppao(1-\eta^2)}\times\\[12pt]
\left[3(1-3\theta^2)
\left(1+\dfrac32\eta^2-2e_o\eta-\dfrac12e_o\eta^3\right)
+\dfrac34(1-\theta^2)(2\eta^2-e_o\eta-e_o\eta^3)\cos 2\omega_o\right]\biggr)
\end{array}\]
\[C_5=2(q_o-s)^4\xi^4\ppao\beta_o{}^2(1-\eta^2)^{-\frac72}
\left[1+\dfrac{11}4\eta(\eta+e_o)+e_o\eta^3\right]\]
\[D_2=4\ppao\xi C_1{}^2\]
\[D_3=\dfrac43\ppao\xi^2(17\ppao+s)C_1{}^3\]
\[D_4=\dfrac23\ppao\xi^3(221\ppao+31s)C_1{}^4.\]

The secular effects of atmospheric drag and gravitation are included through
the equations
\[M_{DF}=M_o+\left[1+\dfrac{3k_2(-1+3\theta^2)}{2\ppao{}^2\beta_o{}^3}
+\dfrac{3k_2{}^2(13-78\theta^2+137\theta^4)}{16\ppao{}^4\beta_o{}^7}\right]
\ppno\dt\]
\[\omega_{DF}=\begin{array}[t]{l}\omega_o+\left[-\dfrac{3k_2(1-5\theta^2)}{2\ppao{}^2\beta_o{}^4}
+\dfrac{3k_2{}^2(7-114\theta^2+395\theta^4)}{16\ppao{}^4\beta_o{}^8}\right.\\[12pt]
\left.+\dfrac{5k_4(3-36\theta^2+49\theta^4)}{4\ppao{}^4\beta_o{}^8}\right]\ppno\dt
\end{array}\]
\[\Omega_{DF}=\Omega_o+\left[-\dfrac{3k_2\theta}{\ppao{}^2\beta_o{}^4}
+\dfrac{3k_2{}^2(4\theta-19\theta^3)}{2\ppao{}^4\beta_o{}^8}
+\dfrac{5k_4\theta(3-7\theta^2)}{2\ppao{}^4\beta_o{}^8}\right]\ppno\dt\]
\[\delta\omega=B^*C_3(\cos\omega_o)\dt\]
\[\delta M=-\dfrac23(q_o-s)^4B^*\xi^4\dfrac{a_E}{e_o\eta}
[(1+\eta\cos M_{DF})^3-(1+\eta\cos M_o)^3]\]
\[M_p=M_{DF}+\delta\omega+\delta M\]
\[\omega=\omega_{DF}-\delta\omega-\delta M\]
\[\Omega=\Omega_{DF}-\dfrac{21}2\dfrac{\ppno k_2\theta}{\ppao{}^2\beta_o{}^2}
C_1\dt^2\]
\[e=e_o-B^*C_4\dt-B^*C_5(\sin M_p-\sin M_o)\]
\[a=\ppao[1-C_1\dt-D_2\dt^2-D_3\dt^3-D_4\dt^4]^2\]
\[\L=\begin{array}[t]{l}
M_p+\omega+\Omega+\ppno\left[\dfrac32C_1\dt^2+(D_2+2C_1{}^2)\dt^3\right.\\[12pt]
+\dfrac14(3D_3+12C_1D_2+10C_1{}^3)\dt^4\\[12pt]
\left.+\dfrac15(3D_4+12C_1D_3+6D_2{}^2+30C_1{}^2D_2+15C_1{}^4)\dt^5\right]
\end{array}\]
\[\beta=\sqrt{(1-e^2)}\]
\[n=k_e\bigg/a^{\frac32}\]
where $\dt$ is time since epoch.  It should be noted that when epoch perigee
height is less than 220 kilometers, the equations for $a$ and $\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\]
\[\L_L=\dfrac{A_{3,0}\sin i_o}{8k_2a\beta^2}(e\cos\omega)
\left(\dfrac{3+5\theta}{1+\theta}\right)\]
\[a_{yNL}=\dfrac{A_{3,0}\sin i_o}{4k_2a\beta^2}\]
\[\L_T=\L+\L_L\]
\[a_{yN}=e\sin\omega+a_{yNL}.\]

Solve Kepler's equation for $\Epo$ by defining
\[U=\L_T-\Omega\]
and using the iteration equation
\[\Epo_{i+1}=\Epo_i+\Delta\Epo_i\]
with
\[\Delta\Epo_i =
\dfrac{U-a_{yN}\cos\Epo_i+a_{xN}\sin\Epo_i-\Epo_i}
{-a_{yN}\sin\Epo_i-a_{xN}\cos\Epo_i+1}\]
and
\[\Epo_1=U.\]

The following equations are used to calculate preliminary quantities needed
for short-period periodics.
\[e\cos E=a_{xN}\cos\Epo+a_{yN}\sin\Epo\]
\[e\sin E=a_{xN}\sin\Epo-a_{yN}\cos\Epo\]
\[e_L=(a_{xN}{}^2+a_{yN}{}^2)^{\frac12}\]
\[p_L=a(1-e_L{}^2)\]
\[r=a(1-e\cos E)\]
\[\dr=k_e\dfrac{\sqrt{a}}re\sin E\]
\[r\df=k_e\dfrac{\sqrt{p_L}}r\]
\[\cos u=\dfrac{a}r\left[\cos\Epo-a_{xN}+
\dfrac{a_{yN}(e\sin E)}{1+\sqrt{1-e_L{}^2}}\right]\]
\[\sin u=\dfrac{a}r\left[\sin\Epo-a_{yN}-
\dfrac{a_{xN}(e\sin E)}{1+\sqrt{1-e_L{}^2}}\right]\]
\[u=\tan^{-1}\left(\dfrac{\sin u}{\cos u}\right)\]
\[\Delta r=\dfrac{k_2}{2p_L}(1-\theta^2)\cos 2u\]
\[\Delta u=-\dfrac{k_2}{4p_L{}^2}(7\theta^2-1)\sin 2u\]
\[\Delta\Omega=\dfrac{3k_2\theta}{2p_L{}^2}\sin 2u\]
\[\Delta i=\dfrac{3k_2\theta}{2p_L{}^2}\sin i_o\cos 2u\]
\[\Delta\dr=-\dfrac{k_2n}{p_L}(1-\theta^2)\sin 2u\]
\[\Delta r\df=\dfrac{k_2n}{p_L}\left[(1-\theta^2)\cos 2u-\dfrac32(1-3
\theta^2)\right]\]

The short-period periodics are added to give the osculating quantities
\[r_k=r\left[1-\dfrac32k_2\dfrac{\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\]
\[\dr_k=\dr+\Delta\dr\]
\[r\df_k=r\df+\Delta r\df.\]
Then unit orientation vectors are calculated by
\[\Uvec=\Mvec\sin u_k+\Nvec\cos u_k\]
\[\Vvec=\Mvec\cos u_k-\Nvec\sin u_k\]
where
\[\Mvec=\left\{\begin{array}{l}
                          M_x=-\sin\Omega_k\cos i_k\\
                          M_y=\cos\Omega_k\cos i_k\\
                          M_z=\sin i_k
                         \end{array}\right\}\]
\[\Nvec=\left\{\begin{array}{l}
                          N_x=\cos\Omega_k\\
                          N_y=\sin\Omega_k\\
                          N_z=0
                         \end{array}\right\}.\]

Then position and velocity are given by
\[\rvec=r_k\Uvec\]
and
\[\drvec=\dr_k\Uvec+(r\df)_k\Vvec.\]

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.
\newpage
\input{SGP4.FOR}
\newpage
\section[The SDP4 Model]{THE SDP4 MODEL}
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 ($\ppno$) and semimajor axis ($\ppao$) are first
recovered from the input elements by the equations
\[a_1=\left(\dfrac{k_e}{n_o}\right)^{\frac23}\]
\[\delta_1=\dfrac32\dfrac{k_2}{a_1{}^2}\dfrac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^{\frac32}}\]
\[a_o=a_1\left(1-\dfrac13\delta_1-\delta_1{}^2-\dfrac{134}{81}\delta_1{}^3\right)\]
\[\delta_o=\dfrac32\dfrac{k_2}{a_o{}^2}\dfrac{(3\cos^2i_o-1)}{(1-e_o{}^2)^{\frac32}}\]
\[\ppno=\dfrac{n_o}{1+\delta_o}\]
\[\ppao=\dfrac{a_o}{1-\delta_o}.\]
For perigee between 98 kilometers and 156 kilometers, the value of the
constant $s$ used in SDP4 is changed to
\[s^*=\ppao(1-e_o)-s+a_E.\]
For perigee below 98 kilometers, the value of $s$ is changed to
\[s^*=20/\mbox{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]^{\frac14}+s-s^*\right]^4.\]
Then calculate the constants (using the appropriate values of $s$ and
$(q_o-s)^4$)
\[\theta=\cos i_o\]
\[\xi=\dfrac{1}{\ppao-s}\]
\[\beta_o=(1-e_o{}^2)^{\frac12}\]
\[\eta=\ppao e_o\xi\]
\[C_2=\begin{array}[t]{l}(q_o-s)^4\xi^4\ppno(1-\eta^2)^{-\frac72}
\left[\ppao(1+\dfrac32\eta^2+4e_o\eta+e_o\eta^3)\right.\\[12pt]
\left.+\dfrac32\dfrac{k_2\xi}{(1-\eta^2)}
\left(-\dfrac12+\dfrac32\theta^2\right)(8+24\eta^2+3\eta^4)\right]
\end{array}\]
\[C_1=B^*C_2\]
\[C_4=\begin{array}[t]{l}
2\ppno(q_o-s)^4\xi^4\ppao\beta_o{}^2(1-\eta^2)^{-\frac72}
\biggl(\left[2\eta(1+e_o\eta)+\dfrac12e_o+\dfrac12\eta^3\right]
-\dfrac{2k_2\xi}{\ppao(1-\eta^2)}\times\\[12pt]
\left[3(1-3\theta^2)\left(1+\dfrac32\eta^2-2e_o\eta-\dfrac12e_o\eta^3\right)
+\dfrac34(1-\theta^2)(2\eta^2-e_o\eta-e_o\eta^3)\cos 2\omega_o\right]\biggr)
\end{array}\]
\[\dM=\left[1+\dfrac{3k_2(-1+3\theta^2)}{2\ppao{}^2\beta_o{}^3}
+\dfrac{3k_2{}^2(13-78\theta^2+137\theta^4)}{16\ppao{}^4\beta_o{}^7}\right]\ppno\]
\[\dom = \left[-\dfrac{3k_2(1-5\theta^2)}{2\ppao{}^2\beta_o{}^4}
+\dfrac{3k_2{}^2(7-114\theta^2+395\theta^4)}{16\ppao{}^4\beta_o{}^8}
+\dfrac{5k_4(3-36\theta^2+49\theta^4)}{4\ppao{}^4\beta_o{}^8}\right]\ppno\]
\[\dOm_1=-\dfrac{3k_2\theta}{\ppao{}^2\beta_o{}^4}\ppno\]
\[\dOm=\dOm_1+\left[\dfrac{3k_2{}^2(4\theta-19\theta^3)}{2\ppao{}^4\beta_o{}^8}
+\dfrac{5k_4\theta(3-7\theta^2)}{2\ppao{}^4\beta_o{}^8}\right]\ppno.\]
At this point SDP4 calls the initialization section of DEEP which calculates
all initialized quantities needed for the deep-space perturbations (see
Section Ten).

The secular effects of gravity are included by
\[M_{DF}=M_o+\dM\dt\]
\[\omega_{DF}=\omega_o+\dom\dt\]
\[\Omega_{DF}=\Omega_o+\dOm\dt\]
where $\dt$ is time since epoch.  The secular effect of drag on longitude of
ascending node is included by
\[\Omega=\Omega_{DF}-\dfrac{21}2\dfrac{\ppno k_2\theta}{\ppao{}^2\beta_o{}^2}
C_1\dt^2.\]

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).

The secular effects of drag are included in the remaining elements by
\[a=a_{DS}[1-C_1\dt]^2\]
\[e=e_{DS}-B^*C_4\dt\]
\[\L=M_{DS}+\omega_{DS}+\Omega_{DS}+\ppno \left[
\dfrac32C_1\dt^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.

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.

Add the long-period periodic terms
\[a_{xN}=e\cos\omega\]
\[\beta=\sqrt{(1-e^2)}\]
\[\L_L=\dfrac{A_{3,0}\sin i_o}{8k_2a\beta^2}(e\cos\omega)
\left(\dfrac{3+5\theta}{1+\theta}\right)\]
\[a_{yNL}=\dfrac{A_{3,0}\sin i_o}{4k_2a\beta^2}\]
\[\L_T=\L+\L_L\]
\[a_{yN}=e\sin\omega+a_{yNL}.\]

Solve Kepler's equation for $\Epo$ by defining
\[U=\L_T-\Omega\]
and using the iteration equation
\[\Epo_{i+1}=\Epo_i+\Delta\Epo_i\]
with
\[\Delta\Epo_i =
\dfrac{U-a_{yN}\cos\Epo_i+a_{xN}\sin\Epo_i-\Epo_i}
{-a_{yN}\sin\Epo_i-a_{xN}\cos\Epo_i+1}\]
and
\[\Epo_1=U.\]

The following equations are used to calculate preliminary quantities needed
for short-period periodics.
\[e\cos E=a_{xN}\cos\Epo+a_{yN}\sin\Epo\]
\[e\sin E=a_{xN}\sin\Epo-a_{yN}\cos\Epo\]
\[e_L=(a_{xN}{}^2+a_{yN}{}^2)^{\frac12}\]
\[p_L=a(1-e_L{}^2)\]
\[r=a(1-e\cos E)\]
\[\dr=k_e\dfrac{\sqrt{a}}re\sin E\]
\[r\df=k_e\dfrac{\sqrt{p_L}}r\]
\[\cos u=\dfrac{a}r\left[\cos\Epo-a_{xN}+\dfrac{a_{yN}(e\sin E)}{1+\sqrt{1-e_L{}^2}}\right]\]
\[\sin u=\dfrac{a}r\left[\sin\Epo-a_{yN}-\dfrac{a_{xN}(e\sin E)}{1+\sqrt{1-e_L{}^2}}\right]\]
\[u=\tan^{-1}\left(\dfrac{\sin u}{\cos u}\right)\]
\[\Delta r=\dfrac{k_2}{2p_L}(1-\theta^2)\cos 2u\]
\[\Delta u=-\dfrac{k_2}{4p_L{}^2}(7\theta^2-1)\sin 2u\]
\[\Delta\Omega=\dfrac{3k_2\theta}{2p_L{}^2}\sin 2u\]
\[\Delta i=\dfrac{3k_2\theta}{2p_L{}^2}\sin i_o\cos 2u\]
\[\Delta\dr=-\dfrac{k_2n}{p_L}(1-\theta^2)\sin 2u\]
\[\Delta r\df=\dfrac{k_2n}{p_L}\left[(1-\theta^2)\cos 2u-\dfrac32(1-3\theta^2)\right]\]

The short-period periodics are added to give the osculating quantities
\[r_k=r\left[1-\dfrac32k_2\dfrac{\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\]
\[\dr_k=\dr+\Delta\dr\]
\[r\df_k=r\df+\Delta r\df.\]
Then unit orientation vectors are calculated by
\[\Uvec=\Mvec\sin u_k+\Nvec\cos u_k\]
\[\Vvec=\Mvec\cos u_k-\Nvec\sin u_k\]
where
\[\Mvec=\left\{\begin{array}{l}
                          M_x=-\sin\Omega_k\cos i_k\\
                          M_y=\cos\Omega_k\cos i_k\\
                          M_z=\sin i_k
                         \end{array}\right\}\]
\[\Nvec=\left\{\begin{array}{l}
                          N_x=\cos\Omega_k\\
                          N_y=\sin\Omega_k\\
                          N_z=0
                         \end{array}\right\}.\]

Then position and velocity are given by
\[\rvec=r_k\Uvec\]
and
\[\drvec=\dr_k\Uvec+(r\df)_k\Vvec.\]

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.
\newpage
\input{SDP4.FOR}
\newpage
\section[The SGP8 Model]{THE SGP8 MODEL}
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 ($\ppno$) and semimajor axis ($\ppao$) are first
recovered from the input elements by the equations
\[a_1=\left(\dfrac{k_e}{n_o}\right)^{\frac23}\]
\[\delta_1=\frac32\frac{k_2}{a_1{}^2}\frac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^\frac32}\]
\[a_o=a_1\left(1-\dfrac13\delta_1-\delta_1{}^2
-\dfrac{134}{81}\delta_1{}^3\right)\]
\[\delta_o=\dfrac32\dfrac{k_2}{a_o{}^2}\dfrac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^{\frac32}}\]
\[\ppno =\dfrac{n_o}{1+\delta_o}\]
\[\ppao=\dfrac{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})\mbox{ XKMPER kg/m$^2$/Earth radii}\]
is a reference value of atmospheric density.

Then calculate the constants
\[\beta^2=1-e^2\]
\[\theta=\cos i\]
\[\dM_1=-\dfrac32\dfrac{\ppn k_2}{\ppa^2\beta^3}(1-3\theta^2)\]
\[\dom_1=-\dfrac32\dfrac{\ppn k_2}{\ppa^2\beta^4}(1-5\theta^2)\]
\[\dOm_1=-3\dfrac{\ppn k_2}{\ppa^2\beta^4}\theta\]
\[\dM_2=\dfrac3{16}\dfrac{\ppn k_2{}^2}{\ppa^4\beta^7}(13-
78\theta^2+137\theta^4)\]
\[\dom_2=\dfrac3{16}\dfrac{\ppn k_2{}2}{\ppa^4\beta^8}(7-
114\theta^2+395\theta^4)+\dfrac54\dfrac{\ppn k_4}{\ppa^4\beta^8}(3-
36\theta^2+49\theta^4)\]
\[\dOm_2=\dfrac32\dfrac{\ppn k_2{}^2}{\ppa^4\beta^8}\theta(4-19\theta^2)
+\dfrac52\dfrac{\ppn k_4}{\ppa^4\beta^8}\theta(3-7\theta^2)\]
\[\dot\ell=\ppn+\dM_1+\dM_2\]
\[\dom=\dom_1+\dom_2\]
\[\dOm=\dOm_1+\dOm_2\]
\[\xi=\dfrac1{\ppa\beta^2-s}\]
\[\eta=es\xi\]
\[\psi=\sqrt{1-\eta^2}\]
\[\alpha^2=1+e^2\]
\[C_o=\dfrac12B\rho_o(q_o-s)^4\ppn\ppa\xi^4\alpha^{-1}\psi^{-7}\]
\[C_1=\dfrac32\ppn\alpha^4C_o\]
\[D_1=\xi\psi^{-2}/\ppa\beta^2\]
\[D_2=12+36\eta^2+\dfrac92\eta^4\]
\[D_3=15\eta^2+\dfrac52\eta^4\]
\[D_4=5\eta+\dfrac{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=\dfrac{A_{3,0}}{k_2}\sin i\]
\[C_2=D_1D_3B_2\]
\[C_3=D_4D_5B_3\]
\[\dn_o=C_1\left(2+3\eta^2+20e\eta+5e\eta^3+\dfrac{17}2e^2+34e^2\eta^2+D_1D_2B_1
+C_2\cos 2\omega+C_3\sin\omega\right)\]
\[C_4=D_1D_7B_2\]
\[C_5=D_5D_8B_3\]
\[D_6=30\eta+\dfrac{45}2\eta^3\]
\[D_7=5\eta+\dfrac{25}2\eta^3\]
\[D_8=1+\dfrac{27}4\eta^2+\eta^4\]
\[\de_o=-C_o\left(4\eta+\eta^3+5e+15e\eta^2+\dfrac{31}2e^2\eta+7e^2\eta^3+D_1D_6B_1
+C_4\cos 2\omega+C_5\sin\omega\right)\]
\[\dal/\alpha=e\de\alpha^{-2}\]
\[C_6=\dfrac13\dfrac{\dn}{\ppn}\]
\[\dxi/\xi=2\ppa\xi(C_6\beta^2+e\de)\]
\[\deta=(\de+e\dxi/\xi)s\xi\]
\[\dpsi/\psi=-\eta\deta\psi^{-2}\]
\[\dC_o/C_o=C_6+4\dxi/\xi-\dal/\alpha-7\dpsi/\psi\]
\[\dC_1/C_1=\dn/\ppn+4\dal/\alpha+\dC_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+\dfrac{45}4\eta^2\]
\[D_{14}=\dxi/\xi-2\dpsi/\psi\]
\[D_{15}=2(C_6+e\de\beta^{-2})\]
\[\dD_1=D_1(D_{14}+D_{15})\]
\[\dD_2=\deta D_{11}\]
\[\dD_3=\deta D_{12}\]
\[\dD_4=\deta D_{13}\]
\[\dD_5=D_5D_{14}\]
\[\dC_2=B_2(\dD_1D_3+D_1\dD_3)\]
\[\dC_3=B_3(\dD_5D_4+D_5\dD_4)\]
\[\dom=-\dfrac32\dfrac{\ppn k_2}{\ppa^2\beta^4}(1-5\theta^2)\]
\[D_{16}=D_9\deta+D_{10}\de+B_1(\dD_1D_2+D_1\dD_2)
+\dC_2\cos 2\omega+\dC_3\sin\omega+\dom(C_3\cos\omega-2C_2\sin 2\omega)\]
\[\ddn_o=\dn\dC_1/C_1+C_1D_{16}\]
\[\dde_o=\begin{array}[t]{l}
\de\dC_o/C_o-C_o\left\{\left(4+3\eta^2+30e\eta+\dfrac{31}2e^2+21e^2\eta^2\right)\deta
+(5+15\eta^2+31e\eta+14e\eta^3)\de\right.\\[12pt]
+B_1\left[\dD_1D_6+D_1\deta\left(30+\dfrac{135}2\eta^2\right)\right]
+B_2\left[\dD_1D_7+D_1\deta\left(5+\dfrac{75}2\eta^2\right)\right]\cos\omega\\[12pt]
\left.+B_3\left[\dD_5D_8+D_5\eta\deta\left(\dfrac{27}2+4\eta^2\right)\right]\sin\omega
+\dom(C_5\cos\omega-2C_4\sin 2\omega)\right\}
\end{array}\]
\[D_{17}=\ddn/\ppn-(\dn/\ppn)^2\]
\[\ddxi/\xi=2(\dxi/\xi-C_6)\dxi/\xi+2\ppa\xi\left(\dfrac13D_{17}\beta^2-2C_6e\de+\de^2+e\dde\right)\]
\[\ddeta=(\dde+2\de\dxi/\xi)s\xi+\eta\ddxi/\xi\]
\[D_{18}=\ddxi/\xi-(\dxi/\xi)^2\]
\[D_{19}=-(\dpsi/\psi)^2(1+\eta^{-2})-\eta\ddeta\psi^{-2}\]
\[\ddD_1=\dD_1(D_{14}+D_{15})+D_1\left(D_{18}-2D_{19}
+\dfrac23D_{17}+2\alpha^2\de^2\beta^{-4}+2e\dde\beta^{-2}\right)\]
\[\dddn_o=\begin{array}[t]{l}
\dn\left[\dfrac43D_{17}+3\de^2\alpha^{-2}+3e\dde\alpha^{-2}
-6(\dal/\alpha)^2+4D_{18}-7D_{19}\right]\\[12pt]
+\ddn\dC_1/C_1+C_1\biggl\{D_{16}\dC_1/C_1+D_9\ddeta+D_{10}\dde+\deta^2(6+30e\eta+68e^2)\\[12pt]
+\deta\de(40+30\eta^2+272e\eta)+\de^2(17+68\eta^2)\\[12pt]
+B_1[\ddD_1D_2+2\dD_1\dD_2+D_1(\ddeta D_{11}+\deta^2(72+54\eta^2))]\\[12pt]
+B_2[\ddD_1D_3+2\dD_1\dD_3+D_1(\ddeta D_{12}+\deta^2(30+30\eta^2))]\cos 2\omega\\[12pt]
+B_3\left[(\dD_5D_{14}+D_5(D_{18}-2D_{19}))D_4+2\dD_4\dD_5+D_5\left(\ddeta D_{13}+\dfrac{45}2\eta\deta^2\right)\right]\sin\omega\\[12pt]
+\dom[(7C_6+4e\de\beta^{-2})(C_3\cos\omega-2C_2\sin 2\omega)+2C_3\cos\omega\\[12pt]
-4C_2\sin 2\omega-\dom(C_3\sin\omega+4C_2\cos 2\omega)]\biggr\}
\end{array}\]
\[p=\dfrac{2\ddn_o{}^2-\dn_o\dddn_o}{\ddn_o{}^2-\dn_o\dddn_o}\]
\[\gamma=-\dfrac{\dddn_o}{\ddn_o}\dfrac1{(p-2)}\]
\[n_D=\dfrac{\dn_o}{p\gamma}\]
\[q=1-\dfrac{\dde_o}{\de_o\gamma}\]
\[e_D=\dfrac{\de_o}{q\gamma}\]
where all quantities are epoch values.

The secular effects of atmospheric drag and gravitation are included by
\[n=\ppno+n_D[1-(1-\gamma\dt)^p]\]
\[e=e_o+e_D[1-(1-\gamma\dt)^q]\]
\[\omega=\omega_o+\dom_1\left[\dt+\dfrac73\dfrac1{\ppno}Z_1\right]+\dom_2\dt\]
\[\Omega=\Omega_o''+\dOm_1\left[\dt+\dfrac73\dfrac1{\ppno}Z_1\right]+\dOm_2\dt\]
\[M=M_o+\ppno\dt+Z_1+\dM_1\left[\dt+\dfrac73\dfrac1{\ppno}Z_1\right]+\dM_2\dt\]
where
\[Z_1=\dfrac{\dn_o}{p\gamma}\left\{\dt+\dfrac1{\gamma(p+1)}[(1-\gamma\dt)^{p+1}-1]\right\}.\]
If drag is very small ($\dfrac{\dn}{\ppno }$ less than $1.5 \times 10^{-6}$/min) then the
secular equations for $n$, $e$, and $Z_1$ should be replaced by
\[n=\ppno+\dn\dt\]
\[e=e''_o+\de\dt\]
\[Z_1=\dfrac12\dn_o\dt^2\]
where $\dt$ is time since epoch and where
\[\de=-\dfrac23\dfrac{\dn_o}{\ppno}(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=\dfrac{M+e\sin E_i-E_i}{1-e\cos E_i}\]
and
\[E_1=M+e\sin M+\dfrac12e^2\sin 2M.\]

The following equations are used to calculate preliminary quantities needed for
the short-period periodics.
\[a=\left(\dfrac{k_e}n\right)^{\frac23}\]
\[\beta=(1-e^2)^{\frac12}\]
\[\sin f=\dfrac{\beta\sin E}{1-e\cos E}\]
\[\cos f=\dfrac{\cos E-e}{1-e\cos E}\]
\[u=f+\omega\]
\[r''=\dfrac{a\beta^2}{1+e\cos f}\]
\[\dr''=\dfrac{nae}\beta\sin f\]
\[(r\df)''=\dfrac{na^2\beta}r\]
\[\delta r=\dfrac12\dfrac{k_2}{a\beta^2}[(1-\theta^2)\cos 2u
+3(1-3\theta^2)]-\dfrac14\dfrac{A_{3,0}}{k_2}\sin i_o\sin u\]
\[\delta\dr=-n\left(\dfrac ar\right)^2\left[\dfrac{k_2}{a\beta^2}(1-\theta^2)\sin 2u
+\dfrac14\dfrac{A_{3,0}}{k_2}\sin i_o\cos u\right]\]
\[\delta I=\theta\left[\dfrac32\dfrac{k_2}{a^2\beta^4}\sin i_o\cos 2u
-\dfrac14\dfrac{A_{3,0}}{k_2a\beta^2}e\sin\omega\right]\]
\[\delta(r\df)=-n\left(\dfrac ar\right)^2\delta r+na\left(\dfrac ar\right)\dfrac{\sin i_o}\theta\delta I\]
\[\delta u=\begin{array}[t]{l}
\dfrac12\dfrac{k_2}{a^2\beta^4}\left[\dfrac12(1-7\theta^2)\sin 2u-3(1-5\theta^2)(f-M+e\sin f)\right]\\[12pt]
-\dfrac14\dfrac{A_{3,0}}{k_2a\beta^2}\left[\sin i_o\cos u(2+e\cos f)
+\dfrac12\dfrac{\theta^2}{\sin i_o/2\:\cos i_o/2}e\cos\omega\right]
\end{array}\]
\[\delta\lambda=\begin{array}[t]{l}
\dfrac12\dfrac{k_2}{a^2\beta^4}\left[\dfrac12(1+6\theta-7\theta^2)\sin 2u
-3(1+2\theta-5\theta^2)(f-M+e\sin f)\right]\\[12pt]
+\dfrac14\dfrac{A_{3,0}}{k_2a\beta^2}\sin i_o
\left[\dfrac{e\theta}{1+\theta}\cos\omega-(2+e\cos f)\cos u\right]
\end{array}\]

The short-period periodics are added to give the osculating quantities
\[r=r''+\delta r\]
\[\dr=\dr''+\delta\dr\]
\[r\df=(r\df)''+\delta(r\df)\]
\[y_4=\sin \dfrac{i_o}2\sin u+\cos u\sin \dfrac{i_o}2\delta u+\dfrac12\sin u\cos \dfrac{i_o}2\delta I\]
\[y_5=\sin \dfrac{i_o}2\cos u-\sin u\sin \dfrac{i_o}2\delta u+\dfrac12\cos u\cos \dfrac{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 \dfrac{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 \dfrac{I}2\]
where
\[\cos \dfrac{I}2=\sqrt{1-y_4{}^2-y_5{}^2}.\]
Position and velocity are given by
\[\rvec=r\Uvec\]
\[\drvec=\dr\Uvec+r\df\Vvec.\]

A FORTRAN IV computer code listing of the subroutine SGP8 is given below.
\newpage
\input{SGP8.FOR}
\newpage
\section[The SDP8 Model]{THE SDP8 MODEL}
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 ($\ppno$) and semimajor axis ($\ppao$) are first
recovered from the input elements by the equations
\[a_1=\left(\dfrac{k_e}{n_o}\right)^{\frac23}\]
\[\delta_1=\frac32\frac{k_2}{a_1{}^2}\frac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^\frac32}\]
\[a_o=a_1\left(1-\dfrac13\delta_1-\delta_1{}^2
-\dfrac{134}{81}\delta_1{}^3\right)\]
\[\delta_o=\dfrac32\dfrac{k_2}{a_o{}^2}\dfrac{(3\cos^2i_o-1)}
{(1-e_o{}^2)^{\frac32}}\]
\[\ppno =\dfrac{n_o}{1+\delta_o}\]
\[\ppao=\dfrac{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})\mbox{ XKMPER kg/m$^2$/Earth radii}\]
is a reference value of atmospheric density.

Then calculate the constants
\[\beta^2=1-e^2\]
\[\theta=\cos i\]
\[\dM_1=-\dfrac32\dfrac{\ppn k_2}{\ppa^2\beta^3}(1-3\theta^2)\]
\[\dom_1=-\dfrac32\dfrac{\ppn k_2}{\ppa^2\beta^4}(1-5\theta^2)\]
\[\dOm_1=-3\dfrac{\ppn k_2}{\ppa^2\beta^4}\theta\]
\[\dM_2=\dfrac3{16}\dfrac{\ppn k_2{}^2}{\ppa^4\beta^7}(13-
78\theta^2+137\theta^4)\]
\[\dom_2=\dfrac3{16}\dfrac{\ppn k_2{}2}{\ppa^4\beta^8}(7-
114\theta^2+395\theta^4)+\dfrac54\dfrac{\ppn k_4}{\ppa^4\beta^8}(3-
36\theta^2+49\theta^4)\]
\[\dOm_2=\dfrac32\dfrac{\ppn k_2{}^2}{\ppa^4\beta^8}\theta(4-19\theta^2)
+\dfrac52\dfrac{\ppn k_4}{\ppa^4\beta^8}\theta(3-7\theta^2)\]
\[\dot\ell=\ppno+\dM_1+\dM_2\]
\[\dom=\dom_1+\dom_2\]
\[\dOm=\dOm_1+\dOm_2\]
\[\xi=\dfrac1{\ppa\beta^2-s}\]
\[\eta=es\xi\]
\[\psi=\sqrt{1-\eta^2}\]
\[\alpha^2=1+e^2\]
\[C_o=\dfrac12B\rho_o(q_o-s)^4\ppn\ppa\xi^4\alpha^{-1}\psi^{-7}\]
\[C_1=\dfrac32\ppn\alpha^4C_o\]
\[D_1=\xi\psi^{-2}/\ppa\beta^2\]
\[D_2=12+36\eta^2+\dfrac92\eta^4\]
\[D_3=15\eta^2+\dfrac52\eta^4\]
\[D_4=5\eta+\dfrac{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=\dfrac{A_{3,0}}{k_2}\sin i\]
\[C_2=D_1D_3B_2\]
\[C_3=D_4D_5B_3\]
\[\dn_o=C_1\left(2+3\eta^2+20e\eta+5e\eta^3+\dfrac{17}2e^2+34e^2\eta^2+D_1D_2B_1
+C_2\cos 2\omega+C_3\sin\omega\right)\]
\[\de_o=-\dfrac23\dfrac{\dn}{\ppn}(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\dt\]
and the secular effects of gravity and atmospheric drag are included in
argument of perigee and longitude of ascending node by
\[\omega=\omega_o+\dom\dt+\dom_1Z_7\]
\[\Omega=\Omega_o+\dOm\dt+\dOm_1Z_7\]
where
\[Z_7=\dfrac73Z_1/\ppno\]
with
\[Z_1=\dfrac12\dn_o\dt^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}+\dn_o\dt\]
\[e=e_{DS}+\de_o\dt\]
\[M=M_{DS}+Z_1+\dM_1Z_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=\dfrac{M+e\sin E_i-E_i}{1-e\cos E_i}\]
and
\[E_1=M+e\sin M+\dfrac12e^2\sin 2M.\]

The following equations are used to calculate preliminary quantities needed for
the short-period periodics.
\[a=\left(\dfrac{k_e}n\right)^{\frac23}\]
\[\beta=(1-e^2)^{\frac12}\]
\[\sin f=\dfrac{\beta\sin E}{1-e\cos E}\]
\[\cos f=\dfrac{\cos E-e}{1-e\cos E}\]
\[u=f+\omega\]
\[r''=\dfrac{a\beta^2}{1+e\cos f}\]
\[\dr''=\dfrac{nae}\beta\sin f\]
\[(r\df)''=\dfrac{na^2\beta}r\]
\[\delta r=\dfrac12\dfrac{k_2}{a\beta^2}[(1-\theta^2)\cos 2u
+3(1-3\theta^2)]-\dfrac14\dfrac{A_{3,0}}{k_2}\sin i_o\sin u\]
\[\delta\dr=-n\left(\dfrac ar\right)^2\left[\dfrac{k_2}{a\beta^2}(1-\theta^2)\sin 2u
+\dfrac14\dfrac{A_{3,0}}{k_2}\sin i_o\cos u\right]\]
\[\delta I=\theta\left[\dfrac32\dfrac{k_2}{a^2\beta^4}\sin i_o\cos 2u
-\dfrac14\dfrac{A_{3,0}}{k_2a\beta^2}e\sin\omega\right]\]
\[\delta(r\df)=-n\left(\dfrac ar\right)^2\delta r+na\left(\dfrac ar\right)\dfrac{\sin i_o}\theta\delta I\]
\[\delta u=\begin{array}[t]{l}
\dfrac12\dfrac{k_2}{a^2\beta^4}\left[\dfrac12(1-7\theta^2)\sin 2u-3(1-5\theta^2)(f-M+e\sin f)\right]\\[12pt]
-\dfrac14\dfrac{A_{3,0}}{k_2a\beta^2}\left[\sin i_o\cos u(2+e\cos f)
+\dfrac12\dfrac{\theta^2}{sin i_o/2\:\cos i_o/2}e\cos\omega\right]
\end{array}\]
\[\delta\lambda=\begin{array}[t]{l}
\dfrac12\dfrac{k_2}{a^2\beta^4}\left[\dfrac12(1+6\theta-7\theta^2)\sin 2u
-3(1+2\theta-5\theta^2)(f-M+e\sin f)\right]\\[12pt]
+\dfrac14\dfrac{A_{3,0}}{k_2a\beta^2}\sin i_o
\left[\dfrac{e\theta}{1+\theta}\cos\omega-(2+e\cos f)\cos u\right]
\end{array}\]

The short-period periodics are added to give the osculating quantities
\[r=r''+\delta r\]
\[\dr=\dr''+\delta\dr\]
\[r\df=(r\df)''+\delta(r\df)\]
\[y_4=\sin \dfrac{I}2\sin u+\cos u\sin \dfrac{i_o}2\delta u+\dfrac12\sin u\cos \dfrac{i_o}2\delta I\]
\[y_5=\sin \dfrac{I}2\cos u-\sin u\sin \dfrac{i_o}2\delta u+\dfrac12\cos u\cos \dfrac{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 \dfrac{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 \dfrac{I}2\]
where
\[\cos \dfrac{I}2=\sqrt{1-y_4{}^2-y_5{}^2}.\]
Position and velocity are given by
\[\rvec=r\Uvec\]
\[\drvec=\dr\Uvec+r\df\Vvec.\]

A FORTRAN IV computer code listing of the subroutine SDP8 is given below.
\newpage
\input{SDP8.FOR}
\newpage
\section[The Deep-Space Subroutine]{THE DEEP-SPACE SUBROUTINE}
The two deep-space models, SDP4 and SDP8, both access the subroutine DEEP to
obtain the deep-space perturbations of the six classical orbital elements.
The perturbation equations are quite extensive and will not be repeated here.
Rather, this section will concentrate on a general description of the flow
between the main program and the deep-space subroutines.  A specific listing
of the equations is available in Hujsak (1979) or Hujsak and Hoots (1977).

The first time the deep-space subroutine is accessed is during the
initialization portion of SDP4/ SDP8 and is via the entry DPINIT.  Through this
entry, certain constants already calculated in SDP4/ SDP8 are passed to the
deep-space subroutine which in turn calculates all initialized (time
independent) quantities needed for prediction in deep space.  Additionally, a
determination is made and flags are set concerning whether the orbit is
synchronous and whether the orbit experiences resonance effects.

The next access to the deep-space subroutine occurs during the secular update
portion of SDP4/ SDP8 and is via the entry DPSEC.  Through this entry, the
current secular values of the ``mean'' orbital elements are passed to the
deep-space subroutine which in turn adds the appropriate deep-space secular
and long-period resonance effects to these mean elements.

The last access to the deep-space subroutine occurs at the beginning of the
osculation portion (periodics application) of SDP4/SDP8 and is via the entry
DPPER.  Through this entry, the current values of the orbital elements are
passed to the deep-space subroutine which in turn adds the appropriate deep-
space lunar and solar periodics to the orbital elements.

During initialization the deep-space subroutine calls the function subroutine
THETAG to obtain the location of Greenwich at epoch and to convert epoch to
minutes since 1950.  All physical constants which are unique to the deep-space
subroutine are set via data statements in DEEP rather than being passed
through a COMMON.

A FORTRAN IV computer code listing of the subroutine DEEP is given below.
These equations contain all currently anticipated changes to the SCC
operational program.  These changes are scheduled for implementation in March,
1981.
\newpage
\input{DEEP.FOR}
\newpage
\section[Driver and Function Subroutines]{DRIVER AND FUNCTION SUBROUTINES}
The DRIVER controls the input and output function and the selection of the
model.  The input consists of a program card which specifies the model to be
used and the output times and either a G-card or T-card element set.

The DRIVER reads and converts the input elements to units of radians and
minutes.  These are communicated to the prediction model through the COMMON
E1.  Values of the physical and mathematical constants are set and
communicated through the COMMONs C1 and C2, respectively.

The program card indicates the mathematical model to be used and the start and
stop time of prediction as well as the increment of time for output.  These
times are in minutes since epoch.

In the interest of efficiency the DRIVER sets a flag (IFLAG) the first time
the model is called.  This flag tells the model to calculate all initialized
(time independent) quantities.  After initialization, the model subroutine
turns off the flag so that all subsequent calls only access the time dependent
part of the model.  This mode continues until another input case is
encountered.

The DRIVER takes the output from the mathematical model (communicated through
the COMMON E1) and converts it to units of kilometers and seconds for
printout.

The function subroutine ACTAN is passed the values of sine and cosine in that
order and it returns the angle in radians within the range of 0 to $2\pi$.
The function subroutine FMOD2P is passed an angle in radians and returns the
angle in radians within the range of 0 to $2\pi$.  The function subroutine
THETAG is passed the epoch time exactly as it appears on the input element
cards.\footnote[1]{If only one year digit is given (as on standard G-cards) the
program assumes the 80 decade.  This may be overridden by putting a 2 digit
year in columns 30--31 of the first G-card.}  The routine converts this time
to days since 1950 Jan 0.0 UTC, stores this in the COMMON E1, and returns the
right ascension of Greenwich at epoch (in radians).

FORTRAN IV computer code listings of the routines DRIVER, ACTAN, FMOD2P,
and THETAG are given below.
\newpage
\input{DRIVER.FOR}
\newpage
\input{ACTAN.FOR}
\newpage
\input{FMOD2P.FOR}
\newpage
\input{THETAG.FOR}
\newpage
\section[Users Guide, Constants, and Symbols]{USERS GUIDE, CONSTANTS, AND SYMBOLS}
The first input card is the program card.  The format is as follows:

\begin{verbatim}
               Column     Format     Description

                1         I1         Ephemeris program desired
                                         1 = SGP
                                         2 = SGP4
                                         3 = SDP4
                                         4 = SGP8
                                         5 = SDP8
                2-11      F10.0      Prediction start time
               12-21      F10.0      Prediction stop time
               22-31      F10.0      Time increment
\end{verbatim}

\noindent
All times are in minutes since epoch and can be positive or negative.  The
second and third input cards consist of either a 2-card transmission or 2-card
G type element set.  Either type can be used with the only condition being
that the two cards must be in the correct order.  For reference a format sheet
for the T-card and G-card element sets follows this section.

The values of the physical and mathematical constants used in the program are
given below.

\begin{tabbing}
Variable name \quad \= distance units/Earth radii \quad \= \kill
\underline{Variable name} \> \underline{Definition} \> \underline{Value}\\[12pt]
CK2 \> $\dfrac12J_2a_E{}^2$ \> 5.413080E-4\\[12pt]
CK4 \> $-\dfrac38J_4a_E{}^4$ \> .62098875E-6\\[12pt]
E6A \> 10$^{-6}$ \> 1.0 E-6\\[12pt]
QOMS2T \> $(q_o-s)^4\mbox{ (er)}^4$ \> 1.88027916E-9\\[12pt]
S \> $s\mbox{ (er)}$ \> 1.01222928\\[12pt]
TOTHRD \> 2/3 \> .66666667\\[12pt]
XJ3 \> $J_3$ \> -.253881E-5\\[12pt]
XKE \> $k_e\left(\dfrac{\mbox{er}}{\mbox{min}}\right)^{\frac32}$ \> .743669161E-1\\[12pt]
XKMPER \> kilometers/Earth radii \> 6378.135\\[12pt]
XMNPDA \> time units/day \> 1440.0\\[12pt]
AE \> distance units/Earth radii \> 1.0\\[12pt]
DE2RA \> radians/degree \> .174532925E-1\\[12pt]
PI \> $\pi$ \> 3.14159265\\[12pt]
PIO2 \> $\pi/2$ \> 1.57079633\\[12pt]
TWOPI \> $2\pi$ \> 6.2831853\\[12pt]
X3PIO2 \> $3\pi/2$ \> 4.71238898
\end{tabbing}

\noindent
where er = Earth radii.  Except for the deep-space models, all ephemeris
models are independent of units.  Thus, units input or output as well as
physical constants can be changed by making the appropriate changes in only
the DRIVER program.

Following is a list of symbols commonly used in this report.

\begin{quote}
$n_o$ = the SGP type ``mean'' mean motion at epoch \\[12pt]
$e_o$ = the ``mean'' eccentricity at epoch \\[12pt]
$i_o$ = the ``mean'' inclination at epoch \\[12pt]
$M_o$ = the ``mean'' mean anomaly at epoch \\[12pt]
$\omega_o$ = the ``mean'' argument of perigee at epoch \\[12pt]
$\Omega_o$ = the ``mean'' longitude of ascending node at epoch \\[12pt]
$\dn_o$ = the time rate of change of ``mean'' mean motion at epoch \\[12pt]
$\ddn_o$ = the second time rate of change of ``mean'' mean motion at epoch \\[12pt]
$B^*$ = the SGP4 type drag coefficient \\[12pt]
$k_e$ = \parbox[t]{5.25in}{$\sqrt{GM}$ where $G$ is Newton's universal
gravitational constant and $M$ is the mass of the Earth} \\[12pt]
$a_E$ = the equatorial radius of the Earth \\[12pt]
$J_2$ = the second gravitational zonal harmonic of the Earth \\[12pt]
$J_3$ = the third gravitational zonal harmonic of the Earth \\[12pt]
$J_4$ = the fourth gravitational zonal harmonic of the Earth \\[12pt]
$\dt$ = time since epoch \\[12pt]
$k_2 = \dfrac12J_2a_E{}^2$ \\[12pt]
$k_4 = -\dfrac38J_4a_E{}^4$ \\[12pt]
$A_{3,0} = -J_3a_E{}^3$ \\[12pt]
$q_o$ = parameter for the SGP4/SGP8 density function \\[12pt]
$s$ = parameter for the SGP4/SGP8 density function \\[12pt]
$B$ = \parbox[t]{5.25in}{$\dfrac12C_D\dfrac Am$, the ballistic coefficient for SGP8
where $C_D$ is a dimensionless drag coefficient and $A$ is the average cross-%
sectional area of the satellite of mass $m$} \\[12pt]
\end{quote}
\newpage
\section[Sample Test Cases]{SAMPLE TEST CASES}
For reference a sample test case is given for each of the five
models.\footnote[1]{The test cases were generated on a machine with 8
digits of accuracy.  After a one day prediction, the test cases have only 5 to
6 digits of accuracy.}  The input used was standard T-cards and the output is
given at 360 minute intervals in units of kilometers and seconds.

When implemented on a given computer, the accuracies with which the test cases
are duplicated will be dominated by the accuracy of the epoch mean motion.
If, after reading and converting, the epoch mean motion has an error $\Delta n
= j \times 10^{-k}$ radians/time, then the predicted positions at time t may
differ from the test cases by numbers on the order of

\[\Delta r=\Delta n\dt(6,378.135) \mbox{ kilometers}\]
\newpage
\begin{verbatim}
1 88888U          80275.98708465  .00073094  13844-3  66816-4 0    8
2 88888  72.8435 115.9689 0086731  52.6988 110.5714 16.05824518  105

SGP  TSINCE              X                Y                Z

      0.            2328.96594238   -5995.21600342    1719.97894287
    360.00000000    2456.00610352   -6071.94232177    1222.95977784
    720.00000000    2567.39477539   -6112.49725342     713.97710419
   1080.00000000    2663.03179932   -6115.37414551     195.73919105
   1440.00000000    2742.85470581   -6079.13580322    -328.86091614


                      XDOT             YDOT             ZDOT

                       2.91110113      -0.98164053      -7.09049922
                       2.67852119      -0.44705850      -7.22800565
                       2.43952477       0.09884824      -7.31889641
                       2.19531813       0.65333930      -7.36169147
                       1.94707947       1.21346101      -7.35499924
\end{verbatim}
\newpage
\begin{verbatim}
1 88888U          80275.98708465  .00073094  13844-3  66816-4 0    8
2 88888  72.8435 115.9689 0086731  52.6988 110.5714 16.05824518  105

SGP4 TSINCE              X                Y                Z

      0.            2328.97048951   -5995.22076416    1719.97067261
    360.00000000    2456.10705566   -6071.93853760    1222.89727783
    720.00000000    2567.56195068   -6112.50384522     713.96397400
   1080.00000000    2663.09078980   -6115.48229980     196.39640427
   1440.00000000    2742.55133057   -6079.67144775    -326.38095856


                      XDOT             YDOT             ZDOT

                       2.91207230      -0.98341546      -7.09081703
                       2.67938992      -0.44829041      -7.22879231
                       2.44024599       0.09810869      -7.31995916
                       2.19611958       0.65241995      -7.36282432
                       1.94850229       1.21106251      -7.35619372
\end{verbatim}
\newpage
\begin{verbatim}
1 11801U          80230.29629788  .01431103  00000-0  14311-1
2 11801  46.7916 230.4354 7318036  47.4722  10.4117  2.28537848

SDP4 TSINCE              X                Y                Z

      0.            7473.37066650     428.95261765    5828.74786377
    360.00000000   -3305.22537232   32410.86328125  -24697.17675781
    720.00000000   14271.28759766   24110.46411133   -4725.76837158
   1080.00000000   -9990.05883789   22717.35522461  -23616.89062501
   1440.00000000    9787.86975097   33753.34667969  -15030.81176758


                      XDOT             YDOT             ZDOT

                       5.10715413       6.44468284      -0.18613096
                      -1.30113538      -1.15131518      -0.28333528
                      -0.32050445       2.67984074      -2.08405289
                      -1.01667246      -2.29026759       0.72892364
                      -1.09425066       0.92358845      -1.52230928
\end{verbatim}
\newpage
\begin{verbatim}
1 88888U          80275.98708465  .00073094  13844-3  66816-4 0    8
2 88888  72.8435 115.9689 0086731  52.6988 110.5714 16.05824518  105

SGP8 TSINCE              X                Y                Z

      0.            2328.87265015   -5995.21289063    1720.04884338
    360.00000000    2456.04577637   -6071.90490722    1222.84086609
    720.00000000    2567.68383789   -6112.40881348     713.29282379
   1080.00000000    2663.49508667   -6115.18182373     194.62816810
   1440.00000000    2743.29238892   -6078.90783691    -329.73434067


                      XDOT             YDOT             ZDOT

                       2.91210661      -0.98353850      -7.09081554
                       2.67936245      -0.44820847      -7.22888553
                       2.43992555       0.09893919      -7.32018769
                       2.19525236       0.65453661      -7.36308974
                       1.94680957       1.21500109      -7.35625595
\end{verbatim}
\newpage
\begin{verbatim}
1 11801U          80230.29629788  .01431103  00000-0  14311-1
2 11801  46.7916 230.4354 7318036  47.4722  10.4117  2.28537848

SDP8 TSINCE              X                Y                Z

      0.            7469.47631836     415.99390792    5829.64318848
    360.00000000   -3337.38992310   32351.39086914  -24658.63037109
    720.00000000   14226.54333496   24236.08740234   -4856.19744873
   1080.00000000  -10151.59838867   22223.69848633  -23392.39770508
   1440.00000000    9420.08203125   33847.21875000  -15391.06469727


                      XDOT             YDOT             ZDOT

                       5.11402285       6.44403201      -0.18296110
                      -1.30200730      -1.15603013      -0.28164955
                      -0.33951668       2.65315416      -2.08114153
                      -1.00112480      -2.33532837       0.76987664
                      -1.11986055       0.85410149      -1.49506933
\end{verbatim}

\newpage
\section[Sample Implementation]{SAMPLE IMPLEMENTATION}
These FORTRAN IV routines have been implemented on a Honeywell-6000 series
computer.  This machine has a processing speed in the 1MIPS range and a 36 bit
floating point word providing 8 significant figures of accuracy in single
precision.  The information in the following table is provided to allow a
comparison of the relative size and speed of the different models\footnote[1]%
{The timing results are for the test cases in Section Thirteen with a one day
prediction.  Times may vary slightly with orbital characteristics and, for
deep-space satellites, with prediction interval.}.
\begin{verbatim}
                    core used     CPU time per call (milliseconds)
           Model     (words)          Initialize       Continue

           SGP          541               .8              2.7
           SGP4       1,041              1.9              2.5
           SDP4       3,095              5.1              3.6
           SGP8       1,601              1.8              2.2
           SDP8       3,149              5.4              3.2
\end{verbatim}

\newpage
\addcontentsline{toc}{section}{ACKNOWLEDGEMENTS}
\section*{ACKNOWLEDGEMENTS}
The authors would like to express our appreciation to Bob Morris for his
helpful comments, to Linda Crawford for supplying one of the function
subroutines, and to Don Larson for his help in validating the computer
programs.  Finally, we wish to acknowledge the skill and patience of Maria
Gallegos who accomplished the unenviable task of typing this report.

\newpage
\addcontentsline{toc}{section}{REFERENCES}
\begin{thebibliography}{99}
\bibitem{}
Brouwer, D., ``Solution of the Problem of Artificial Satellite Theory without
Drag'', {\em Astronomical Journal\/} {\bf 64}, 378--397, November 1959.
\bibitem{}
Hilton, C.G. and Kuhlman, J.R., ``Mathematical Models for the Space Defense
Center'', Philco-Ford Publication No. U-3871, 17--28, November 1966.
\bibitem{}
Hoots, F.R., ``A Short, Efficient Analytical Satellite Theory''.  {\em AIAA\/}
Paper No.\ 80-1659, August 1980.
\bibitem{}
Hoots, F.R., ``Theory of the Motion of an Artificial Earth Satellite'',
accepted for publication in {\em Celestial Mechanics}.
\bibitem{}
Hujsak, R.S., ``A Restricted Four Body Solution for Resonating Satellites with
an Oblate Earth'', {\em AIAA\/} Paper No.\ 79-136, June 1979.
\bibitem{}
Hujsak, R.S. and Hoots, F.R., ``Deep Space Perturbations Ephemeris
Generation'', Aerospace Defense Command Space Computational Center Program
Documentation, DCD 8, Section 3, 82--104, September 1977.
\bibitem{}
Kozai, Y., ``The Motion of a Close Earth Satellite'', {\em Astronomical
Journal\/} {\bf 64}, 367--377, November 1959.
\bibitem{}
Lane, M.H. and Cranford, K.H., ``An Improved Analytical Drag Theory for the
Artificial Satellite Problem'', {\em AIAA\/} Paper No.\ 69-925, August 1969.
\bibitem{}
Lane, M.H., Fitzpatrick, P.M., and Murphy, J.J., ``On the Representation of
Air Density in Satellite Deceleration Equations by Power Functions with
Integral Exponents'', Project Space Track Technical Report No.\
APGC-TDR-62-15, March 1962, Air Force Systems Command, Eglin AFB, FL.
\bibitem{}
Lane, M.H. and Hoots, F.R., ``General Perturbations Theories Derived from the
1965 Lane Drag Theory'', Project Space Track Report No.\ 2, December 1979,
Aerospace Defense Command, Peterson AFB, CO.
\end{thebibliography}
\end{document}
