Universal variable formulation

A common problem in orbital mechanics is the following: given a body in an orbit and a time t₀, find the position of the body at any other given time t. For elliptical orbits with a reasonably small eccentricity, solving Kepler's Equation by methods like Newton's method gives adequate results. However, as the orbit becomes more and more eccentric, the numerical iteration may start to converge slowly or not at all.[1][2] Furthermore, Kepler's equation cannot be applied to parabolic and hyperbolic orbits, since it specifically is tailored to elliptic orbits.

Although equations similar to Kepler's equation can be derived for parabolic and hyperbolic orbits, it is more convenient to introduce a new independent variable to take the place of the eccentric anomaly E, and having a single equation that can be solved regardless of the eccentricity of the orbit. The new variable s is defined by the following differential equation:

${\displaystyle {\frac {ds}{dt}}={\frac {1}{r}}}$

where ${\displaystyle r=r(t)}$ is the time-dependent distance to the center of attraction. The fundamental equation ${\displaystyle {\frac {d^{2}\mathbf {r} }{dt^{2}}}+\mu {\frac {\mathbf {r} }{r^{3}}}=\mathbf {0} }$ is regularized by applying this change of variables to yield:[2]

${\displaystyle {\frac {d^{2}\mathbf {r} }{ds^{2}}}+\alpha \ \mathbf {r} =-\mathbf {P} }$

where P is a constant vector and ${\displaystyle \alpha }$ is defined by

${\displaystyle \alpha ={\frac {\mu }{a}}}$

The equation is the same as the equation for the harmonic oscillator, a well-known equation in both physics and mathematics. Taking the derivative again, we get a third-degree differential equation:

${\displaystyle {\frac {d^{3}\mathbf {r} }{ds^{3}}}+\alpha {\frac {d\mathbf {r} }{ds}}=\mathbf {0} }$

The family of solutions to this differential equation[2] are written symbolically as the functions ${\displaystyle 1,\ s\ c_{1}(\alpha s^{2}),\ s^{2}\ c_{2}(\alpha s^{2}),}$ where the functions ${\displaystyle \ c_{k}(x)}$, called Stumpff functions, are generalizations of sine and cosine functions. Applying this results in:[3]

${\displaystyle t-t_{0}=r_{0}\ s\ c_{1}(\alpha s^{2})+r_{0}{\frac {dr_{0}}{dt}}\ s^{2}\ c_{2}(\alpha s^{2})+\mu \ s^{3}\ c_{3}(\alpha s^{2})}$

which is the universal variable formulation of Kepler's Equation. This equation can now be solved numerically using a root-finding algorithm such as Newton's method or Laguerre's method for a given time ${\displaystyle t}$ to yield ${\displaystyle s}$, which in turn is used to compute the f and g functions:

${\displaystyle {\begin{aligned}f(s)&=1-\left({\frac {\mu }{r_{0}}}\right)s^{2}c_{2}(\alpha s^{2}),\\g(s)&=t-t_{0}-\mu s^{3}c_{3}(\alpha s^{2}),\\{\frac {df}{dt}}&={\dot {f}}(s)=-\left({\frac {\mu }{rr_{0}}}\right)sc_{1}(\alpha s^{2}),\\{\frac {dg}{dt}}&={\dot {g}}(s)=1-\left({\frac {\mu }{r}}\right)s^{2}c_{2}(\alpha s^{2})\end{aligned}}}$

The values of the f and g functions determine the position of the body at the time ${\displaystyle t}$:

${\displaystyle \mathbf {r} =\mathbf {r} _{0}\ f(s)+\mathbf {v} _{0}\ g(s)}$

In addition the velocity of the body at time ${\displaystyle t}$ can be found using ${\displaystyle {\dot {f}}(s)}$ and ${\displaystyle {\dot {g}}(s)}$ as follows:

${\displaystyle \mathbf {v} =\mathbf {r} _{0}\ {\dot {f}}(s)+\mathbf {v} _{0}\ {\dot {g}}(s)}$

where ${\displaystyle \mathbf {r} }$ and ${\displaystyle \mathbf {v} }$ are the position and velocity respectively at time ${\displaystyle t}$, and ${\displaystyle \mathbf {r} _{0}}$ and ${\displaystyle \mathbf {v} _{0}}$ are the position and velocity, respectively, at arbitrary initial time ${\displaystyle t_{0}}$.