The Universal-Variable Lambert Algorithm

Stumpff functions, one equation for every conic, Newton iteration and verification by propagation

Cambridge, Massachusetts, 1960s

At MIT's Instrumentation Laboratory, Richard Battin's group is writing the guidance equations for a computer that will fly to the Moon with 2,048 words of erasable memory. A navigation routine that branches — if ellipse do this, if hyperbola do that, if nearly parabolic do something else — is a routine that eventually takes the wrong branch at the wrong moment. What they want is a formulation in which the type of conic is never asked about. One unknown, one equation, one iteration, valid everywhere. The tool that makes it possible had been sitting in the celestial-mechanics literature for years: the universal anomaly, and a pair of functions named after the German astronomer Karl Stumpff.

Watch first · presented by Shoaib Aziz

The numbers in this film are the labs' own: Curtis Example 5.2 (r₁ = 11,376 km, r₂ = 16,383 km, 1 h apart, Δθ = 100.29°) converges from z = 1.5400 to 1.53986 to 1.5398579 in three Newton steps, giving a = 20,003 km and e = 0.4335; holding the two points, 2 h gives z = 6.40 and 6.15 km/s, 30 min gives z = −1.92 (a hyperbola) and 12.07 km/s, and 15 min gives z = −3.56 and 23.87 km/s; the closed-form C(z) is off by 4 × 10⁻⁸ at z = 10⁻¹⁰ and 4 × 10⁻⁴ at 10⁻¹⁴; and 300 random problems (239 ellipses, 61 hyperbolas) all land within a metre and a millimetre per second.

Recall from the unit on Kepler's equation and orbit propagation that Kepler's equation comes in three versions — elliptic (

M = E - e\sin E

), parabolic (Barker's equation) and hyperbolic (

M = e\sinh F - F

). The lesson on hyperbolic orbits and universal variables replaced all three with the universal anomaly

\chi

; here it is again, now as the unknown of Lambert's problem. Define

\alpha = 1/a

(positive for ellipses, zero for parabolas, negative for hyperbolas) and

z = \alpha\chi^2

. Then every trigonometric or hyperbolic function in the three equations can be written with two series:

C(z) = \frac{1}{2!} - \frac{z}{4!} + \frac{z^2}{6!} - \cdots, \qquad S(z) = \frac{1}{3!} - \frac{z}{5!} + \frac{z^2}{7!} - \cdots

The Stumpff functions. For z > 0, C = (1 − cos√z)/z and S = (√z − sin√z)/√z³; for z < 0 the same with cosh and sinh. At z = 0 they are simply 1/2 and 1/6 — no special case required.

⚠The closed form destroys itself near z = 0

1 - \cos\sqrt{z}

for small

z

subtracts two numbers that agree in almost every digit. At

z = 10^{-10}

the textbook formula for

C

is already wrong in the eighth decimal; at

10^{-14}

it is wrong in the fourth. Near-parabolic transfers live exactly here. Every production implementation switches to the series below some threshold — you will measure the damage yourself in the first lab.

Deriving the Lambert equation in z. Write the propagation

\vec r_2 = f\vec r_1 + g\vec v_1

with the universal-variable Lagrange coefficients, and eliminate everything except

z

. The geometry enters only through a constant

A = \sin\Delta\theta\,\sqrt{\frac{r_1 r_2}{1 - \cos\Delta\theta}}

A carries the transfer angle, and its sign carries the direction: positive for the short way, negative for the long way. A = 0 at Δθ = 180°, which is the same singularity in new clothes.

together with an auxiliary function of

z

y(z) = r_1 + r_2 + A\,\frac{z\,S(z) - 1}{\sqrt{C(z)}}

y must be positive for a real solution — the reason the solver has to start in the right place.

and the time of flight becomes one equation in one unknown:

F(z) = \left[\frac{y(z)}{C(z)}\right]^{3/2} S(z) + A\sqrt{y(z)} - \sqrt{\mu}\,\Delta t = 0

z > 0 ellipse, z = 0 parabola, z < 0 hyperbola. For a single revolution the root lies below z = 4π², where C(z) = 0 and the flight time becomes infinite. F increases monotonically with z, so Newton converges once it starts on the right side.

Loading interactive…

💡The algorithm (Bate–Mueller–White; Curtis Algorithm 5.2)

Compute $r_1$ , $r_2$ and $\Delta\theta$ (with the direction rule), then $A$ .
Start-up. Find a $z$ with $y(z) > 0$ and $F(z) \ge 0$ — for example, march upward from a large negative value. Because $F$ is increasing, this places you on the correct side of the root.
Newton. $z \leftarrow z - F(z)/F'(z)$ , with the analytic derivative

F' = \left(\tfrac{y}{C}\right)^{3/2}\!\left[\tfrac{1}{2z}\!\left(C - \tfrac{3S}{2C}\right) + \tfrac{3S^2}{4C}\right] + \tfrac{A}{8}\!\left[\tfrac{3S}{C}\sqrt{y} + A\sqrt{\tfrac{C}{y}}\right]

(and a separate limit at

z = 0

Velocities. $f = 1 - y/r_1$ , $g = A\sqrt{y/\mu}$ , $\dot g = 1 - y/r_2$ ; $\vec v_1 = (\vec r_2 - f\vec r_1)/g$ , $\vec v_2 = (\dot g\,\vec r_2 - \vec r_1)/g$ .

📐Worked Example: Curtis Example 5.2

\vec r_1 = (5000, 10000, 2100)

km,

\vec r_2 = (-14600, 2500, 7000)

km,

\Delta t = 1

h, prograde.

r_1 = 11{,}376

km,

r_2 = 16{,}383

km.

(\vec r_1\times\vec r_2)_z > 0

, so

\Delta\theta = 100.29°

and

A = 12{,}372

km. Newton from the start-up point:

z = 1.540007

1.5398579

1.5398578966

— converged in three steps. Since

z > 0

the transfer is an ellipse.

y = 13{,}523

km, giving

\vec v_1 = (-5.9925,\ 1.9254,\ 3.2456)

km/s and

\vec v_2 = (-3.3125,\ -4.1966,\ -0.3853)

km/s, with

a = 20{,}003

km,

e = 0.4335

— exactly the answer the Lagrange solver of the previous lesson produced from a completely different equation.

💡Watching z change sign

Hold the endpoints and shorten the flight time. At 2 h,

z = 6.40

(a lazy ellipse,

|\vec v_1| = 6.15

km/s). At 1 h,

z = 1.54

. At 30 min,

z = -1.92

: the transfer is now a hyperbola, with positive energy, and

|\vec v_1| = 12.07

km/s. The algorithm never asked which conic it was solving for. That is the entire point of universal variables.

🛰Trust, then verify

Operational solvers are never trusted on their own word. The standard check is to fly the answer: propagate

(\vec r_1, \vec v_1)

for

\Delta t

with an independent propagator and measure the miss at

\vec r_2

. In the third lab you do this for 300 random problems — and discover that the propagator needs its own care: the textbook starting guess for the universal anomaly diverges on fast hyperbolas, and Vallado's hyperbolic starting value fixes it.

Beware of bugs in the above code; I have only proved it correct, not tried it.

– Donald Knuth, memo to Peter van Emde Boas, 1977

Limitations. The formulation handles single-revolution transfers; multiple revolutions need

z

beyond

4\pi^2

and a second root-finding branch per revolution. The start-up march costs evaluations. And near

\Delta\theta = 180°

A \to 0

and the problem becomes ill-conditioned — no reformulation removes a geometric singularity, but a better one can stop it from hurting everywhere else.

↓But there is a problem

One equation, one unknown, every conic. The universal-variable solver is the workhorse of textbooks and a great deal of real software. But look at what it still needs: a start-up search to find a safe initial

z

, Newton steps that slow down far from the root, and a separate treatment for each extra revolution. In 2015 Dario Izzo at the European Space Agency revisited the problem with a different variable — one in which the time-of-flight curve for every Lambert problem looks almost the same. Same curve, same initial guess, two iterations.