Crawford (1995): Kepler's Equations in C

Why This Matters for SGP4

Paul Crawford is a co-author of the Vallado et al. (2006/2008) Rev-1 paper and the maintainer of the “Dundee code” — the most comprehensive independent SGP4 implementation. This 1995 paper provides the specific Kepler’s equation solver whose bounded Newton-Raphson technique was incorporated into the Rev-1 modernization.

Every SGP4 propagation must solve Kepler’s equation to convert mean anomaly $M$ (which advances linearly with time) into true anomaly $T$ (the satellite’s actual angular position). The original STR#3 FORTRAN used an unbounded Newton-Raphson iteration that could diverge near perigee for high-eccentricity orbits. Crawford’s paper provides the fix.

Rev-1 Section 6 references this paper explicitly for the correction bounding technique.

Kepler’s Equation

The inverse problem — given $M$ , find $E$ and then $T$ — requires solving:

$E = M + e \sin E$

which has no closed-form solution. Newton-Raphson iteration on $f(E) = M + e\sin E - E = 0$ gives:

$E_{N+1} = E_N - \frac{M + e\sin E_N - E_N}{e\cos E_N - 1}$

This converges quadratically near the root (roughly doubling significant figures each iteration), but can diverge wildly near perigee at high eccentricities when the initial estimate is poor.

The Bounded Correction

Crawford’s key insight is Equation 9:

$|E - M| \leq e$

Since $|\sin E| \leq 1$ , the eccentric anomaly can never differ from the mean anomaly by more than the eccentricity. Before each Newton-Raphson step, the correction is clamped:

if (fabs(F) >= fabs(ecc * DF))   // Bound exceeded?
    DE = -SIGN(ecc, F);           // Limit correction to +/- e
else
    DE = F / DF;                   // Normal Newton-Raphson

This simple guard prevents divergence without slowing convergence. For the demanding test case $e = 0.995$ , $M = 0.1$ rad, the iteration converges to $E = 0.842731$ rad in under 10 steps.

Design Principles

Crawford’s implementation satisfies six design objectives that remain relevant for any SGP4 Kepler solver:

Objective	How Achieved
Fast evaluation	Newton-Raphson with precomputed $\sqrt{(1+e)/(1-e)}$
Bounded accuracy	Convergence to `EPS = 1e-8` rad (final error $\approx$ EPS $^2$ )
No singularities	`atan2(sin, cos)` replaces `atan(tan)`
Smooth output	Modulo $2\pi$ reduction with remainder restoration
Derivatives available	Optional $dT/dM$ and $dM/dT$ via chain rule
Portable code	ANSI C, no platform dependencies

The Dundee Code Connection

Crawford’s affiliation (University of Dundee) and his explicit mention of “conversion of some older programmes (generally written in FORTRAN) to ‘C’” is the genesis of what would become the Dundee code. This independent C translation of SGP4 — maintained by Crawford from the mid-1990s through the 2006 modernization — tracked bug fixes and corrections more carefully than any DoD release.

The Vallado (2006) paper identifies three independent SGP4 fix histories: GSFC, Dundee (Crawford), and AFSPC standalone. Crawford’s was the most comprehensive, pre-dating many corrections that the DoD later confirmed.

Convergence Performance

Crawford tested the solver across eccentricities from 0.001 to 0.95, evaluating 4096 mean anomaly points per eccentricity:

Eccentricity	Average Iterations	Worst Case
0.001	2	2
0.01	2	3
0.1	3	4
0.5	4	6
0.9	5	8
0.95	6	9

For typical LEO satellites ( $e < 0.1$ ), Kepler’s equation is solved in 2-3 iterations. Even highly eccentric orbits like Molniya ( $e \approx 0.74$ ) need only 4-5.

Vallado Rev-1 (2006) The SGP4 modernization that incorporated Crawford's Kepler solver

Vallado (2006): Original Paper Documents Crawford's Dundee code and its role in the modernization

Rev-1: Expected Code Updates Section 6 references Crawford (1995) for the correction bounding

STR#3: SGP4 Model The original FORTRAN with unbounded Kepler iteration