Expected Code Updates

Although we searched many locations to obtain the latest openly available documentation on official AFSPC practice, a few topics remain unknown. The primary areas of discussion are those giving the largest differences in results — specifically negative inclinations, integrator problems, and solution of Kepler’s equation. If the reader is aware of other corrections, we would appreciate learning about them. The intention is to produce a new baseline that is as close as possible to the current operational version to enhance compatibility for the external user. While we could not verify these, we felt the changes were so obvious that AFSPC has already made them, thus we have included the options in the code. We used a comment (keyword “sgp4fix” in the codes) by each change to make any future retraction or addition easier. For official users who are constrained by the AFSPCI 33-105 restrictions and have only an executable version of the current code, it should be a simple matter to confirm these fixes.[^10]

A. Error Checking

We increased the amount of error checking in our code to handle cases such as decayed satellites, or satellites having inconsistent values. CelesTrak employs a significant amount of error checking on the TLE data, but programs allowing the user to enter data could result in values that would cause errors. Inclination values near 180.0 degrees can cause divide-by-zero problems in the initialization and the routine operation. This is fixed by setting a tolerance in both routines. The decay condition simply checks the position magnitude on each step.

B. Constants

Kaya et al. (2001, 2004) focuses on the difficulties encountered when mixing WGS-72 and WGS-84 constants. Because the SGP4 codes contain references to WGS-72, AFSPC may have updated the constants to WGS-84, but there is no other documentation supporting this so we present the development in case new official documentation is released. However, because many operational sites may still have embedded software containing a version of SGP4 using WGS-72, and the fact that the accuracy of the theory would not really be impacted, AFSPC may well have chosen to retain the older set of constants to better maintain interoperability with its internal resources. We use WGS-72 as the default value. As with other changes we discuss, this is only necessary to interface with external programs, but it will cause a difference in ephemeris results. The proper sequence to form the constants for WGS-72 is shown below. Note that we determined $\mu$ from the SGP4 code value of XKE because it is not specified directly in the code, and this makes future revisions easier. We also provide TUMin because XKE is simply the reciprocal of this quantity. TUMin is possibly more familiar as it is the number of minutes in one time unit — a necessary conversion when using canonical constants.

Table 2. WGS-72 Constants. The fundamental and derived constants are shown below. Notice that XKE and TUMin are reciprocal values. The original STR#3 listed XKE as 0.074 366 916.

Symbol	Calculation	Value
$\mu$		398,600.8 km^3/s^2
$R_\oplus$		6378.135 km
$J_2$		0.001 082 616
$J_3$		-0.000 002 538 81
$J_4$		-0.000 001 655 97
XKE	$60/\sqrt{R_\oplus^3/\mu}$	0.074 366 916 133 17 /min
TUMin	$\sqrt{R_\oplus^3/\mu}\,/\,60$	13.446 839 696 959 31 min

If we use WGS-84 values, we find the following values.

Table 3. WGS-84 Constants. The fundamental and derived constants are shown below. The zonal harmonic values are converted from the normalized values.

Symbol	Calculation	Value
$\mu$		398,600.4418 km^3/s^2
$R_\oplus$		6378.137 km
$J_2$	$C_{2,0} = -0.000\ 484\ 166\ 850\ 00$	0.001 082 629 989 05
$J_3$	$C_{3,0} = 0.000\ 000\ 957\ 063\ 90$	-0.000 002 532 153 06
$J_4$	$C_{4,0} = 0.000\ 000\ 536\ 995\ 87$	-0.000 001 610 987 61
XKE	$60/\sqrt{R_\oplus^3/\mu}$	0.074 366 853 168 71 /min
TUMin	$\sqrt{R_\oplus^3/\mu}\,/\,60$	13.446 851 082 044 98 min

Other constants may not be familiar at first. For example, XPDOTP is a conversion from rev/day to rad/min.

XPDOTP = 1440.0 / 2\pi = 229.183\ 118\ 052\ 329\ 3

RPTIM is simply the rotational velocity of the earth in rad/min. Note this does not use the GRS-80 defining parameter for the rotation of the Earth, $2\pi / (86{,}400/1.002\ 737\ 909\ 350\ 795) \times 60.0 = 7.292\ 115\ 855\ 3 \times 10^{-5}$ , but rather the GRS-67 value that Aoki et al. (1982) used in the definition of time.

RPTIM = 7.292\ 115\ 146\ 7 \times 10^{-5} \times 60.0 = 0.004\ 375\ 269\ 088\ 02\ \text{rad/min}

Other constants are combined with other values, or use the values mentioned previously in their formulation. We do not believe any update has occurred to any of the embedded constants in the deep-space portions as no documentation has ever suggested this.

C. Negative Inclination Orbits (Satellite 25954, 28626)

Deep-space orbits with low-inclination values (typically geosynchronous orbits) can, due to the effects of lunar and solar gravity, result in a negative inclination with time. This can create a step-function discontinuity in the positional components. Normally this is resolved by shifting the ascending node longitude by 180deg. In the computer code, we corrected this by removing the quadrant check from DSINIT before the ‘initialize resonance terms’ section, but kept the check in SGP4 before the ‘long period periodics’ section.

Satellites 25954 (at times beyond 274 minutes) and 28626 (at times beyond 1130 minutes) illustrate the effect of correcting negative inclination prematurely. The bottom graph in Fig. 5 of z-position reveals a discontinuity around this time in incorrect implementations of the code.

Figure 5. Negative Inclination Performance for Satellite 25954. The inclination and z-component of the position vector show the step function discontinuity of the previous SGP4 versions. The STR#3 and GSFC versions exhibits the problem while the Dundee version does not.

(Two plots: Inclination (deg) vs Min from Epoch, -1440 to +1440 min, showing Corrected and GSFC curves; z-Position Component (km) vs Min from Epoch, -1440 to +1440 min, showing Corrected and GSFC curves.)

D. Integrator Problems (Satellite 09880)

The original FORTRAN codes contained a generous mix of GOTOs and other structures that made accurate debugging nearly impossible. One area that appears to have suffered from this practice was the secular integrator used for 12^h and 24^h resonant cases. In particular, several satellites show difficulties when propagated ‘backwards,’ that is going to some time away from epoch (either positive or negative time) and then taking time steps towards the epoch again. The problem seems to be with the setup of the positive and negative steps (stepp and stepn) with values of 720 minutes in the DSPACE routine (SREZ in the older programs). It appears that ‘cleaning up’ the code has fixed the problem. The original STR#3 style of logic would integrate from epoch to the required time using a Taylor series approximation:

F(x+h) = F(x) + \frac{h}{1!}F'(x) + \frac{h^2}{2!}F''(x) + \cdots

$\hfill (3)$

where the integral at epoch $F(0)$ is defined as zero. If the time ( $h$ ) from epoch was greater than 720 minutes, it would step in 720 minute intervals ( $x = 720, 1440, \ldots$ ), recalculating the 1st and 2nd derivatives each time and saving the current (multiple of 720 minute) values for future use. Integrator resets occurred only when crossing the epoch. This was correct and efficient provided the model was only called with increasing time steps (either positive or negative), but it gives inconsistent results if you go to a time far from the epoch and return ‘backwards’ towards the epoch. The original code for this paper always integrated from the epoch to the required time, and restarted each time the model is called with a ‘backwards’ step. This was slightly less CPU efficient but led to repeatable results.[^11] Satellite 09998 demonstrates this symptom quite clearly as it is propagated from one day before the epoch until one day after the epoch, though all of the resonant and synchronous cases show it to some extent. Consider Fig. 6. We have made additional improvements to the DSPACE routine to enhance the CPU performance. For satellites with epochs far (decades) from the propagation time, the performance improvement is significant, and the accuracy is the same.

E. Solving Kepler’s Equation (Satellite 23333)

The partial fix discussed earlier handles a majority of the problem cases one would encounter in operations. However, additional robustness could be handled via several alternative methods. A simple but very effective fix for this is covered in Crawford (1995) where it is noted that the difference between mean and eccentric anomaly is never more than $\pm e$ radians, so if you limit the first Newton-Raphson correction to somewhere around this, it converges reliably for all cases. As this problem only applies to very high eccentricity orbits, an even simpler option fixes a limit of 0.9 — 1.0 for the maximum correction. Another option is that of Nijenhuis (1991), who examines the problem for eccentricities of 0.999 and 0.9999 and also examines the overall CPU load as well as the iteration count. Note that this iteration is not the ‘traditional’ iteration to find eccentric anomaly discussed in the literature (e.g., Vallado, 2004: 72—85). Results for this change were shown in Fig. 2.

A series of tests were run to determine the number of iterations for a complete satellite catalog, and the satellite tests we have included with this paper. For an example case of $e = 0.9$ the STR#3 version took an average of 5.685 iterations with a maximum of 8. The corrected Dundee version had an average of 3.984 iterations, with a maximum of 5. As with the Lyddane choice mentioned earlier, this change affects only a very small number of satellites.

Figure 6. Propagation Problems for Satellite 09998. The integrator problem is shown by looking at the semimajor axis (top) and the positional component differences (bottom). The scale is small, but the semimajor axis clearly shows the jump caused by incorrect integrator performance near 720 minutes prior to the epoch. The problem appears in all older versions.

(Two plots: Semimajor Axis (km) vs Min from Epoch, -1440 to -600 min, showing corrected and gsfc curves; Position Component Difference (km) vs Min from Epoch, -1440 to -600 min, showing x-, y-, z-components.)