Using SGP4 to obtain a spacecraft position & velocity from the NORAD 
(US SpaceCom) two-line orbit element (TLE) sets.

B. Knapp, 2006-01-19
Laboratory for Atmospheric and Space Physics
University of Colorado, Boulder
 
This is a rather large topic, and a good software system to provide orbit 
propagations (or interpolations) is a non-trivial development project.  
You really do not want to use SGP4 "as is", that is, as you obtain it from 
Celestrak or elsewhere.  I have developed an orbit position Fortran 
subroutine (spacecraft_pv.f) which uses my improved (Fortran 77) version
of SGP4, and a number of subsidiary routines which I also developed.  
Here is an outline of the call tree:

spacecraft_pv        // returns the position and velocity of a spacecraft
  get_lun            // get a logical unit number to open the TLE file
  nextTle            // read, check, and filter the TLE sets
     tleChk          // TLE check-sum check
  free_lun           // free the logical unit number
  TLE_Subset         // select a consistent subset of the TLE sets
     TLE_Filter      // internal consistency checking and filtering
        Neighbor_Consistency  // is a TLE set consistent with its neighbors?
           Pair_Difference    // "difference" (distance) between two TLE sets
              SGP4            // propagate a position and velocity
        Remove_set   // remove a TLE set from the array of all sets
  PV_from_Elts       // compute position, velocity from four adjacent TLE sets
     SGP4            // propagate a position and velocity
     ut2td           // convert universal time to dynamic time
       deltat        // the difference between UT and dynamic time
     Nutate_Matrix   // obtain the nutation in longitude and obliquity
        Obliq        // mean obliquity of the ecliptic
        Nutate       // compute nutation

There are several important features of this system, as compared with using 
SGP4 by itself.  First, this system is designed to supply *definitive* 
orbit positions, not *predictive*.  The latter are propagations forward in 
time (extrapolations) from the single most recent TLE set.  The uncertainties
in these extrapolations increase rapidly as the distance from the TLE epoch 
increases.  Another feature of spacecraft_pv is that it returns estimates 
of the uncertainties in the six elements of the spacecraft state vector 
(position, velocity).  Finally, spacecraft_pv provides some necessary 
management and filtering of the TLE sets--the raw TLE sets as retrieved from 
Space Track (www.space-track.org) have occasional problems that need to be 
detected and the problematic TLE sets discarded.

In order to provide "definitive" spacecraft ephemeris for a given time, 
we require the four bracketing TLE sets--that is, the two before and the
two following the given time.  Hence with this method predictive
ephemeris is not possible.  Then, of the four bracketing TLE sets, the
two TLE sets with adjacent epochs that are nearest in time (Euclidean
distance) to the given time are used to produce the propagation 
(interpolation).  Usually, but not always, the pair used is the bracketing 
pair.  For example, if there is a large gap (several days) between the 
middle two of the four TLE sets, then the first two or the last two may 
actually be closer in time to the requested time, than the bracketing pair.

To estimate the uncertainty in the result, we will propagate to the given
time (JD) from each of the two TLE sets with adjacent epochs that are 
nearest in time to the given JD. The error is a function of the difference 
Delta between these two propagations, and the distance in time from the 
nearest TLE set.  The best case is that the given JD coincides exactly 
with one of the two TLE epochs, in which case we estimate the uncertainty 
to be Delta/2 (the minimum).  If the given JD is equidistant from two 
bracketing TLE sets, we estimate the uncertainty to be Delta (a local 
maximum).  These conditions define a quartic: dy = (1-4*x^2*(1-2*x^2))*Delta, 
where x is the ratio of the distance in time from the midpoint of the two 
TLE epoch times to the difference in those epoch times, and y is one of the 
six position, velocity magnitudes.

I also provide an IDL call_external interface to spacecraft_pv (named
spacecraft_pv.pro), and if you are an IDL user, this is by far the easiest 
way to use this software.  I have a make file (for use on Sun Solaris with
the Sparc compilers) which will build the necessary shareable object library 
(orbit.so) for the IDL call_external call.  It will also build a Fortran 
demo program "xspacecraft_pv," which is a simple console program for one 
query at a time.

To install these routines on your computer, some amount of additional
work may be required.  With the exception of the driver program
xspacecraft_pv, all of the Fortran code should compile as is, and should 
not require any modification.  The program xspacecraft_pv.f should be 
edited to set the values of the variables path and pathlen to the name 
and name length of the path to your TLE file.

For Sun Solaris and the Sparc compilers, the only change needed to the 
Makefile should be the name of the environment variable (my_root) which 
identifies the root directory where your /lib directory resides where 
the orbit.so file will be built.  For other systems and compilers, more 
extensive modifications of the Makefile will be required to set the 
appropriate compiler switches, and to identify the locations of the math 
libraries needed by Fortran.  The help provided by RSI on building 
shareable objects (located in one of the samples directories in the IDL 
distribution) may be useful.

The IDL procedure spacecraft_pv.pro should be modified if you want the
default paths to the TLE file and the DLL (.so) files to be different
from the current working directory (so that they need not be specified
using the tle_path and dll_path keyword input parameters).

Please report any problems or comments to barry.knapp@lasp.colorado.edu.