Constant Chain of Custody

This is the single most critical design constraint in pg_orbit. Get it wrong and positions silently drift by kilometers. There is no runtime check that can detect this class of error after the fact.

The problem

Two-Line Elements are not raw orbital measurements. They are mean elements produced by a differential correction process that fits observed positions against an SGP4 propagator running with a specific set of geopotential constants --- WGS-72. The mean elements absorb geodetic model biases: the eccentricity, inclination, and mean motion encode corrections that only make physical sense when propagated with the same constants used to generate them.

Substituting WGS-84 constants into the propagator does not “upgrade” accuracy. It breaks the internal consistency of the element set. The resulting position error can exceed the natural prediction error of the TLE by an order of magnitude.

This is not theoretical

The WGS-72 equatorial radius is 6378.135 km. The WGS-84 equatorial radius is 6378.137 km. The 2-meter difference looks negligible, but it enters every altitude computation, every semi-major axis derivation, and every GiST index key. Mixing the two constants corrupts results in ways that pass casual inspection but fail against reference implementations.

The rules

Four rules govern constant usage across the entire codebase. No exceptions.

Rule 1: WGS-72 for SGP4/SDP4 propagation

All propagation uses WGS-72 constants: $\mu$, $a_e$, $J_2$, $J_3$, $J_4$, $k_e$. These flow through sat_code’s norad_in.h defines and are never overridden. The functions SGP4_init(), SGP4(), SDP4_init(), and SDP4() operate entirely in the WGS-72 domain.

Rule 2: WGS-84 for coordinate output

Geodetic latitude, longitude, and altitude use the WGS-84 ellipsoid ($a = 6378.137$ km, $f = 1/298.257223563$). This is the modern standard for ground-station positioning, GPS receivers, and mapping services. The conversion from ECEF to geodetic in coord_funcs.c:ecef_to_geodetic() uses WGS-84.

Rule 3: Reduced TEME nutation

The SGP4 output frame (True Equator, Mean Equinox) uses only 4 of the 106 IAU-80 nutation terms. Applying the full nutation model would “correct” for effects that SGP4 already accounts for internally, introducing error rather than removing it.

Rule 4: No other combination is valid

WGS-72 for propagation, WGS-84 for output. Perigee and apogee altitudes use WGS-72 because they derive from mean elements. Geodetic altitude uses WGS-84 because it converts a physical position. There is no scenario where mixing these is correct.

Constant inventory

The complete set of constants, with provenance and location in both pg_orbit and sat_code.

WGS-72 constants (propagation domain)

Source: Hoots & Roehrich, “Models for Propagation of NORAD Element Sets,” Spacetrack Report No. 3, 1980.

Constant	Symbol	Value	types.h	norad_in.h
Gravitational parameter	$\mu$	$398600.8\ \text{km}^3/\text{s}^2$	WGS72_MU	(implicit in $k_e$)
Equatorial radius	$a_e$	$6378.135\ \text{km}$	WGS72_AE	earth_radius_in_km
Zonal harmonic $J_2$	$J_2$	$0.001082616$	WGS72_J2	xj2
Zonal harmonic $J_3$	$J_3$	$-2.53881 \times 10^{-6}$	WGS72_J3	xj3
Zonal harmonic $J_4$	$J_4$	$-1.65597 \times 10^{-6}$	WGS72_J4	xj4
Derived rate constant	$k_e$	$0.0743669161\ \text{min}^{-1}$	WGS72_KE	xke

The rate constant $k_e$ is derived from $\mu$ and $a_e$:

$$k_e = \frac{\sqrt{\mu} \times 60}{a_e^{3/2}}$$

The factor of 60 converts from seconds to minutes, matching the SGP4 convention of radians per minute for mean motion.

WGS-84 constants (output domain)

Source: NIMA TR8350.2, “Department of Defense World Geodetic System 1984.”

Constant	Symbol	Value	types.h
Equatorial radius	$a$	$6378.137\ \text{km}$	WGS84_A
Flattening	$f$	$1/298.257223563$	WGS84_F
Eccentricity squared	$e^2$	$f(2 - f)$	WGS84_E2

Why two copies of AE?

types.h carries a parallel copy of the WGS-72 constants even though sat_code defines them in norad_in.h. This is intentional.

types.h is the single header for all pg_orbit C sources. norad_in.h is an internal sat_code header not meant for external consumers. The GiST index (gist_tle.c) and TLE accessor functions (tle_type.c) need $k_e$ and $a_e$ without pulling in sat_code internals. The values must be identical.

The perigee and apogee altitude computations derive from mean elements:

$$a_{er} = \left(\frac{k_e}{n}\right)^{2/3} \quad \text{[earth radii]}$$ $$\text{perigee}_\text{km} = a_{er} \cdot (1 - e) \cdot a_e - a_e$$ $$\text{apogee}_\text{km} = a_{er} \cdot (1 + e) \cdot a_e - a_e$$

These must use WGS-72 $a_e$ (6378.135 km), not WGS-84 (6378.137 km), because $n$ is a mean motion fitted against the WGS-72 geopotential. Using the wrong radius shifts every altitude by 2 meters. The error compounds in GiST index operations where altitude-band overlap determines whether two orbits are candidates for conjunction screening.

Where the boundary lives

The WGS-72/WGS-84 boundary is crossed in exactly two places in the codebase:

TEME to Geodetic

coord_funcs.c:ecef_to_geodetic() converts ECEF Cartesian coordinates (derived from WGS-72 propagation through GMST rotation) to geodetic latitude, longitude, and altitude on the WGS-84 ellipsoid. This is the correct boundary --- the ECEF position is a physical location in space, and WGS-84 is the standard for expressing that location as geodetic coordinates.

Observer to ECEF

coord_funcs.c:observer_to_ecef() converts a ground station’s geodetic coordinates (on WGS-84, as entered by the user) to ECEF Cartesian for the topocentric transform. The observer’s position is a real-world location defined in WGS-84; converting it to ECEF puts it in the same Cartesian frame as the satellite.

Everything upstream of these functions operates in the WGS-72 domain. Everything downstream operates on physical positions that have already been converted. The boundary is narrow, explicit, and documented in the source comments.

The GMST question

The GMST computation uses the IAU 1982 formula (Vallado Eq. 3-47):

$$\text{GMST} = 67310.54841 + (876600 \times 3600 + 8640184.812866) \cdot T_{UT1} + 0.093104 \cdot T_{UT1}^2 - 6.2 \times 10^{-6} \cdot T_{UT1}^3$$

where $T_{UT1} = (JD - 2451545.0) / 36525.0$, and the result is in seconds of time, converted to radians by multiplying by $\pi / 43200$ and normalized to $[0, 2\pi)$.

pg_orbit deliberately does not use a higher-precision GMST model (e.g., IAU 2000A). The SGP4 output is only accurate to the precision of its own GMST model. Applying a more precise rotation would not improve the final position and could introduce a systematic offset between the propagated TEME position and the Earth-fixed frame.

This is the constant chain of custody in action: match the precision of the input, not the precision available in the literature.

Verification

The chain of custody is verified through the Vallado 518 test vectors --- 518 reference propagations across a range of orbit types (LEO, MEO, GEO, deep-space, high-eccentricity). Every test vector must match to machine epsilon before any other development proceeds.

If a code change causes even one vector to drift, the constant chain has been broken somewhere. The test suite is the enforcement mechanism for the design constraint.

-- Verify against Vallado test vector (ISS-like orbit)
-- Expected: position match within 1e-8 km
SELECT eci_x(sgp4_propagate(tle, epoch + interval '720 minutes'))
FROM vallado_test_vectors
WHERE norad_id = 25544;

For maintainers

If you are modifying any function in tle_type.c, gist_tle.c, or coord_funcs.c, check which constant set you are using. If you find yourself reaching for WGS84_A in a function that computes from mean elements, stop. You are about to break the chain.