Four numbers that carry a state without an orbit

The universal variable answers a question nobody asks directly. It is a number, χ, related to the time since some starting moment by an equation that does not care whether the orbit is an ellipse, a parabola or a hyperbola, and solving that equation is cheap and robust. But χ is not a position. What a navigator, an ephemeris or an orbit determination wants is where the body is and how fast it is moving, and the step that turns a solved χ into those is the one that makes the method worth using.

It is four numbers. Given the position $\mathbf r_0$ and velocity $\mathbf v_0$ at the start, the state at a later time is

$$\mathbf r = f\,\mathbf r_0 + g\,\mathbf v_0, \qquad \mathbf v = \dot f\,\mathbf r_0 + \dot g\,\mathbf v_0,$$

with $f$, $g$, $\dot f$ and $\dot g$ simple functions of χ, the starting distance and the time. No semi-major axis is formed, no eccentricity, no anomaly, no node, no argument of periapsis. The orbit is carried forward as a linear map on the state, and the map’s four coefficients are all that has to be computed.

Why the state is linear in the state #

The linearity looks like a trick and is a consequence of the geometry. A two-body orbit lies in a plane, and the plane contains the starting position and the starting velocity. Any later position is in that plane too, so it can be written as some combination of the two starting vectors — they span the plane, unless they happen to be parallel, which happens only on a radial trajectory with no angular momentum at all. The coefficients of that combination are what $f$ and $g$ are. The same holds for the later velocity.

So the content of the coefficients is not the linearity, which is free, but their values, and those come from the universal-variable solution:

$$f = 1 - \frac{\chi^2}{r_0} C(z), \quad g = \Delta t - \frac{\chi^3}{\sqrt\mu} S(z), \quad \dot g = 1 - \frac{\chi^2}{r} C(z), \quad \dot f = \frac{\sqrt\mu}{r\,r_0}\,\chi\,\big(z S(z) - 1\big),$$

with $z = \chi^2/a$ and $C$ and $S$ the Stumpff functions that become cosines and sines for an ellipse, hyperbolic cosines and sines for a hyperbola, and a polynomial for the parabola. Nothing in the four formulas distinguishes the three cases; the sign of $1/a$ inside $z$ does it all.

That is why the coefficients are how the universal variable is actually used. The six classical elements are undefined for orbits that are perfectly ordinary — a circular orbit has no periapsis, an equatorial one no node — and converting a state to elements and back passes through those undefined angles whenever the orbit is near one of those cases. The coefficients never form them. A circular, equatorial, parabolic orbit is propagated by the same four lines as any other.

The identity that is angular momentum #

The four numbers are not independent. They satisfy

$$f\,\dot g - \dot f\,g = 1$$

exactly, for every orbit and every time step. The reason is a conservation law. The angular momentum per unit mass is $\mathbf r \times \mathbf v$. Substituting the linear map,

$$\mathbf r \times \mathbf v = (f\dot g - \dot f g)\,(\mathbf r_0 \times \mathbf v_0),$$

because the cross terms $\mathbf r_0\times\mathbf r_0$ and $\mathbf v_0\times\mathbf v_0$ vanish. Angular momentum is conserved, so the bracket is one. Geometrically, the map from the starting position and velocity to the later ones preserves the area of the parallelogram they span — which is Kepler’s equal areas written as a determinant.

The identity is the smallest case of a much more general property. The map that carries any mechanical system’s positions and momenta forward in time preserves a particular measure of area in its state space — the map is symplectic — and for a single body moving in a plane that measure reduces to exactly this two-by-two determinant. Numerical integrators built to preserve the property, rather than merely to be accurate over each step, are the ones that hold planetary orbits for billions of years without their energy drifting, because they are wrong about where a body is while staying right about how much energy it has. The Lagrange coefficients are the exact symplectic map for the two-body problem, and the identity is the one-line test of whether a computed map still is one.

In the first figure the identity holds to a few parts in $10^{15}$ at every one of four hundred instants round the orbit, which is the limit of double-precision arithmetic. That is a check on the computation rather than a property of the figure. It becomes useful when the computation is not exact.

Where g passes through zero #

The first figure has a feature worth stopping at. The scaled coefficient $g$ starts at zero, rises, and comes back through zero exactly half a period later, at apoapsis — and on any orbit it crosses zero whenever the body is directly in line with where it started, on the same side of the centre or the opposite one. At those instants the later position is a pure multiple of the starting position, $\mathbf r = f\mathbf r_0$, and the starting velocity contributes nothing to it.

That is harmless going forwards and fatal going backwards. Lambert’s problem asks for the starting velocity given both positions, and the coefficients give it as $(\mathbf r - f\mathbf r_0)/g$. When the two positions are on a line through the centre, $g$ is zero, the formula divides by it, and the starting velocity is undetermined: any orbit in any plane containing that line connects the two points. That is the singularity at a transfer angle of 180 degrees that no change of variable removes, and that appears in mission design as a wall across the porkchop plot along exactly the line where the transfer goes half-way round the Sun. Seen through the coefficients it is not a special pathology of Lambert’s problem. It is the moment in every orbit where the position stops depending on the velocity it started with.

The coefficients degenerate in one other place, and it is the same kind of place. If the starting position and velocity are parallel — a body falling straight in or rising straight out, with no angular momentum — they no longer span a plane, the later position cannot be written as a combination of two vectors that point the same way unless it too points that way, and the identity’s right-hand side, the angular momentum, is zero. A radial fall is the limit of an ellipse of no width, and in that limit the linear map collapses onto a line. The universal variable itself is fine there; it is the representation of the state as a combination of two vectors that runs out of vectors.

A series that says how wrong it is #

The coefficients have a second form, older than the universal variable and still in use. Expand $f$ and $g$ as power series in the time step, using only the starting state: Gauss did this in 1809, in the method that recovered the dwarf planet Ceres from a few weeks of observations, and every orbit determination from three sightings since has begun the same way.

The series is what makes an orbit from three observations tractable at all. With only directions measured and no distances, the unknown positions at the three times can be written in terms of the middle one by truncated $f$ and $g$, and the problem reduces to a single equation for the middle distance. The price is the truncation. The error of a truncated series grows as a power of the time step, one more than the order kept, and the drawn slopes come out at exactly three, four and five. A series kept to the third power of the step and applied over a hundredth of an orbit misplaces the body by three parts in a million — fine for a first orbit of an asteroid seen over a few nights, ruinous over a few months.

The episode that made the series famous is worth recalling for what it shows about their range. Ceres was discovered on the first night of 1801 and followed for forty-one days before it was lost in the Sun’s glare, having moved through only about three degrees of its orbit. The astronomers of the day expected to find it again by assuming a circular orbit, and failed. Gauss, twenty-four, fitted an orbit to three of the observations by the series method and an iteration, predicted where Ceres would reappear at the end of the year, and it was found within half a degree of his position. Forty-one days is about a fortieth of Ceres’ period — inside the range where the figure shows the third-order series accurate to a part in a hundred thousand, which is where the iteration needs it to be.

The dashed curves are the point of the figure. Each is how badly the truncated set of four numbers misses the identity, computed from those four numbers alone, and each rises at the same rate as the true error beside it. For the third- and fourth-order series it sits within a factor of a few of the error; for the second order it sits a few hundred times above it, because that truncation drops terms that matter to the identity more than to the position. Either way the ratio is fixed for a given method. Calibrated once, the identity is a meter that reads the error of every subsequent use without a reference solution to compare against. It tells a calculation when its sightings are too far apart for its own approximations.

The identity as an error meter #

The same reading applies to the exact formulas when they are computed inexactly. The universal-variable solve is an iteration, and stopping it early leaves χ slightly wrong; the four coefficients computed from that χ then describe a slightly wrong state.

A hundred and twenty propagations — three step lengths, forty tolerances — lie along one line of unit slope. However the solve was cut short and however far the step reached, the identity’s violation equals the position error to within a factor of order one. A propagator that computes the bracket after every step, at the cost of two multiplications and a subtraction, knows how accurate that step was. In a long chain of steps — an ephemeris integrated for a century, a navigation filter running for the length of a mission — that is the difference between trusting a number and checking it.

There is a limit to what the meter sees. An error in the starting state, or in the gravitational parameter, propagates faithfully into a state that satisfies the identity perfectly: the wrong orbit, carried forward exactly. The identity checks the arithmetic of the map, not the physics fed into it. An ephemeris is a fit, and the errors in its fitted starting conditions are invisible to any check inside the propagation.

One formula through the parabola #

The coefficients carry the universal variable’s main virtue with them.

In the classical treatment these are three calculations. The ellipse uses Kepler’s equation, the hyperbola its hyperbolic analogue, and the parabola Barker’s cubic, and each formula approaches its neighbour badly: Kepler’s equation loses digits as the eccentricity approaches one, and the hyperbolic form does the same from the other side. A comet whose orbit is within a part in a thousand of parabolic — and most long-period comets are — sits exactly where both classical forms are worst. In the coefficients there is nothing to approach: the three curves are one function evaluated at three values of $1/a$, the middle one is simply the middle one, and a body that changes from bound to unbound in the middle of a calculation, as a comet’s orbit can under a planetary perturbation, does not change formula.

Many small steps lose what one long step keeps #

The linearity has one more consequence, and it bears on how a propagation should be organised.

A numerical integrator must take small steps, because it approximates the motion over each step by a polynomial. A two-body propagation by the coefficients need not, because each step is exact up to the accuracy of its own solve: one step of forty orbits is as good as one step of an hour. Chaining steps then has no benefit and a cost. Each step’s small error places the body slightly wrong, the next step starts from the misplaced state and inherits the error, and the errors accumulate rather than resetting. In the figure, ten times as many steps, each solved as accurately as before, end twelve times further off.

That is the place where, as the regularised form of the same problem showed for collisions, the choice of variables decides the accuracy rather than merely the convenience. A Kepler propagation that can jump straight to the time wanted pays for one solve’s error. A perturbed propagation cannot jump — the planets’ pulls change the orbit continuously — but the best methods for nearly-Keplerian motion use exactly this structure: they propagate the two-body part by the coefficients over long steps and treat the perturbations as small corrections between them, which is why such methods hold an orbit for millions of years where a general-purpose integrator drifts in thousands.

What the coefficients are for, in practice #

Three uses follow from the properties drawn above, and between them they cover most of what the universal variable is used for.

Orbit determination begins from the series, because the series needs no orbit yet: with three directions and two time intervals, the truncated coefficients turn the problem into one for a single unknown distance, and the identity says when the intervals are too long for the truncation. The answer is then refined by replacing the series with the exact coefficients and iterating.

Lambert’s problem — finding the orbit that takes a body from one position to another in a given time — is solved with the exact coefficients run backwards: given both positions and the time, the coefficients give the starting velocity directly, $\mathbf v_0 = (\mathbf r - f\mathbf r_0)/g$, once χ is found. The twenty-line solvers that every mission design uses are this formula with a good iteration for χ.

And state propagation in navigation and tracking uses the coefficients because they also give the sensitivity of the later state to the earlier one. The matrix of partial derivatives that says how an error in the starting position becomes an error later — the quantity that decides how a trajectory must be corrected — can be written in closed form from the same four numbers and their derivatives with respect to χ, with no numerical differentiation at all.

What the picture leaves out #

Everything here is the two-body problem in a plane. The coefficients extend to three dimensions unchanged, because the orbit’s plane is still spanned by the starting position and velocity; that is part of their appeal. They do not extend to perturbed motion, where the orbit’s plane and shape change continuously, and there the linear map holds only over steps short enough that the perturbation can be neglected, which is the regime the chain figure deliberately explores.

The error meter has a blind spot already named: it sees errors in the propagation, not in its inputs. And the per-step errors in the chain figure come from stopping a solve early, which is a choice; with every solve taken to the limit of the arithmetic, the remaining errors are rounding, and how those accumulate over a long chain depends on details of floating-point arithmetic that no figure at this scale can show.

Still open: errors that the identity cannot see #

The identity is one scalar condition on four numbers, so a set of coefficients can satisfy it and still be wrong — three degrees of freedom remain unchecked. For the errors that arise in practice, from a truncated series or an incomplete solve, the figures show that the identity catches them, because those errors break it in proportion to their size. Whether that holds for every source of error — rounding in nearly parabolic orbits, where the Stumpff functions lose digits, or errors in the time itself — has not been shown in general, and a propagation that passes its own check while drifting in a direction the check is blind to is the failure that long-running navigation has to guard against by other means.