Computer Code Development

The revised computer code developed in this paper is provided in C++, FORTRAN, MATLAB, and Pascal to permit reasonable flexibility for applications (C++ is given in the appendix as this language is becoming commonplace). Conversion to other languages should be aided by the re-structuring effort that has been performed on the code.

There can be large variations between the numerous implementations of SGP4 — hence the need to establish a newer baseline that is compatible with CMOC as closely as possible to provide enhanced compatibility. Where obvious updates and corrections have been made and verified, we account for each in our revised code. For other improvements that appeared “obvious” to us, we tried to determine if these changes might be present in today’s standalone version.

Our starting point was the 1980 version from STR#3. From this point, STR#6, the Dundee modifications, and the GSFC code release verified several suspected code changes. There were too many changes in this update to describe them all — we list a few of the major ones below. Note that all satellite numbers refer to examples in the test case file (sgp4-all.tle). Also, all element plots are osculating values.

Merging of SGP4 and SDP4. A primary change from STR#3 was the merging of SGP4 and SDP4 code. A large number of researchers had noticed the commonality of the two models and simplified the code in this manner, however, not all had recognized and simplified the initialization code. Due to this simplification, most now refer to the merged SGP4/SDP4 models simply as ‘SGP4’.
Double-precision throughout. Although ultimately STR#6 has little relevance for TLE use, one notable change was the move to double-precision code throughout (rather than the mix of single and double in the original) and corresponding increase in accuracy for certain astrodynamic constants, all made practical by the improvement in computing power since STR#3. Such changes do not improve the “accuracy” of the model as such, but they lead to “smoother” behavior which helps with some tasks (such as differential correction), and to greater consistency of results on differing computer systems and/or compilers.
Kepler’s equation update. Solving Kepler’s equation was updated, but not completely fixed. The solution of Kepler’s equation continues to present challenges in astrodynamics hundreds of years after its introduction. The original 1980 version of SGP4 had a fixed limit of 10 iterations and a tolerance of $10^{-6}$ , but contained no code to prevent certain high-eccentricity orbits from failing to converge. Spacetrack Report Number 6 changed the tolerance to $10^{-12}$ (commensurate with double-precision work) and the GSFC version tried to solve the convergence problem by removing the iteration limit. However, these are incomplete approaches and can still result in infinite loops. The revised version code includes an updated SGP4 routine following the Dundee version that allows realistic controls on the iterations. Figure 1 shows the impact of this practice.
Lunar-solar term computation. The practice of only computing the lunar-solar terms if propagation time changes by more than 30 minutes to save CPU effort was dropped, thus resulting in smoother behavior for deep-space orbits with small time steps. This was the only function of the SAVTSN variable in the original DPPER subroutine. This resulted in ‘choppy’ behavior in some ephemerides from the STR#3 version.
Periodic lunar-solar perturbation updates. The application of periodic lunar-solar perturbations was updated. There are actually three problems relating to the application of the periodic lunar-solar perturbations:
1. The “Lyddane bug.” The first of these, sometimes known as the “Lyddane bug” (because it was first noted in independent investigations of the Lyddane modifications in DPPER), is due to the jump in the actan/atan2 output where the perturbed value due to the discontinuity of this function at either 90deg/270deg or +/-180deg (respectively). In the STR#3 code, the actan discontinuity occurred at 270deg. Spacetrack Report Number 6 tried using the atan2 function, but that simply moved the discontinuity to 180deg. The GSFC code (the IF statements at the end of the ‘apply periodics’ section in DPPER) confirmed the suspicions of several researchers about the need to evaluate the relative quadrant of the resulting angle and to correct accordingly. A similar problem exists with the modulo $2\pi$ reduction of the XNODE variable. The effect of not correcting the quadrant is illustrated in Fig. 2. This problem also occurs when intrinsic functions (mod, atan, etc.) are used instead of the STR#3 versions. We feel intrinsic functions are better suited for the program, but that full envelope testing of the Lyddane implementation is probably in order.
2. Initialization of deep-space terms on perturbed values. The second difficulty with the lunar-solar perturbations was the initialization of deep-space terms based on perturbed values. This was corrected in the DPPER and SGP4 routines of the Dundee and GSFC versions. In STR#3, the terms computed during initialization assumed fixed epoch values for inclination, etc., but of course they are perturbed by the deep-space terms. The approach used by the Dundee and GSFC versions includes any terms based on the Keplerian orbit being re-computed based on the new perturbed values.
3. The “Lyddane choice.” The third area of confusion with the lunar-solar perturbations is the decision for when to use the Lyddane modification, and we refer to it as the “Lyddane choice” (Satellites 14128, 20413). Lyddane (1963) reformulated the Brouwer expansions (done in Delaunay variables) in Poincare variables. Both are canonical, and the Poincare variables were intended to be non-singular for small eccentricity and inclination values. Because this was a reformulation, its use was intended for all computations, and remains that way today in the Navy Position Partials and Time (PPT3) (Hoots et al., 2004). During the development of SDP4 equations, some SIN(inclination) divisor problems were noted and the Lyddane formulation was examined. Because it also exhibited singularities, an alternate formulation was sought, ultimately resulting in new parameter choices. The decision was made to use these variables when the inclination was less than 11.4592deg (0.2 rad).

Figure 1. Solving Kepler’s Equation for Satellite 23333. The mean anomaly (bottom) illustrates the severe discrepancy in incorrectly solving Kepler’s equation after about 200 minutes. The effect also shows up in the inclination (top). The problem existed in the STR#3 version, shown here, but corrections were attempted in STR#6 and the GSFC version, with a better approach in the Dundee version. The inclination plot also shows the choppy (but smaller) behavior of the 30 minute updating of the lunar-solar terms in the STR#3 version before about 200 minutes. Notice that the effect goes away after about 1400 minutes.

(Two plots: Inclination (deg) vs Min from Epoch, 0—1600 min; Mean Anomaly (deg) vs Min from Epoch, 0—1600 min. Both show STR#3 and Corrected curves.)

Figure 2. Lunar-solar modifications for Satellite 23599. The argument of perigee and positional components illustrate the discrepancy in incorrectly updating the lunar-solar perturbations, and not accounting for the proper quadrant in the periodic calculations. The problem existed in the STR#3 version, shown here, and an attempted correction was made in STR#6, but it was mostly corrected in the GSFC version.

(Two plots: Argument of Perigee (deg) vs Min from Epoch, 250—750 min; Position Components (km) vs Min from Epoch, 250—750 min showing x-, y-, z-components with STR#3 and Corrected curves.)

Thus, the code implementation introduced two methods of applying the lunar-solar perturbations in the deep-space code, with the Lyddane modification used with smaller inclinations to avoid a divide-by-zero type of computation problem. In the STR#3 version, the choice was based on the unperturbed epoch inclination proximity to 11.4592deg. In STR#6 (and the GSFC code subroutine DPPER), the test used the perturbed inclination (XIP, with the secular term applied, unlike the XQNCL common term used in STR#3). This approach leads to a potential for the model switching lunar-solar methods as a function of propagation time, which is clearly undesirable. Note that the difference is usually small and relies on positional differences rather than the actual positional values. However, the results can be greatly magnified in some orbits (satellite 20413 which is a multi-day orbit). The basic effect can be demonstrated by satellite 14128 after about 2000 minutes from epoch when the perturbed inclination emerges above 11.4592deg as shown in Fig. 3.

Figure 3. Lyddane Choice Modification for Satellite 14128. The difference in positional components illustrates the discrepancy in applying the Lyddane modification using the perturbed inclination rather than the original inclination. The effect exists when comparing the GSFC and STR#3 versions. Note the differences diminish once the inclination crosses the 11.4592deg threshold (after 2000 min).

(Plot: Position Component Difference (km) vs Min from Epoch, 0—3000 min, showing y-component and z-component.)

A few comments are necessary on the Lyddane choice. This case occurs exceedingly infrequently as the satellite inclination must be very close to the 11.4592deg limit, and it must be a deep-space satellite. Without carefully crafted test cases, this (like several of the problems we discuss) is not easy to detect in normal operation. At the time of coding the STR#3 version, this form of software testing was not a commonly taught practice.[^8] We can consider several methods of resolving the situation, such as (a) going back to the STR#3 practice of testing the unperturbed epoch inclination at $t = 0$ and using that as a fixed decision for the TLE, (b) testing the perturbed inclination at each propagation time (as in the GSFC version), or (c) making more significant code changes to find a smoother way of blending the two lunar-solar perturbation methods. Although a rare case, we think some fix should be included and hope that AFSPC will confirm the current state of the code so users can be compatible in all cases. For the present paper, we have included Option (b) in the code with qualifiers to assist in location and potential future resolution. However, we note that there is probably a better “crossover” point to apply the Lyddane modification (Option c) that will not result in such large discrepancies in the two ephemerides, but time did not permit a thorough investigation and recommendation for such a change.

Next, the following changes were made to comply with modern programming standards, and to facilitate any changes in the future. With the exception of the variable precision and the integrator issues (discussed later), none of these changes affected the technical performance of the program and could be considered “cosmetic”.

Comprehensive variable declarations. Implicit typing in FORTRAN was replaced by comprehensive variable declarations. This was a critical step before conversion to C++ and others. Modern compilers can generally sort out the variable names, but the possibility of mistaken variables, variables being set to zero and used in calculations, etc., was too great. In addition, knowing which variables were calculated and set assisted the process of forming structures. Finally, this also eliminated much of the need for the FORTRAN SAVE command to hold values between function calls with certain compilers.
Structures for data passing. Structures were created to pass the large amounts of data between functions. Numerous variables were passed between functions in the original code. With no typing in the original code, this approach proved relatively easy, but it was difficult to gain an understanding of the underlying structure. The structures were set up to support integrated near-earth and deep-space functionality provided in the code. This change also supported processing multiple satellites at one time. While processing a single satellite is illustrative for simple scenarios, it is unrealistic for many modern applications. For instance, the SOCRATES effort uses TLE data to generate potential conjunction information. During these runs, one must have two or more satellites in memory at one time.
GOTO elimination. GOTO statements in FORTRAN were eliminated, using more modern constructs. This old programming construct is often seen in legacy programs, but completely unnecessary with modern programming techniques and tools. Looping and decision constructs were inserted, as appropriate.
Intrinsic functions. Intrinsic functions replace user-written routines. Trigonometric and exponential routines should use intrinsic calls within the programming language. The only exception should be in cases where a specific quadrant, ordering, etc. is required. The only case where this appears to be necessary is in DPPER where the MOD function modifies a variable that is then used outside a trigonometric expression. In this case, we opted to retain the modern MOD function, but simply add an IF statement to check if the result is less than zero. The ATAN2 function a few lines later may also require a similar modification, but this has a much smaller effect in the end results.
Initialization separation. Initialization functions were separated for better organization. The code was modularized, keeping initialization functions separate from routine function use. Although modern compilers can generally sort these differences out, the code is easier to maintain if the functions are isolated for a particular operation. The reorganization of the computer code simplified the processing flow. In addition, simple timing studies performed during the original development demonstrated increased processing speed of about 10%. The basic program structure is illustrated in Fig. 4.
Variable name improvements. Variable names were changed to better conform to the variables they represent. Many variable names were limited to conform to the former FORTRAN limitation of six (6) characters. This is no longer necessary and has been dropped. Variable names were changed to match “standard” nomenclature, such as that used throughout Vallado (2004). Constants were kept as constants in the code, and not assigned as variables with limited precision.

Figure 4. SGP4 Structural Organization. The computer program structure is shown for the original and derivative programs (top) and the revised version for this paper (bottom). Note that the initialization was interspersed throughout the original program, while it is better isolated in the revised code.

(Top flowchart: SGP4 -> INITL -> Near Earth/Deep Space branches with DPPER, DSCOM, DSPACE, SREZ subroutines. Bottom flowchart: SGP4init -> INITL -> DSCOM -> DPPER -> DSINIT; SGP4 -> DSPACE/DPPER branches.)