Kepler's Equation

Mean, eccentric and true anomaly, and the iterative solution of a transcendental equation

Prague, 1605

Johannes Kepler has spent the better part of a decade on the orbit of Mars, working from Tycho Brahe's observations — the most precise measurements ever made without a telescope, and Kepler had to fight the Brahe family in court to keep them. He has the answer. The orbit is an ellipse with the Sun at one focus. He has the second law too: the line from Sun to planet sweeps equal areas in equal times. And he is stuck. Because an astronomer does not want to know the shape of Mars's orbit. An astronomer wants to know where Mars will be on the fourteenth of next month, and the ellipse alone will not say. Kepler writes down the relation that ought to answer it — a single line, one unknown, no ambiguity anywhere. And he cannot solve it. Not because he lacks the algebra. Because there is no solution to find. He says so himself, in the Astronomia Nova: he asks the reader's indulgence and offers an approximation instead. Four hundred and twenty years later the situation is unchanged. Every spacecraft in the solar system carries software that solves Kepler's equation by iteration, thousands of times a second, because after four centuries and some of the best mathematicians who ever lived, nobody has found anything better.

Watch first · presented by Shoaib Aziz

The numbers in this film are the lesson's and the labs' own: for a = 10,000 km and e = 0.3 the period is 9,952 s, and 45 minutes after perigee M = 97.7°, E = 113.4° and ν = 128.6°, with the satellite 11,193 km out — 31° (about 6,000 km) from where steady motion would put it. In the labs: E = 78.71° for M = 45° at e = 0.6; 4 Newton passes at e = 0.6 and 8 at e = 0.98 near perigee; perigee to 90° in 1,552 s (15.6% of the period); and an equation-of-centre peak of 23.0° at e = 0.2.

The trouble is visible in something you already proved. Unit 1 established that the areal velocity is constant:

\frac{dA}{dt} = \frac{h}{2}

Equal areas in equal times. But the position on the orbit is an angle, and on an ellipse a fixed area does not correspond to a fixed angle. Near perigee the radius is short, so sweeping a given area requires swinging through a large angle quickly. Near apogee the radius is long, and the same area is covered by a small, lazy angle.

\nu

is not linear in time. It is not linear in anything convenient. And the more eccentric the orbit, the worse the distortion gets.

💡Deriving the orbital period

Before the hard part, collect a debt from Unit 1. We promised to derive Kepler's Third Law properly, and constant areal velocity delivers it in two lines. Sweeping the entire ellipse takes one period, and the area of an ellipse is

\pi a b

T = \frac{\text{total area}}{\text{areal rate}} = \frac{\pi a b}{h/2} = \frac{2\pi a b}{h}

Now substitute

b = a\sqrt{1-e^2}

and

h = \sqrt{\mu p} = \sqrt{\mu a(1-e^2)}

T = \frac{2\pi a \cdot a\sqrt{1-e^2}}{\sqrt{\mu a (1-e^2)}} = \frac{2\pi a^{3/2}}{\sqrt{\mu}}

The eccentricity cancels completely. The period depends only on $a$ — a fact Kepler discovered empirically in 1619 by staring at planetary data, and which falls out here as a consequence of the two-body equation.

T = 2\pi\sqrt{\frac{a^3}{\mu}}, \qquad n = \frac{2\pi}{T} = \sqrt{\frac{\mu}{a^3}}

Kepler's Third Law and the mean motion n — the average angular rate, in radians per second. Kepler's Third Law is now a theorem, not an observation.

Now the hard part. We want

\nu(t)

, and we cannot get it directly. So we go around, through two invented angles that are easier to work with.

The mean anomaly. Imagine a fictitious satellite that goes round a circle of radius

a

at the constant rate

n

, starting from perigee at the same moment as the real one and completing each lap at the same instant. Its angle is the mean anomaly:

M = n(t - t_{p}) = \sqrt{\frac{\mu}{a^3}}\,(t - t_{p})

The mean anomaly is not an angle you can point at — no physical object is at that angle. It is time, wearing the costume of an angle. $t_{p}$ is the epoch of perigee passage.

M

is trivially easy: it is linear in time by construction. It is also physically meaningless on its own. The whole problem is the bridge from

M

back to the real

\nu

The eccentric anomaly. The bridge runs through a geometric trick that goes back to Apollonius. Draw the auxiliary circle — the circle of radius

a

that circumscribes the ellipse. From the satellite's true position, drop a perpendicular to the major axis and continue it until it meets the circle. The angle from the centre of the ellipse to that point is the eccentric anomaly

E

Figure 5.1

The auxiliary-circle construction. Three angles, three different vertices: ν is measured at the focus, E at the centre, and M is not measured anywhere at all — it belongs to a fictitious body on the circle. At perigee and apogee all three coincide exactly; in between they diverge, and the more eccentric the orbit, the wider the spread.

The eccentric anomaly earns its place because it makes both of the relations we need clean. Position on the orbit is simply

r = a(1 - e\cos E), \qquad \tan\frac{\nu}{2} = \sqrt{\frac{1+e}{1-e}}\,\tan\frac{E}{2}

Radius and true anomaly, both from E, both exact and both closed-form. Getting to E is the only difficulty in the entire problem.

And now, connecting

E

back to time. Work out the swept area in terms of

E

, set it equal to

(h/2)(t - t_p)

, and after a page of geometry the entire time problem collapses to one line:

\boxed{\;M = E - e\sin E\;}

Kepler's equation. One unknown, no ambiguity, and no solution in elementary functions. This is the wall.

⚠Why Kepler's equation has no closed-form solution

It is tempting to assume that

E

appearing both alone and inside a sine is merely inconvenient. It is worse than that. The equation is transcendental: it mixes an algebraic term with a transcendental one, and there is no finite combination of elementary functions that inverts it. This is not a gap in our cleverness. In 1824 Bessel constructed an exact infinite-series solution using what are now called Bessel functions — invented, in part, for this problem. Series, yes. Closed form, no. The list of people who attacked it is not a list of second-raters: Newton, Lagrange, Gauss, Bessel, Laplace, Cauchy. Between them they produced dozens of methods, all of which are iterative or infinite. So every orbit propagator ever written, in every language, on every spacecraft, contains a loop. There is no way around it, and there never will be.

The practical answer is Newton–Raphson, and for this particular equation it is superb. Write

f(E) = E - e\sin E - M

and iterate:

E_{k+1} = E_{k} - \frac{E_k - e\sin E_{k} - M}{1 - e\cos E_k}

Kepler by Newton–Raphson. Starting from E₀ = M, a typical satellite converges to machine precision in three or four passes.

Figure 5.2

Newton–Raphson closing on the root. Each step slides down the tangent line to the axis, and the error is squared every iteration: 10⁻¹, then 10⁻², then 10⁻⁴, then 10⁻⁸, then 10⁻¹⁶. Quadratic convergence is why nobody minds that the equation is unsolvable.

⚠Convergence near e → 1

Look at the denominator:

1 - e\cos E

. As

e \to 1

and

E \to 0

— that is, a nearly-parabolic orbit near perigee — it approaches zero, and Newton's method starts dividing by almost nothing. Steps become enormous, iterates jump wildly, and convergence stalls or fails. This is not exotic. Comets have

e = 0.999

and worse, and perigee is exactly where you most want accuracy. The standard repairs, in escalating order: Better starting guess.

E_0 = M + e\sin M

instead of

E_0 = M

, which absorbs most of the first correction. For high

e

, use

E_0 = \pi

. Damping. Cap the step size so a wild derivative cannot throw the iterate across the orbit. Higher order. Halley's or Laguerre's method, which use the second derivative and are far more robust for

e > 0.95

. Change the equation. For near-parabolic trajectories, switch to Barker's equation (

e = 1

exactly, and it does have a closed-form solution) or to universal variables, which handle every conic with a single formulation. The pattern is the same one from Unit 3: the physics is fine, and the parameterisation is what is failing.

The gap between where the satellite really is and where the fictitious mean satellite thinks it is has a name: the equation of centre,

\nu - M

. For small

e

it has a tidy series:

\nu - M \approx 2e\sin M + \frac{5}{4}e^2\sin 2M + O(e^3)

The equation of centre. To first order the peak error is 2e radians — for e = 0.2 that is about 23°, which is a very long way to be wrong about where a satellite is.

Figure 5.3

M, E and ν plotted against time for one revolution. M is a straight line by definition. E bows away from it, and ν bows further still. All three meet exactly at perigee and apogee and diverge most in between — which is precisely where the ground station needs to be pointing.

📐Worked Example: where is it, 45 minutes after perigee?

A satellite has

a = 10{,}000

km and

e = 0.3

. Find its position 45 minutes after perigee passage. Mean motion.

n = \sqrt{398600/10000^3} = \sqrt{3.986\times10^{-7}} = 6.3135\times10^{-4}

rad/s. Period.

T = 2\pi/n = 9952

s = 165.9 min. So 45 minutes is 27% of the way round. Mean anomaly.

M = nt = 6.3135\times10^{-4} \times 2700 = 1.7046

rad

= 97.67°

. Solve Kepler. Start at

E_0 = M = 1.7046

E_1 = 1.7046 - \dfrac{1.7046 - 0.3\sin(1.7046) - 1.7046}{1 - 0.3\cos(1.7046)} = 1.9905

E_2 = 1.97990

E_3 = 1.97989

E_4 = 1.97989

— converged, four iterations,

E = 113.44°

. True anomaly.

\tan(\nu/2) = \sqrt{1.3/0.7}\tan(56.72°) = 1.3628 \times 1.5235 = 2.0762

, so

\nu/2 = 64.28°

and

\nu = \mathbf{128.6°}

. Radius.

r = a(1 - e\cos E) = 10000(1 - 0.3\cos 113.44°) = 10000(1 + 0.1193) = \mathbf{11{,}193}

km. Read the numbers.

M = 97.7°

E = 113.4°

\nu = 128.6°

. The fictitious mean satellite thinks it has gone 98° round. The real one has gone 129°. Assume uniform motion and you would be 31° wrong — with the satellite 11,193 km out, that is a position error of roughly six thousand kilometres.

🛰Kepler's equation inside a GPS receiver

Every GPS satellite broadcasts a navigation message containing its own orbital elements. Not its position — its elements, including the mean anomaly at a reference epoch and the mean motion. Which means the receiver in your pocket is solving Kepler's equation. Continuously. For every satellite it can see, several times a second, before it can even begin the trilateration that produces a fix. The published GPS interface specification even mandates the algorithm — iterate until convergence — because a receiver that got this wrong would produce position errors of kilometres. Kepler's unsolvable equation is being solved, right now, in a few billion pockets simultaneously.

↓But there is a problem

So the wall has a door, and the door is a loop. Give me an orbit and a time, and I can now tell you exactly where the satellite is: compute

M

, iterate to

E

, convert to

\nu

and

r

. Four hundred years of mathematics, reduced to about eight lines of code. Run it in a loop, stepping time forward, and you can watch the future unfold. Ground tracks that actually track. Passes over a station. A conjunction warning three days before the conjunction. But look again at the equation.

M = E - e\sin E

is about an ellipse: it needs a circle drawn around the orbit, and a period to wrap around. A spacecraft leaving for Mars has neither. Nor does a comet falling in at almost exactly escape speed, where

e

sits a hair from 1 and the formulas of this lesson start dividing by numbers close to zero. And a tracking station never hands you an orbit anyway. It hands you a position and a velocity, and the shape of the orbit is something you find out afterwards. What is needed is one equation that works for every orbit, without ever asking which kind it is.