The position that has no formula, and is computed anyway

Kepler’s first two laws describe an orbit completely. The path is an ellipse with the primary at a focus, and the radius vector sweeps equal areas in equal times. Between them they contain every fact about where a body goes and how fast.

They do not answer the question anybody actually asks, which is: where will it be on Tuesday.

Getting from a time to a position requires solving one equation, and that equation has no solution in elementary functions. Kepler wrote it down in 1609, said he could not solve it, and added that he did not think anyone would. He was right, and four centuries of orbital mechanics have been conducted by iterating instead.

Three angles for one position #

The trouble begins with a definitional awkwardness: there are three different angles in use for “how far round the orbit the body is”, and none of them can be dispensed with.

The true anomaly $\nu$ is the honest one. It is the angle at the focus, from periapsis to the body, and it is what a telescope would measure. Everything physical is expressed in it. It is also the one that varies most awkwardly with time, racing through the close approach and crawling at the far end, which is the second law in angular form.

The mean anomaly $M$ is the fictitious one. It is defined as $2\pi t/P$ — time, rescaled so that one orbit is $2\pi$. It advances perfectly uniformly by construction and corresponds to no position at all, being the angle a body on a circular orbit of the same period would have reached. It is the input to the problem.

The eccentric anomaly $E$ is the intermediary, and it is the one that makes the problem tractable. It is the angle at the centre of the ellipse — not the focus — to the point on the circumscribing circle directly above the body.

That third one looks arbitrary until the construction is drawn. An ellipse is a circle compressed by the factor $b/a$ in one direction, so every point of the ellipse sits directly below a point of the circle. Taking the angle to the circle’s point rather than to the orbit’s is a change of variable, and it is the change of variable that turns an intractable integral into an equation with two terms.

Where the equation comes from #

The area swept is proportional to the time, which is the second law. Expressing that swept area in terms of $E$ is a short piece of geometry.

The area of the elliptical sector from periapsis to the body is, by the compression, exactly $b/a$ times the area of the corresponding circular sector-plus-triangle. That circular region has area $\tfrac{1}{2}a^2 E$ for the sector at the centre, minus $\tfrac{1}{2}a^2 e\sin E$ for the triangle cut off by the focus’s offset. Multiplying by $b/a$ and setting the result proportional to time gives

$$M = E - e\sin E,$$

which is Kepler’s equation. Two terms, one transcendental function, and no way to invert it.

The direction that works is trivial: give $E$, and $M$ is one subtraction. The direction that is needed is the other one, because time is what is known.

Why there is no formula #

“Transcendental” is a precise claim rather than a complaint, and it is worth stating what it rules out.

An algebraic equation relates its unknown by a finite sequence of additions, multiplications and roots. Kepler’s equation does not: it mixes $E$ with $\sin E$, and no finite combination of the elementary operations produces the inverse. This can be made rigorous — the function $M \mapsto E$ has singularities in the complex plane, at $M = \pm i\,\text{arccosh}(1/e) \mp \sqrt{1-e^2}$, and their positions rule out any closed elementary form.

Those singularities also explain a fact that looks like an accident. Lagrange found a series solution in 1770,

$$E = M + \sum_{n=1}^{\infty} \frac{2}{n}J_n(ne)\sin nM,$$

with $J_n$ the Bessel functions — which are, historically, functions Bessel developed precisely for this problem. The series converges for every $M$ only when $e < 0.6627434\ldots$, and diverges above it. That number is not a rounding of anything; it is the Laplace limit, the eccentricity at which the nearest complex singularity crosses the real axis, and it is the reason nineteenth-century tables of planetary positions stop at the planets and do not handle comets.

The practical significance is small — every planet is far below it — and the theoretical significance is not. A problem that is easy for every case anyone cared about turns out to have a hard boundary sitting just outside, put there by the analytic structure of a function nobody was thinking about.

How it is solved instead #

Newton’s method converges on this equation extremely fast, and the reason is worth knowing because it is not automatic.

Rearranged as a root-finding problem, $f(E) = E - e\sin E - M$, the derivative is $f'(E) = 1 - e\cos E$, which for $e < 1$ is strictly positive. The function is monotonic, there is exactly one root, and the iteration

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

converges quadratically: the number of correct digits doubles at every step.

Four iterations are ample for any orbit a planet is on. That is the whole reason the intractability does not matter: a problem with no closed-form solution and a four-step numerical one is, for every practical purpose, solved.

The remaining difficulty is entirely in the starting guess. Near $e = 1$ the derivative $1 - e\cos E$ approaches zero at periapsis, the iteration’s step becomes enormous, and a bad start can be thrown clear across the orbit or into a divergent cycle. The standard fixes are a starting guess of $E_0 = \pi$ for high eccentricity, which is safely on the correct side of the inflection, and for the very worst cases a bisection to bracket the root before Newton is allowed near it.

What was actually computed #

The equation has a long history of being solved by people who could not afford to solve it many times, and that constraint shaped what the answers looked like.

Before mechanical computation, a planetary position was not calculated on demand. It was looked up. The ephemeris — a table of positions at fixed intervals, published years in advance — existed because the arithmetic was expensive enough that doing it once centrally and printing the result was the economical arrangement. The tables were produced by a small number of people doing an enormous number of iterations by hand, and the interpolation between tabulated entries was itself a studied art. The same tables carried the equation of time, for the same reason and by the same people.

The precision demanded was set by navigation. Finding longitude at sea by lunar distances required the Moon’s position on the sphere to about an arcsecond, and the Nautical Almanac, published continuously since 1767, is a direct descendant of that requirement. Its computers — the word meant a person — were paid by the lunation and their work was independently duplicated, with the two results compared, because a single arithmetic slip propagated into a ship’s position.

The modern situation is the reverse in every respect. A single core solves Kepler’s equation several tens of millions of times a second. A satellite navigation receiver does it continuously, for every satellite in view, because the broadcast ephemeris gives each satellite’s orbital elements and the receiver must convert them into positions before it can do anything else. The equation Kepler could not solve is evaluated, by every phone on the planet, several times a second, forever.

The accuracy required has risen to match. GPS positions are computed to centimetres in the best applications, which means the satellite positions must be known to comparable precision, which means Kepler’s equation must be solved to a relative accuracy of about $10^{-12}$ — comfortably inside what four Newton iterations deliver in double precision, and far outside anything a table could have held.

Running it backwards, several million times #

Everything above treats the elements as known and the time as the input. The commoner situation in practice is the opposite: the observations are known and the orbit is not, and Kepler’s equation then sits in the innermost loop of a fit.

The shape of the problem is a least-squares one. Propose a set of six orbital elements, predict what would have been observed at each of the times something was observed, compare with what was, and adjust. Every one of those predictions requires solving Kepler’s equation once. A fit with a few hundred observations and a few hundred iterations of the outer loop therefore solves it tens of thousands of times, and that is the easy case.

The hard case is exoplanet detection. A radial-velocity curve is fitted by sampling the posterior distribution over the orbital parameters — typically by Markov chain Monte Carlo, with millions of trial parameter sets, each evaluated at every one of the observation times. Solving Kepler’s equation is the single most expensive operation in the whole calculation, and a noticeable fraction of the literature on fast solvers exists because of it rather than because of anything astronomical. Methods that would be pointless for computing one position matter a great deal when the count reaches $10^{9}$: piecewise-polynomial starters that get within $10^{-4}$ before Newton begins, table-driven initial guesses, and formulations that avoid evaluating a sine at all.

There is a pleasing symmetry in that. The equation was intractable in the seventeenth century because arithmetic was scarce, and it is a bottleneck in the twenty-first because the arithmetic is being done a billion times over. What changed is not the difficulty of the equation.

The generalisation: the same equation, three times over #

Kepler’s equation is usually presented as a fact about ellipses. It is one member of a family, and the other members are needed by anything that does not return.

For a parabolic orbit the eccentric anomaly does not exist — there is no circumscribing circle, because there is no far end. Its place is taken by Barker’s equation, a cubic in the parabolic anomaly, which is solvable in closed form by Cardano’s formula. That is a strange inversion: the one case that has an exact solution is the case that never physically occurs, since no real orbit has exactly zero energy.

For a hyperbolic orbit the circular functions become hyperbolic ones,

$$M_h = e\sinh H - H,$$

with $H$ the hyperbolic anomaly. It is transcendental in the same way and is solved the same way, and it is what an interstellar object’s trajectory is propagated with.

Near $e = 1$, both forms become numerically vile — the ellipse’s terms nearly cancel, the hyperbola’s nearly cancel, and double precision loses most of its digits. The standard remedy is a universal variable formulation which handles all three conic types with one equation and one anomaly, with the eccentricity appearing only through a well-behaved parameter. It is what a modern propagator actually uses, and it exists because the near-parabolic case is common: long-period comets live there.

The starting guess, which is where the difficulty actually lives #

Newton’s method converges quadratically once it has caught, and the qualifier is doing more work than it appears to.

The iteration’s step is the residual divided by $1 - e\cos E$, and near periapsis on a high-eccentricity orbit that denominator approaches $1 - e$ — a very small number. A starting guess a little to the wrong side produces an enormous first correction, which can throw the iterate several revolutions away or into a cycle that never settles.

The standard remedies are all about the first guess rather than the iteration. For $e < 0.8$, $E_0 = M$ is safe. Above it, $E_0 = \pi$ is the usual choice: it sits on the correct side of the function’s inflection, so the first step is guaranteed to move toward the root rather than past it. Better starters use a cubic in $M$ fitted to the solution’s shape, reaching $10^{-4}$ before Newton begins and cutting the iteration count by half.

For the hardest cases the fix is to give up on pure Newton and bracket the root first by bisection, which cannot diverge, then switch. That combination is what a production solver uses, and its existence is the honest measure of how badly the pure method can behave.

The direction that is not a problem #

It is worth noticing that only one of the two directions is hard, because the asymmetry says what the difficulty actually is.

Given a time, finding the position requires solving $M = E - e\sin E$ for $E$, which is the whole subject of this essay. Given a position, finding the time requires evaluating $E - e\sin E$, which is one line of arithmetic with no iteration in it at all.

So the equation is not hard in itself; it is hard to invert. And the reason is that it mixes an angle with the sine of the same angle — a transcendental combination that cannot be untangled by any finite sequence of the operations that elementary functions are built from. Kepler’s own remark, that the arc and the sine are of different natures, is exactly this observation two centuries before the language for it existed.

That asymmetry decided the form of every printed ephemeris. A table of positions against time was expensive to produce, because each entry cost an iteration by hand, and cheap to use. A navigator or an observer needed the answer in the hard direction and was given a table computed once, by somebody else, at the cost of many thousands of iterations — which is why national almanac offices existed, employing human computers, from the eighteenth century until the middle of the twentieth.

A published table is an inversion performed in advance and sold, and the whole institution grew out of one equation’s refusal to be rearranged.

One more setting sits at the eccentricity where the iteration is hardest.

Where the model stops #

Two bodies. The equation is exact for a two-body orbit and for nothing else. In a real system the elements themselves drift, and propagating a position over long spans means re-solving with elements that are functions of time.

Constant elements. Even between perturbations, atmospheric drag and radiation pressure change the orbit continuously. A low satellite’s mean anomaly is the least predictable of its elements for exactly this reason, and it is why re-entry dates are quoted with enormous error bars.

Newtonian. The relativistic two-body problem does not close — the orbit precesses — so there is no periodic $M$ to invert.

Numerical, not analytic. Every position produced by this equation carries the error of the iteration and of the arithmetic, which is small and is not zero. Over the millions of revolutions involved in a long-term integration of the solar system, that residue accumulates, and controlling it is most of the work of such an integration.

The figures have a limitation worth stating explicitly, because it is the same one that makes the three anomalies hard to hold in mind. Each figure shows one instant: a single value of $E$, with the three angles marked. What the equation is actually about is how those three angles run at different rates as time passes, and a static picture cannot show a rate. Reading the three figures as a sequence in eccentricity is not the same as watching one orbit go round, and the second is the thing the equation describes.

The construction is the same at every eccentricity and the difficulty is not, so it is worth drawing at two more values that bracket the ones above.

The ladder from here #

Later rungs on this anchor: the series solution and the Laplace limit derived. Starter formulae for high eccentricity, and the ones that fail. Universal variables and the Stumpff functions. Barker’s equation and the parabolic case. The hyperbolic anomaly. Gauss’s method of orbit determination, which runs the whole machinery backwards from three observations. Equinoctial elements, which remove the singularities at zero eccentricity and zero inclination. And the numerical propagators that abandon elements altogether and integrate the equations of motion directly, which is what anything with a thruster requires.

Kepler wrote that he was “sufficiently satisfied that it cannot be solved a priori, on account of the different nature of the arc and the sine”. That is an accurate statement of the transcendence of the equation, made two hundred years before the concept existed.