Orbits

The one solve that does not ask which conic it is

Kepler's equation is for ellipses, Barker's cubic for parabolas, and a hyperbolic sine for the rest — three parameterisations of one motion, each worst exactly where its neighbour takes over. The universal variable removes the question, and the removal is not a convenience.

Assumes Angular momentum, Conic sections and Vis-viva.

An inverse-square force allows four curves and the eccentricity alone decides which. That is a statement about shape, and it is clean. The statement about time is not clean at all. Position against time on an ellipse comes from solving M=EesinEM = E - e\sin E; on a parabola from a cubic in tan(ν/2)\tan(\nu/2); on a hyperbola from M=esinhHHM = e\sinh H - H. Three equations, three variables that mean three different things, and a program that computes an orbit has to ask which one it is holding before it can say where the body will be.

The awkward part is not the asking. It is that each of the three is at its worst exactly where the next one begins.

One analytic function each, and a case split that is only in the writing. Left: the Stumpff functions c₂ and c₃ across z = αχ², from z = −30 to 30. Positive z is an ellipse, negative z a hyperbola, and z = 0 exactly a parabola — and there is nothing at the origin: c₂ passes through 0.5 and c₃ through 0.166667 with every derivative continuous, because each is the single power series Σ(−z)^j/(k+2j)!. Written as (1 − cos√z)/z and (√z − sin√z)/z^3/2 they need three cases; written as their series they need none, and the case split in the code above is about floating point rather than about the mathematics. Right: how much floating point. The closed form subtracts two numbers that agree ever more closely as z → 0, so it loses digits — 11.9 of the sixteen it started with at z = 10⁻¹², measured against the series rather than estimated. That is why the drawn curve switches to the series inside |z| < 0.001, and why an implementation that uses the closed form everywhere is worst exactly at the parabola it was written to include.
Fig. 1 The functions the whole method rests on, and what the usual way of writing them costs. Left: c2c_2 and c3c_3 across z=αχ2z = \alpha\chi^2, where positive zz is an ellipse, negative a hyperbola and zero exactly a parabola. There is nothing at the origin — c2c_2 passes through 1/21/2 and c3c_3 through 1/61/6 with every derivative continuous — because each is one power series, (z)j/(k+2j)!\sum(-z)^j/(k+2j)!, with no branch in it anywhere. Written instead as (1cosz)/z(1-\cos\sqrt z)/z and (zsinz)/z3/2(\sqrt z - \sin\sqrt z)/z^{3/2} they need three cases, and the cases are about arithmetic rather than about mathematics. Right: how much arithmetic. The closed form subtracts two numbers that agree ever more closely as z0z\to 0, losing 11.9 of its sixteen digits at z=1012z = 10^{-12} — measured against the series rather than estimated.

Three equations that are one equation

Kepler’s equation is a statement about area. The eccentric anomaly EE is an angle on the circle that circumscribes the ellipse, the mean anomaly MM is proportional to the area swept, and M=EesinEM = E - e\sin E is the equal-area law written as a relation between two angles. It is transcendental, has no closed form, and is solved anyway — in four Newton steps for a planet, in seven for a comet. Push ee towards 1 and the construction loses its meaning by degrees. The circumscribing circle grows without bound; the ellipse it circumscribes becomes a long thin sliver; and EE, the angle on that circle, comes to be a quantity of which almost all values describe the far end of an orbit nobody is looking at. At e=0.9999e = 0.9999 the whole of the interesting motion — the part near periapsis, where the comet is bright and the observations are — occupies a fraction of a degree of EE. The equation still holds. The parameterisation has become a bad set of coordinates for the question.

At e=1e = 1 exactly there is a different equation and it is better than either neighbour: Barker’s,

μ2q3t=D+D33,D=tanν2,\sqrt{\frac{\mu}{2q^3}}\,t = D + \frac{D^3}{3}, \qquad D = \tan\frac{\nu}{2},

a depressed cubic with one real root, solvable in closed form by Cardano. No iteration, no transcendence, no difficulty of any kind. And it applies to precisely one value of the eccentricity, which no real orbit has.

Above 1 the hyperbolic form takes over, M=esinhHHM = e\sinh H - H, and it inherits its own version of the trouble from below. Just above e=1e = 1 the two terms are nearly equal and their difference is the answer; subtracting them throws away digits in proportion to how nearly parabolic the orbit is.

36 steps round one orbit, spaced two ways. The same orbit at e = 0.9, marked at 36 equal steps of time on the left and 36 equal steps of the fictitious time χ on the right. Equal time steps are what a fixed-step integrator takes, and they put the marks where the body is slow: the longest gap is 11.8 times the shortest, and the pericentre passage — the only part of the orbit where anything interesting happens to the trajectory — is crossed in 1 of them. Stepping in χ instead makes the time step proportional to r, and the marks come out 2.25-to-one even, with 3 inside the same pericentre region. Nobody chose that. dt/dχ = r/√μ is the substitution itself, and since vr equals the angular momentum at both apses, a fixed step in χ is very nearly a fixed step in arc length — an adaptive integrator's behaviour obtained from a change of variable rather than from a controller.
Fig. 2 And what it buys an integrator. At an eccentricity of 0.9 a fixed step in time is far too coarse at periapsis and wastefully fine at apoapsis; a fixed step in the regularised variable is uniform in neither and appropriate in both, because the new independent variable advances with the true anomaly rather than with the clock. That is the practical reason the universal formulation is used, and it has nothing to do with hyperbolic orbits: it is about where the steps go.

The substitution that has no branch in it

The universal formulation begins by refusing to use aa and using α=1/a\alpha = 1/a instead.

That is the whole trick, and it is worth pausing on because it looks like nothing. The semi-major axis is finite and positive on an ellipse, infinite on a parabola, and finite and negative on a hyperbola. Any formula written in terms of aa therefore has a singularity in the middle of the family, and every branch in the classical treatment is that singularity being stepped around. The reciprocal has none: α\alpha runs smoothly from positive through zero to negative as the eccentricity passes 1, and zero is not a special value of it in any way. Since α=2ε/μ\alpha = -2\varepsilon/\mu with ε\varepsilon the specific orbital energy, what has really happened is that the parameter has been changed from a length to an energy — and energy passes through zero without noticing. The independent variable changes too. In place of an angle, a quantity χ\chi with dimensions of the square root of a length, defined so that dχ=μdt/rd\chi = \sqrt{\mu}\,dt/r. Then the time of flight from a state (r0,vr0)(r_0, v_{r0}) is

μΔt=σ0χ2c2(z)+(1αr0)χ3c3(z)+r0χ,z=αχ2,\sqrt{\mu}\,\Delta t = \sigma_0 \chi^2 c_2(z) + (1 - \alpha r_0)\chi^3 c_3(z) + r_0\chi, \qquad z = \alpha\chi^2,

with σ0=r0vr0/μ\sigma_0 = r_0 v_{r0}/\sqrt\mu, and the radius along the way is

r=χ2c2(z)+σ0χ(1zc3)+r0(1zc2).r = \chi^2 c_2(z) + \sigma_0\chi\,(1 - z c_3) + r_0(1 - z c_2).

The functions c2c_2 and c3c_3 are the Stumpff functions. On an ellipse χ=aΔE\chi = \sqrt a\,\Delta E and the first equation reduces to Kepler’s; on a hyperbola χ=aΔH\chi = \sqrt{|a|}\,\Delta H and it reduces to the hyperbolic form; at α=0\alpha = 0 it reduces to Barker’s cubic, since c2=1/2c_2 = 1/2 and c3=1/6c_3 = 1/6 and the expression becomes a polynomial in χ\chi of degree three. The three classical equations are not analogous. They are the same equation at three values of one parameter.

One squaring, and the singularity is gone. Left, a harmonic oscillator: an ellipse centred on the origin, marked at 24 equal steps of its own phase. Right, the same points after squaring as complex numbers, u ↦ u². The image is an ellipse with the origin at a focus — checked here by the focal property, r₁ + r₂ = 2a to nine figures at four points — with semi-axes (A²+B²)/2 = 2.5000 and AB = 1.5000 and the focus at (A²−B²)/2 = 2.0000, which is ae exactly. Three things follow at once. The angle doubles, so one turn of the oscillator is one whole orbit traversed twice as fast in phase; the equal phase steps on the left arrive as the eccentric anomaly on the right, which is why that anomaly and not the true one is what the equations want; and the collision at r = 0, where the inverse square is infinite, is the point u = 0, where the oscillator has a perfectly ordinary velocity. The singularity was a property of the coordinates.
Fig. 3 What the change of variable does to the singularity, rather than to the case split. Regularising the two-body problem in the Levi-Civita variables turns the collision singularity into a regular point of a harmonic oscillator: the equations become linear, the step size that was collapsing near periapsis becomes uniform, and an integrator that was losing accuracy at closest approach stops noticing it. The universal formulation and the regularisation are the same idea applied to two different problems — one removes a branch on the conic class, the other removes a divergence in the time.
One Newton loop, four conics, 6–11 steps. The residual |√μΔt − f(χ)| after each Newton step of the universal Kepler equation, at e = 0.6, e = 1, e = 1.0001, e = 1.4, all through the same code, the same starting rule and the same convergence test. Each orbit has periapsis at 0.5 AU and the solve is for 315.0 days past it. The step counts are 6 at e = 0.6, 10 at e = 1, 11 at e = 1.0001, 9 at e = 1.4, and the pair that matters is the crossing: 10 steps at e = 1 against 11 at e = 1.0001, a difference of 1. The spread across the whole set is 1.83, and it is a property of the STARTING GUESS rather than of the conic class — the run that begins nearest its own answer finishes first, on every branch. What the figure is about is the absence of a step at e = 1: the derivative of the equation being solved is r itself, which is positive on every conic, so the iteration has nothing to notice there. The classical forms have the opposite property — Kepler's equation degenerates as e → 1 from below, Barker's cubic is exact only AT 1, and e sinh H − H loses digits to cancellation just above it.
Fig. 4 The solve itself, across the classes it does not distinguish between. One Newton iteration on one function of one variable converges for the ellipse, the parabola, the orbit a ten-thousandth past parabolic, and the hyperbola — and the iteration counts are within one of each other. There is no branch anywhere in the code: the conic class enters only as the sign of a number the solver never tests, and the near-parabolic case that defeats every classical formulation is not special here at all.

Why the derivative is the thing to look at

Newton’s method on any of these forms needs the derivative of the left-hand side, and in the universal form that derivative is

d(μΔt)dχ=r.\frac{d(\sqrt\mu\,\Delta t)}{d\chi} = r.

The distance from the primary. Not an approximation to it, not a quantity resembling it: the radius itself, which the formulation is computing anyway. This is why the iteration does not care which conic it is on. A Newton step divides by the derivative, so a solve becomes ill-conditioned where its derivative approaches zero — and rr is strictly positive on every orbit of every class. It is small at periapsis and large far out, so the step size varies, but the sign never changes and the value never vanishes. Compare Kepler’s equation, whose derivative is 1ecosE1 - e\cos E: at e1e \to 1 and E0E \to 0 that goes to zero, and the near-parabolic ellipse near periapsis is precisely the case where a Newton step wants to divide by almost nothing.

The step counts in the figure above make the point without argument. Six steps at e=0.6e = 0.6, ten at the parabola, eleven a ten-thousandth above it, nine well beyond. Whatever varies there is the quality of the first guess. Nothing varies at e=1e = 1.

Four conics through one periapsis, drawn by one solve. Distance from the Sun against time for four orbits sharing a periapsis of 0.5 AU, at e = 0.3, e = 1, e = 1.0001, e = 2, over 900 days. Every point on every curve came from the same universal Kepler solve — no branch on the conic class anywhere in it — and each curve was then checked against the classical solution of its own kind at the midpoint: Kepler's equation at e = 0.3 agrees to machine precision; Barker's cubic at e = 1 agrees to machine precision; e sinh H − H at e = 1.0001 agrees to machine precision; e sinh H − H at e = 2 agrees to machine precision. The curve to read twice is e = 1.0001: over this arc it lies within 0.02% of the parabola and is indistinguishable from it, and it is the only one of the four whose fate the drawing cannot show. What the classical parameterisation costs there is arithmetic rather than impossibility: at the midpoint of this arc, e sinh H − H throws away 3.3 of its sixteen digits to cancellation, against 0.1 at e = 2 — enough to matter to an ephemeris and not enough to stop a plot.
Fig. 5 The same four conics with the eccentricities pushed further apart — 0.3, exactly parabolic, a hair above, and a strongly hyperbolic 2.0. All four are drawn by one function of one variable, and the near-parabolic pair still lie on top of each other for most of the sweep. Nothing in the code branches on which conic it is, and that is the whole claim: the universal formulation has no case analysis in it, so the pair that differ in the fifth decimal place cost it nothing at all.

What the series is for

Between the mathematics and the arithmetic there is one more layer, and it is the only place a case split survives.

Written in closed form, c2(z)=(1cosz)/zc_2(z) = (1 - \cos\sqrt z)/z. As z0z \to 0 the numerator is a difference of two numbers that agree to more and more digits, and floating-point subtraction of nearly equal numbers destroys precision in exact proportion to how nearly equal they are. At z=1012z = 10^{-12}, cosz\cos\sqrt z differs from 1 in the twelfth decimal place; subtracting gives four significant digits out of sixteen, and dividing by zz scales the error up rather than away.

The series has no such problem. c2(z)=12z24+z2720c_2(z) = \tfrac12 - \tfrac z{24} + \tfrac{z^2}{720} - \cdots is a sum of terms of decreasing size with no cancellation between them, and it converges fastest exactly where the closed form is worst.

So an implementation uses the closed form away from the origin and the series near it, and the switch is at a threshold chosen by measuring where the two agree to full precision. It is worth being exact about what that case split is and is not. It is not the branch on the conic class returning by another door: it does not test the eccentricity, it does not distinguish an ellipse from a hyperbola, and it does not change the function being computed. It tests whether a subtraction is about to lose digits. The mathematical object has no branch; the machine that evaluates it has a floating-point format.

What was actually computed, and what was measured

Nothing on this page is an observation. Every number here comes out of arithmetic performed on two given quantities, and the honest description of the whole essay is that it is about a method rather than about the sky.

The observations enter one step earlier and they are worth naming, because the method exists to serve them. An orbit is fitted to a sequence of astrometric positions — two angles at a time, with no distance in any of them — and the fit works by predicting where the body should have been at each observed epoch and adjusting the elements until the predictions land on the observations. That prediction is a solve of exactly the kind above, performed once per observation per iteration of the fit, thousands of times for a well-observed object. A solve that needs a different code path either side of e=1e = 1 is a fit whose derivatives with respect to the elements are discontinuous at e=1e = 1, which is where a marginally-unbound comet’s fit spends its time.

The measurement that shows this is real, rather than tidy, is the population of comets catalogued with eccentricities between 0.999 and 1.001. Their elements were determined by iterating the very quantity that all three classical forms handle badly. The reason a distinguished body of numerical literature exists about near-parabolic motion is not aesthetic.

Four conics through one periapsis, drawn by one solve. Distance from the Sun against time for four orbits sharing a periapsis of 1.2 AU, at e = 0.6, e = 1, e = 1.0001, e = 1.4, over 900 days. Every point on every curve came from the same universal Kepler solve — no branch on the conic class anywhere in it — and each curve was then checked against the classical solution of its own kind at the midpoint: Kepler's equation at e = 0.6 agrees to machine precision; Barker's cubic at e = 1 agrees to machine precision; e sinh H − H at e = 1.0001 agrees to machine precision; e sinh H − H at e = 1.4 agrees to machine precision. The curve to read twice is e = 1.0001: over this arc it lies within 0.01% of the parabola and is indistinguishable from it, and it is the only one of the four whose fate the drawing cannot show. What the classical parameterisation costs there is arithmetic rather than impossibility: at the midpoint of this arc, e sinh H − H throws away 3.7 of its sixteen digits to cancellation, against 0.3 at e = 1.4 — enough to matter to an ephemeris and not enough to stop a plot.
Fig. 6 The same four conics through a periapsis at 1.2 astronomical units rather than 0.5. Everything scales and nothing about the construction changes — which is worth seeing, because the periapsis distance is the quantity the classical formulations divide by. The universal variable’s insensitivity to qq is the same insensitivity that makes it work near a parabola, and the two are one property rather than two.

The case this repository could not draw

There is a specific and slightly embarrassing piece of evidence for all of this, and it is on this site.

The generator that draws hyperbolic orbits here computes distance against time from M=esinhHHM = e\sinh H - H, in the classical way. It carries an assertion of its own — that the drawn curve’s secants agree with the radial energy equation — and at e=1.0001e = 1.0001 that assertion fails: the residual comes out at several per cent of the local radial speed against a two per cent tolerance. The failure was chased, two unrelated defects in the tolerance were found and fixed along the way, and the case still would not draw. It was recorded in the fleet’s consolidation notes as an open item with the cause unknown, and the note ends by guessing that the parameterisation wants the universal-variable form.

That guess was right, and the figure above is what it looks like when it is taken up. The e=1.0001e = 1.0001 curve there is drawn by the universal solve and checked against the classical answer at its midpoint; the two agree, and the classical form’s cancellation is measured rather than asserted — 3.3 digits of sixteen at that point, against 0.2 digits for the comfortable hyperbola beside it. Three lost digits do not stop a plot. They do stop an ephemeris that is differencing positions to get a velocity, which is what an orbit determination does at every step.

The surprising part is which direction the difficulty runs. A reader meeting the three classical forms for the first time would guess that the hard case is the extreme one — the wildly hyperbolic orbit, the nearly circular one — and that the ordinary middle of the range is safe. It is the opposite. The boundary between two well-behaved regimes is where both are worst, and it is worst there because it is a boundary: each form’s variable is defined by a construction that degenerates as the other’s takes over.

Four conics through one periapsis, drawn by one solve. Distance from the Sun against time for four orbits sharing a periapsis of 0.5 AU, at e = 0.6, e = 1, e = 1.0001, e = 1.4, over 900 days. Every point on every curve came from the same universal Kepler solve — no branch on the conic class anywhere in it — and each curve was then checked against the classical solution of its own kind at the midpoint: Kepler's equation at e = 0.6 agrees to machine precision; Barker's cubic at e = 1 agrees to machine precision; e sinh H − H at e = 1.0001 agrees to machine precision; e sinh H − H at e = 1.4 agrees to machine precision. The curve to read twice is e = 1.0001: over this arc it lies within 0.02% of the parabola and is indistinguishable from it, and it is the only one of the four whose fate the drawing cannot show. What the classical parameterisation costs there is arithmetic rather than impossibility: at the midpoint of this arc, e sinh H − H throws away 3.3 of its sixteen digits to cancellation, against 0.2 at e = 1.4 — enough to matter to an ephemeris and not enough to stop a plot.
Fig. 7 Nine hundred days of the same four orbits, propagated by the same call. The ellipse comes back, the parabola and the hyperbolas do not, and nothing in the propagation had to know which was which in advance — the universal variable simply advances further along the ones that are unbound. That is the property the method is used for in practice: a fitted orbit whose eccentricity is 1.0002 ± 0.0004 does not have to be classified before it can be propagated, and classifying it would introduce a discontinuity the data does not support.

The problem where the class is not known in advance

Everything above assumes the orbit is in hand and the question is where the body will be. There is a second problem, at least as common, in which the conic class is not merely inconvenient to branch on — it is unknown until the answer is found.

Given two positions and the time to go between them, find the orbit. That is Lambert’s problem, and it is what a trajectory designer solves to connect a departure to an arrival, and what an orbit determination solves to turn two observations into a first guess.

The awkwardness is immediate. Fix the two positions and shorten the time of flight: the orbit connecting them is an ellipse, then a longer and flatter ellipse, then at one particular flight time exactly a parabola, and for anything faster a hyperbola. The conic class is a function of the input, and it is the input the designer is sweeping across — a porkchop plot is precisely a scan over departure and arrival dates, and its contours run straight through the parabolic case without pausing.

A solver built on the classical forms therefore has to guess the class before it can begin, discover it was wrong, and switch. Worse, the quantity being iterated on differs between the branches, so the switch discards the progress made. And the case where this happens is not exotic: the parabolic boundary sits in the middle of the useful region for fast transfers to the outer planets, which is exactly where the interesting trajectories are.

The universal formulation removes the question rather than answering it. The same variable χ\chi, the same Stumpff functions, and one iteration that walks from the elliptical side to the hyperbolic side without any code noticing which side it is on — because α\alpha passing through zero is not an event. A designer sweeping a grid of departure dates gets a surface that is smooth in the quantity being optimised, and a smooth surface is what every optimiser downstream requires.

That is the practical argument for the whole apparatus, and it is stronger than the argument from tidiness. A branch on the conic class is not merely ugly; it puts a discontinuity into the derivative of a function that something else is about to differentiate.

One analytic function each, and a case split that is only in the writing. Left: the Stumpff functions c₂ and c₃ across z = αχ², from z = −45 to 45. Positive z is an ellipse, negative z a hyperbola, and z = 0 exactly a parabola — and there is nothing at the origin: c₂ passes through 0.5 and c₃ through 0.166667 with every derivative continuous, because each is the single power series Σ(−z)^j/(k+2j)!. Written as (1 − cos√z)/z and (√z − sin√z)/z^3/2 they need three cases; written as their series they need none, and the case split in the code above is about floating point rather than about the mathematics. Right: how much floating point. The closed form subtracts two numbers that agree ever more closely as z → 0, so it loses digits — 11.9 of the sixteen it started with at z = 10⁻¹², measured against the series rather than estimated. That is why the drawn curve switches to the series inside |z| < 0.001, and why an implementation that uses the closed form everywhere is worst exactly at the parabola it was written to include.
Fig. 8 The two Stumpff functions over a wider range of the argument — forty-five either side of zero rather than thirty. Both are entire: they are power series that converge everywhere, and their values at the origin are 1/2 and 1/6 exactly. The removable singularity of the classical formulae is a hole in a function that has no hole in it, and this is the picture of that statement: the functions carry on smoothly through z=0z = 0 because they were defined by series rather than by ratios.

The gain from the change of variable is largest exactly where the ordinary variable is worst, and it is worth seeing at an eccentricity where a fixed time step is hopeless.

36 steps round one orbit, spaced two ways. The same orbit at e = 0.99, marked at 36 equal steps of time on the left and 36 equal steps of the fictitious time χ on the right. Equal time steps are what a fixed-step integrator takes, and they put the marks where the body is slow: the longest gap is 37.7 times the shortest, and the pericentre passage — the only part of the orbit where anything interesting happens to the trajectory — is crossed in 1 of them. Stepping in χ instead makes the time step proportional to r, and the marks come out 6.02-to-one even, with 1 inside the same pericentre region. Nobody chose that. dt/dχ = r/√μ is the substitution itself, and since vr equals the angular momentum at both apses, a fixed step in χ is very nearly a fixed step in arc length — an adaptive integrator's behaviour obtained from a change of variable rather than from a controller.
Fig. 9 The same construction at e = 0.99 with a periapsis at a tenth of the unit distance. On the left, thirty-six equal steps of time crowd into the far end of the orbit and leave the periapsis passage covered by one or two; on the right, thirty-six equal steps of the fictitious time are spread nearly evenly around the path. The integrator that uses the second grid resolves the encounter with the same effort it spends on the rest of the orbit.

That is the whole of what regularisation buys, and it is an adaptive step size obtained without any adaptation. A conventional integrator with a fixed step must be sized for the fastest part of the motion and therefore wastes almost all of its work on the slow part; one with a controller re-sizes the step by measuring an error estimate and paying for the measurement. The change of variable does it by construction, because the fictitious time advances in proportion to the distance and the distance is exactly the thing that sets how fast the motion is.

The cost is that the answer arrives on the wrong axis. A solution parameterised by χ\chi has to be converted back to tt to be useful, and that conversion is itself an equation to be solved. What makes it worthwhile is that the equation to be solved is smooth and monotonic in every regime, which is precisely what the original Kepler equation is not. The trade is therefore one awkward equation for two well-behaved ones, and it is worth taking because the awkwardness in the original is concentrated exactly where trajectories are most interesting — at close approach, at high eccentricity, near the parabolic boundary — while the awkwardness in the replacement is uniform and small.

Where the model stops

The universal formulation solves the two-body problem and nothing else. Everything it removes is a difficulty of parameterisation; nothing it removes is a difficulty of physics.

It does not help with a third body. It does not help with perturbations, which change the elements the solve is performed with. It has nothing to say about the Roche limit, tides, or radiation pressure. An orbit that is being pushed on is not a conic, and no parameterisation of conics will make it one. It also has a failure mode of its own, and it is the mirror image of the classical ones. The starting guess for χ\chi is not as well understood as the starting guess for EE on a nearly circular orbit, where EME \approx M is famously good. For a strongly hyperbolic orbit far from periapsis the Newton iteration on χ\chi can be slow to find its footing, and production implementations bracket it or fall back on a bisection first. The universal form removes a discontinuity and does not remove the need to think.

And there is a limit that is not numerical at all. On a hyperbola, χ\chi grows without bound as the body recedes; the products χ3c3(z)\chi^3 c_3(z) that appear in the equation are then differences of large quantities, and precision has to be watched again — for a different reason, in a different regime, by a different remedy. Every parameterisation is good over some range and no parameterisation is good over all of them. What the universal variable buys is not universality in the strong sense; it is one continuous family where there were three, with the seams moved away from the place the objects actually live.

One Newton loop, four conics, 6–12 steps. The residual |√μΔt − f(χ)| after each Newton step of the universal Kepler equation, at e = 0.6, e = 1, e = 1.0001, e = 1.4, all through the same code, the same starting rule and the same convergence test. Each orbit has periapsis at 0.5 AU and the solve is for 630.0 days past it. The step counts are 6 at e = 0.6, 11 at e = 1, 12 at e = 1.0001, 12 at e = 1.4, and the pair that matters is the crossing: 11 steps at e = 1 against 12 at e = 1.0001, a difference of 1. The spread across the whole set is 2.00, and it is a property of the STARTING GUESS rather than of the conic class — the run that begins nearest its own answer finishes first, on every branch. What the figure is about is the absence of a step at e = 1: the derivative of the equation being solved is r itself, which is positive on every conic, so the iteration has nothing to notice there. The classical forms have the opposite property — Kepler's equation degenerates as e → 1 from below, Barker's cubic is exact only AT 1, and e sinh H − H loses digits to cancellation just above it.
Fig. 10 The Newton iteration solving for a later point on the orbit — seven tenths of the way through the sweep rather than a third. The four conics still converge in six to twelve steps and the number of steps is nearly independent of which one is being solved. A solver whose cost does not depend on the case is the practical payoff, and it is what lets a trajectory integrator step across a near-parabolic passage without noticing it happened.

That is the whole argument for the universal formulation: not that it is faster, and not that it is more accurate, but that it removes a branch from code that has to run without supervision for decades.

What the picture cannot show

The figures here plot a radius against a time and a residual against a step. Three things they cannot carry:

The plane. Every curve is a distance from the primary, and a distance is one number. Turning it into a position needs three further angles, and the whole of the orientation problem sits outside these axes.

Which way the body is going. A radius-against-time curve is symmetric about periapsis for every conic here, so an inbound body and an outbound one give the same drawing. The distinction lives in the sign of χ\chi, which no axis on this page carries.

Whether the orbit is real. The e=1.0001e = 1.0001 curve is drawn to eleven digits and belongs to an object whose eccentricity, if it had one, would be known to perhaps three. The figure asserts a distinction between bound and unbound that the observations behind any such orbit could not support, and the amplification of that uncertainty is the reason marginal comets are catalogued as parabolas by convention.

Where the ladder goes next

Later rungs on this anchor: the ff and gg functions, which turn a solved χ\chi back into a position and velocity without ever forming the elements — the form in which the universal variable is actually used, and the reason a Lambert solver can be written in twenty lines. The Sundman transformation, of which dχ=μdt/rd\chi = \sqrt\mu\,dt/r is one case, and what regularising a collision singularity means. The Kustaanheimo–Stiefel transformation, which takes the same idea into four dimensions and makes the two-body problem into a harmonic oscillator exactly. Starting guesses, as a subject: the literature on how to begin a Newton iteration for χ\chi is longer than the literature on the iteration. And the accuracy of a sequence of solves, which is where an ephemeris lives and where losing three digits per call stops being invisible.

Barker’s equation was published in 1757 and is a cubic that Cardano could have solved two centuries earlier. Kepler’s is from 1609 and has no closed form at all. That the harder-looking of the two is the special case, and the impossible-looking one the general rule, is the sort of thing a single parameterisation makes look like an accident — and a better one makes look like the two ends of a curve.

What this makes readable

Essays that name this one as a prerequisite.

About the same objects

Not linked from either essay — found by the objects both name.

What links here

Essays that link to this one from their own argument.

The objects this essay names

Each one links to every other essay that touches it.

Barkers equationConditioningConic sectionsEccentricityHyperbolic orbitKepler's equationNewton's methodOrbital energySemi-major axisStumpff functionsUniversal variables