Program Interface Issues

A few technical questions and comments are necessary to effectively integrate these analytical solutions into today’s environment.

A. Theoretical Issues

TLE data support a mix of coordinate systems and analytical theories. The SGP theory was largely based on Kozai (1959), while the SGP4 theory was primarily based on Brouwer (1959). The two theories are rather different, but both are still in use today. Neither the Kozai nor Brouwer theory originally included drag effects, so different treatments of atmospheric drag are in use. SGP approximates drag via rate changes of mean motion (Hilton and Kuhlman, 1966), while SGP4 uses power density functions (Lane and Cranford, 1969; Lane and Hoots, 1979) that require a term that encapsulates the ballistic coefficient, Bstar (see Vallado, 2004: 113—116). Simplified force modeling and the batch-least-squares processing of observational data often yield a Bstar that has “soaked up” force model errors. Occasionally, one finds negative Bstar values, indicating erroneously that energy is being added to the system, but this is simply a consequence of the limited SGP4 force modeling with respect to the actual dynamical environment.

B. Configuration Control

TLE data do not reveal which version of SGP4 was employed to estimate the orbital parameters. Different definitions of the so-called True Equator Mean Equinox (TEME) coordinate system and time systems may also have been used at different times. Without a list of dates to synchronize these changes with historical TLE data, the user must decide which version of the SGP4 propagator might be consistent. Because the accuracy of the propagator is generally in the kilometer-level range (Hartman, 1993), this may not be a problem for most cases, but as we’ll see shortly, some of the technical modifications can cause results to differ by hundreds of kilometers. This topic is perhaps the least likely to have a simple solution, but could potentially account for significant differences in ephemeris generation.

C. Data Formats

The TLE format appears to have changed slightly over the years, and numerous TLE data were disseminated with missing or erroneous values. Some of these cases simply test the error handling of the code and its ability to handle premature ending of the propagation.

The TLE Element Type is always set to zero for distributed data, although STR#3 suggests the following assignments: 1 = SGP, 2 = SGP4, 3 = SDP4, 4 = SGP8, 5 = SDP8. The TLE sets also use differing formats (e.g., use of leading zeros, or not). Sometimes, parameters are omitted within the TLE data (e.g., a second time derivative of mean motion or Bstar drag term equal to zero). These variations can confound fixed-read implementations in a computer program. The parsing of the TLE files is a bigger problem in languages such as C where the fixed-position approach (common in FORTRAN) is unusual, and where the ‘NUL’ (zeroth in the ASCII collating sequence) has a special end-of-string significance. Additionally, there are possible differences between DOS-formatted text files (CR/LF for end of line) and UNIX format (LF only). Attention is paid in the conversion utility to account for these discrepancies and the parsing routine is kept separate from the SGP4 routines to permit users the option of tailoring their parsing needs for a particular operation.

The TLE format has a simplistic form of error checking by having a checksum character for each line; however, it is prudent to check for other ‘fixed’ aspects (such as the “1” and “2” for each line, matching satellite numbers on the two lines, variable ranges, etc.) since the modulo-10 checksum only provides a 90% detection rate for uniformly random errors. Even with the checksum, there has been some ambiguity over the value assigned to the characters (the + sign in particular, which we believe should be zero), some additional explanation can be found on the web.[^6]

D. Coordinate System

The actual SGP4 model has little need for any specific coordinate or time system (e.g., the near-Earth part is rotationally symmetric about the pole), but when used for propagating TLE generated by DoD it becomes important to use the same coordinate system as the DoD orbit determination routines use. The commonly accepted output coordinate system is that of the “true equator, mean equinox” (TEME) (Herrick, 1971:325, 338, 341). An exact operational definition of TEME is very difficult to find in the literature, but conceptually its primary direction is related to the “uniform equinox” (Seidelmann, 1992:116, and Atkinson and Sadler, 1951). The intent was to provide an efficient, if approximate, coordinate system for use with the AFSPC analytical theories. Technically, the direction of the uniform equinox resides along the true equator “between” the origin of the intermediate Pseudo Earth Fixed (PEF) and True of Date (TOD) frames (Vallado, 2004:211, 221). It is found by observing that $\theta_{GAST82}$ may be separated into its components. Thus,

\vec{r}_{TOD} = ROT3(-\theta_{GAST82})\vec{r}_{PEF} \quad \text{and} \quad \theta_{GAST82} = \theta_{GMST82} + Eq_{e82}

\vec{r}_{TEME} = ROT3(-\theta_{GMST82})\vec{r}_{PEF}

\vec{r}_{PEF} = ROT3(\theta_{GMST82})\vec{r}_{TEME}

$\hfill (1)$

We recommend converting TEME to a truly standard coordinate frame before interfacing with other external programs. The preferred approach is to rotate to PEF using Greenwich Mean Sidereal Time (GMST), and then rotate to other standard coordinate frames. Conversions are well documented from this point. To implement, you simply apply a sidereal rotation about the Z-axis by GMST (using UT1 as we discuss later). Because polar motion has been historically neglected for General Perturbation (GP) applications, we assume that the pseudo Earth-fixed frame is the closest conventional frame.[^7]

If a rotation is made to TOD using the equation of the equinoxes, several approximations are introduced with the calculation of the nutation of the longitude ( $\Delta\Psi$ ) and the obliquity of the ecliptic ( $\epsilon$ ). There are at least three possible sources of uncertainty with this method: the number of terms to include in the nutation series, the inclusion of the post-1996 “kinematic correction” terms to the equation of the equinoxes, and small angle approximations. After choosing the length of the IAU 1980 nutation series (4, 10, and 106 terms are popular choices with 4 being most common), the transformation is sometimes further reduced by assuming that $\Delta\Psi \approx 0$ , $\epsilon \approx \bar{\epsilon}$ , and $\Delta\epsilon \approx 0$ . This results in a nutation matrix that is significantly simpler than the complete nutation matrix, although the complete form is more common today. The equation of the equinox may be approximated by ignoring the “kinematic correction” terms starting in 1997 [such that $Eq_{e_{1980}} \approx \Delta\Psi\cos(\epsilon)$ ]. Finally, because some of the multiplicative quantities are small, second-order terms may be neglected.

However, you should be aware of an additional nuance, specifically the ‘of date’ and ‘of epoch’ formulations.

TEME of Date — With this option, the epoch of the TEME frame is always the same as the epoch of the associated ephemeris generation time. The transformation to ECEF is done by first finding the conversion from TEME to TOD (third equation in Eq. (1)). Next the standard transformation from TOD to ECF is computed. We could have gone directly to PEF without the TOD frame (second equation in Eq. (1)), but this implementation enables comparison with the TEME of Epoch approach. All transformations are found using the complete IAU-76/FK5 formulae, including nutation.
TEME of Epoch — In this approach, the epoch of the TEME frame is held constant. Subsequent rotation matrices must therefore account for the change in precession and nutation from the epoch of the TEME frame to the epoch of the transformation. This is accomplished by finding a static transformation from TEME to J2000 — this includes the equation of the equinoxes, the nutation, and the precession which are all calculated at the epoch of the TLE. This static transformation is applied at each time requested in an ephemeris generation. Once the J2000 vector is found, standard techniques can convert this to other coordinate systems, at the appropriate time. This is computationally intensive, and introduces error into the subsequent solutions. All transformations, after the initial static calculations, are computed using the complete IAU-76/FK5 formulae, including all terms of the nutation theory.

Researchers generally believe the ‘of date’ option is correct, but confirmation from official sources is uncertain, and others infer that the ‘of epoch’ is correct. To be complete, we provide the equations and an example problem of both in the Appendix.

E. Time System Issues

Time accounting within SGP4 is referenced to the epoch of the TLE data. This practice makes individual satellite ephemeris generation and use relatively easy, although it can complicate multiple satellite analyses. The time system is assumed here to be UTC, but no formal documentation exists and UTC, as currently defined, was only introduced in 1972. UT1 is needed to calculate GMST for the coordinate transformations discussed in the appendix, but it is unknown whether UT1 or UTC is what is required by the software, although we assume UT1 for this paper. The error associated with approximating UT1 with UTC is within the theoretical uncertainty of the SGP4 theory itself. Except for the GMST calculation, this paper and code assumes time to be realized as UTC.

Time accounting also affects how the year of epoch values are handled within a system. This feature is only peripherally related to SGP4, and not part of the mathematical definition. It appears in the epoch calculations and affects how the two-digit year of the TLE is treated. Several possibilities exist. If the year is less than 50, 57, or some other value, one can add 2000, otherwise, 1900 is added. Of course, these are only temporary fixes with the correct option to be the use of a 4-digit year, Julian Date, Modified Julian Date, etc. During the so-called “Y2K” millennial rollover, some attention was focused here although nothing apparently changed.

It is doubtful that a leap-second capability was implemented into the peripheral software for SGP4 since the historical source code uses relative “time since epoch”. Any such addition is clearly outside the mathematical formulation of SGP4, but necessary for programs to interface with other agencies. As some software libraries have no support for the ‘61 second’ minute that is needed to properly represent or convert UTC time at the point of leap-second insertion, we suspect this is simply ignored for the majority of non-critical users, and a 1-second timing difference will occur sometime during the period between TLE updates where the leap second is added or subtracted. Although this is outside the direct scope of SGP4, it is part of many System Acceptance Tests for large programs, and is included here as a reminder of those operations.

F. GHA Calculation

The Greenwich Hour Angle is usually calculated using the Julian Date. However, you can also find expressions using the elapsed time from some epoch. Among the versions of SGP4 that are available today, several epochs arise: 1950 Jan 1 0^h, 1970 Jan 0 0^h, and 1970 Jan 1 12^h, UT1. The various combined constants illustrate the potential for error when using this approach. As new timing systems are developed, the associated timing parameters change slightly. The precision of these parameters also change slightly. Consider the following examples from various versions:

Jan 1, 1950 0 hr (original STR#3)
  THETA = 1.72944494D0 + 6.3003880987D0*DS50

Jan 0, 1970 0 hr
  C1    = 1.72027916940703639D-2
  THGR70 = 1.7321343856509374D0
  FK5R  = 5.07551419432269442D-15
  C1P2P = C1+TWOPI
  THGR  = DMOD(THGR70+C1*DS70+C1P2P*TFRAC+TS70*TS70*FK5R, twopi)

These approaches yield “essentially” the same values. A series of calculations were constructed to test these against the IAU convention (Vallado, 2004:191) using the Julian centuries of UT1 ( $T_{UT1}$ ).

\theta_{GMST1982} = 67{,}310.54841^s + (876{,}600^h + 8{,}640{,}184.812866^s)T_{UT1} + 0.093104 T_{UT1}^2 - 6.2 \times 10^{-6} T_{UT1}^3

$\hfill (2)$

The results showed comparisons of about $10^{-9}$ degrees difference. This is well below the level that the answers would be affected. We have included code for the conventional approach of Eq (2) for users wishing a more modern approach, but retained the older method of calculation for consistency with AFSPC.