Comparison Analyses
Many versions of SGP4 are available in code today, although most are initially from STR#3. Virtually none have been re-worked to restructure the code or to provide multiple computer programming languages and test results. Our aim is to correct that situation. Note that the basic structure of the computer code given in this paper has been available for several years in FORTRAN, Pascal, Ada, and C++ on the following web site (http://CelesTrak.com/software/vallado-sw.asp), although there has been extensive analysis to update the code for this paper. There are only three known “official” versions with which we could make comparisons. These include:
STR#3 (FORTRAN) Both the original single/double mix, and a double-precision version (just by adding the IMPLICIT DOUBLE statement) of STR#3 code were used. The electronic code was released to all users who asked for it. T. S. Kelso released an electronic package of the 1980 report in December 1988.
GSFC (FORTRAN) http://seawifs.gsfc.nasa.gov/SEAWIFS/SOFTWARE/src/bobdays/sgp4sub.f (original) Note this version is no longer available at this website although numerous downloads are known by organizations and countries. In addition, the code is still easily found on archive pages throughout the internet. A current site is similar, but potentially confusing as the subroutine name is the same, but the module is clearly labeled as a Brouwer—Lyddane model. http://www.icess.ucsb.edu/seawifs/seadas/src/utils/bobdays/sgp4sub.f (new, but different file)
JPL (FORTRAN) ftp://naif.jpl.nasa.gov/pub/naif/toolkit/FORTRAN/PC_Linux/packages/toolkit.tar.Z An additional source of SGP4 implementations is the JPL NAIF ‘spicelib’ toolkit, with source files ev2lin.f (basically SGP4.FOR equivalent), dpspce.f (basically SDP4.FOR) and zznrddp.f (basically the DEEP.FOR).
A few other codes were examined to determine what other researchers had done with the code. Examples that were tested but not included in the results presented here were:
ProjectPluto (C++) http://www.projectpluto.com/sat_code.htm Not an official version, but it is one of the more interesting and intelligent conversions to C++ available.
TrakStar (Pascal) http://CelesTrak.com/software/tskelso-sw.asp This is a very well known example, but is essentially a direct conversion of STR#3 and so it can be expected to behave in a manner similar to the STR#3 double-precision case.
Dundee (C) http://www.sat.dundee.ac.uk/~psc/sgp4.html The original translation into C by Paul Crawford and Andrew Brooks was virtually identical in behavior to STR#3, but with much better code structuring. Then many of the other fixes included such as the Kepler’s equation solution and secular integrator were added. The update over the last year for this paper had all of the corrections discussed and agreed with the authors, resulting in virtually identical results to the paper’s versions.
Because our version of SGP4 does not claim to be the official version, it was important to compare the results over a wide range of test conditions, and to compare with the released official versions. Specifically, the verification and stressing test cases provided a technical look at the performance, but these comparison tests were intended to show the robustness of the calculations under full-catalog simulations. Tests were run on several complete catalogs for varying dates. Each satellite was propagated from -1440 minutes to 1440 minutes at 20-minute time steps. The results were then compared between programs. The C++, FORTRAN, MATLAB, and Pascal versions gave virtually the same results, as shown in Fig. 7.
Figure 7. SGP4 Full-Catalog Comparisons. Maximum differences between ephemerides generated for two days are shown. Note the small scale for the C++ and FORTRAN comparison (top). The Pascal comparison (bottom) shows very small additional variations and these are from the 8-byte versus 10-byte precision in the language.
(Two scatter plots: Delta r (m) vs Period (min), 0—1600 min, log scale 0.00001—100000 m. Top: C++ vs FORTRAN (differences ~0.00001 m). Bottom: C++ vs Pascal (differences ~0.001 m).)
Comparisons were then run between the versions. Each figure shows the largest difference between the simulations, and each satellite is plotted against the orbital period. The scales are kept constant within each figure to permit rapid assessment of the differences. Figure 8 shows the results compared to the GSFC version. This was important to illustrate the similarity with the last known release.
Figure 8. SGP4 Full-Catalog Comparisons. Maximum differences between ephemerides generated for 2 days are shown. The top plot shows the paper C++ version against the GSFC code, while the bottom plot shows the comparison to the GSFC code, but only for propagations positive from the epoch. The differences are all related to the integrator problems before 720 minutes prior to epoch with geosynchronous and semi-synchronous orbits.
(Two scatter plots: Delta r (m) vs Period (min), 0—1600 min, log scale 0.00001—100000 m. Top: full propagation range showing up to ~100 m differences for resonant/synchronous orbits. Bottom: positive-only propagation showing sub-meter differences.)
As Fig. 8 shows, the GSFC version is very close to our revised version, and nearly identical to the performance between languages for the revised versions. The minor differences (usually a few meters) in the resonant cases (718-minute Molnyia and 1436-minute geostationary orbits) only show up with time steps that ‘go backwards’ in time (a problem in the secular integrator). In these tests, we begin at -1440 minutes and then step towards zero, before going ‘forwards’ towards +1440 minutes. In our revised version, the direction of propagation is not important. The GSFC version also has larger errors with the direct/Lyddane choice, and the inclination going negative during propagation, but these are not shown in Fig. 8 (it requires rare or ‘difficult’ TLEs to show up). Those differences were discussed earlier.
The comparisons with results from STR#3 show significantly larger differences for almost all satellites. Note that both (mixed) single and double-precision results are given in Fig. 9.
Figure 9. SGP4 Full-Catalog Comparisons. Maximum differences between ephemerides generated for 2 days are shown. The top plot shows the single-precision STR#3 version against the paper C++ version, while the bottom plot shows the STR#3 double-precision version comparison.
(Two scatter plots: Delta r (m) vs Period (min), 0—1600 min, log scale 0.00001—100000 m. Top: single-precision STR#3 showing ~100—10000 m differences. Bottom: double-precision STR#3 showing ~10—10000 m differences.)
Both versions of STR#3 (single/double and double only) show similar results, with agreement to reasonable accuracy (sub-km) for near-Earth orbits (limited by the precision of specifying the astrodynamic constants). This is less in deep space, where some of the limited precision (e.g., Kepler’s equation tolerance), the re-computing of perturbed terms, and the possibly the ‘Lyddane bug’ behavior show up more strongly. The differences between the GSFC version and the STR#3 versions are nearly identical to those in Fig. 9 for this catalog snapshot. This is important because many “correct” implementations of SGP4 in use are based on the STR#3 version (e.g., TrakStar), but this comparison shows the typical additional errors that users can expect as compared to the standalone AFSPC code.
Despite a rigorous attempt to review the fundamental constructs of the code, the JPL case is not particularly good in its original form. For near-Earth satellites, there is clearly some problem with the implementation (drag equations perhaps?) as it is much worse than STR#3 code in mixed precision. The deep-space cases have other problems, one of which is the choice to zero the LS offsets at epoch[^12] (which does not appear to be correctly implemented in any case). It also shares the negative inclination problem. These are unfortunate, as it is an interesting attempt to order and modernize the FORTRAN code, showing some insight into improving things, but missing others (such as the commonality of the SGP4/SDP4 codes) completely. The Project Pluto code (not shown) compared favorably to the revised code version but its lack of ‘official heritage’ made it difficult to accept changes based solely on its presence in this version. Results for the JPL code are shown in Fig. 10.
Figure 10. SGP4 Full-Catalog Comparisons. Maximum differences between ephemerides generated for 2 days are shown. The plot shows the paper C++ version against the original JPL code (plot on the top). Note that by changing the DOPERT variable in the JPL code, the results can be improved by about two orders of magnitude (plot on the bottom).
(Two scatter plots: Delta r (m) vs Period (min), 0—1600 min, log scale 0.00001—100000 m. Top: original JPL showing ~10—10000 m differences. Bottom: modified JPL showing ~0.01—1000 m differences.)