Solve Kepler's equation M = E − e·sin(E) for the eccentric anomaly E,
given the mean anomaly M (radians) and eccentricity e (0 ≤ e < 1) — the
one step in orbital position math that has no closed form. Returns the E
for the M actually given (not just its 0–2π remainder), so
E − e·sin(E) = M holds directly.
Real, observed bug history: this used to be plain Newton-Raphson seeded
with E = M, documented as converging "for any bound orbit". It doesn't:
near periapsis of a very eccentric orbit the tangent is almost flat, the
first step overshoots by whole revolutions, and 15 iterations weren't
enough to come back — 1.3% of mean anomalies failed at e = 0.98, 5.3% at
0.99 (e = 0.99, M = 5.9803 returned E ≈ −3430). Because E only ever goes
through sin/cos, the result was still a point ON the ellipse, just the
wrong one: no NaN, nothing to notice. Halley (e ≈ 0.967) only just
escaped it.
Two standard fixes, together:
seed with E = π above e = 0.8, from where Newton converges
monotonically for any M (and with M + e·sin M below it, a better start
than M itself for low e);
safeguard with bisection: f(E) = E − e·sin E − M is strictly
increasing for e < 1 and brackets its root in [0, 2π] for M in that
range, so any Newton step that would leave the current bracket is
replaced by a bisection step. That can't diverge, whatever e and M.
Solve Kepler's equation M = E − e·sin(E) for the eccentric anomaly E, given the mean anomaly M (radians) and eccentricity e (0 ≤ e < 1) — the one step in orbital position math that has no closed form. Returns the E for the M actually given (not just its 0–2π remainder), so E − e·sin(E) = M holds directly.
Real, observed bug history: this used to be plain Newton-Raphson seeded with E = M, documented as converging "for any bound orbit". It doesn't: near periapsis of a very eccentric orbit the tangent is almost flat, the first step overshoots by whole revolutions, and 15 iterations weren't enough to come back — 1.3% of mean anomalies failed at e = 0.98, 5.3% at 0.99 (e = 0.99, M = 5.9803 returned E ≈ −3430). Because E only ever goes through sin/cos, the result was still a point ON the ellipse, just the wrong one: no NaN, nothing to notice. Halley (e ≈ 0.967) only just escaped it.
Two standard fixes, together: