Universal Variable Formulation

Motivation

The classical anomaly equations differ according to the sign of orbital energy.

For an ellipse,

\[ M=E-e\sin E. \]

For a parabola, Barker’s equation is used.

For a hyperbola,

\[ M_h=e\sinh H-H. \]

This is not because the underlying dynamics change.

All three orbit classes satisfy the same two-body equation,

\[ \ddot{\mathbf r} = -\frac{\mu}{r^3}\mathbf r. \]

The difference arises from the coordinate used to parameterize motion along the conic.

Therefore, from a computational and conceptual point of view, it is natural to ask:

Can all three conic types be propagated with a single time variable and a single set of equations?

There is also a second motivation.

Suppose we are given a Cartesian state directly:

\[ \mathbf r_0,\qquad \mathbf v_0. \]

Do we really need to convert to orbital elements, solve a conic-specific anomaly equation, convert to true anomaly, build a perifocal state, and rotate back to inertial coordinates?

A direct state-vector propagator would be preferable.

This leads to the Lagrange \(f\) and \(g\) functions.

7. The Cartesian initial-value problem

At time \(t_0\), suppose

\[ \mathbf r(t_0)=\mathbf r_0, \]

and

\[ \mathbf v(t_0)=\mathbf v_0. \]

We want

\[ \mathbf r(t_0+\Delta t) \]

and

\[ \mathbf v(t_0+\Delta t). \]

The two-body force is central, so

\[ \mathbf h=\mathbf r\times\mathbf v \]

is conserved.

Thus the motion remains in a fixed plane.

Both

\[ \mathbf r_0 \]

and

\[ \mathbf v_0 \]

lie in this plane.

Assuming they are not collinear, they span the orbital plane.

Therefore every later position vector can be written as a linear combination of them:

\[ \boxed{ \mathbf r = f\mathbf r_0 + g\mathbf v_0. } \]

This is the fundamental Lagrange representation.

Since \(\mathbf r\) and \(\mathbf r_0\) have dimensions of length,

\[ [f]=1. \]

Since \(\mathbf v_0\) has dimensions of length per time,

\[ [g]=T. \]

Differentiating with respect to time,

\[ \boxed{ \mathbf v = \dot f\mathbf r_0 + \dot g\mathbf v_0. } \]

The quantities

\[ f,\qquad g,\qquad\dot f,\qquad\dot g \]

are called the Lagrange coefficients.

8. Initial conditions for \(f\) and \(g\)

At

\[ t=t_0, \]

we require

\[ \mathbf r=\mathbf r_0. \]

Thus,

\[ \mathbf r_0 = f(t_0)\mathbf r_0 + g(t_0)\mathbf v_0. \]

Hence,

\[ \boxed{ f(t_0)=1, \qquad g(t_0)=0. } \]

Similarly,

\[ \mathbf v_0 = \dot f(t_0)\mathbf r_0 + \dot g(t_0)\mathbf v_0. \]

Therefore,

\[ \boxed{ \dot f(t_0)=0, \qquad \dot g(t_0)=1. } \]

9. The identity \(f\dot g-\dot f g=1\)

Start with

\[ \mathbf r = f\mathbf r_0+g\mathbf v_0, \]

and

\[ \mathbf v = \dot f\mathbf r_0+\dot g\mathbf v_0. \]

Then

\[ \mathbf r\times\mathbf v = (f\mathbf r_0+g\mathbf v_0) \times (\dot f\mathbf r_0+\dot g\mathbf v_0). \]

Expanding,

\[ \mathbf r\times\mathbf v = f\dot g\,\mathbf r_0\times\mathbf v_0 + g\dot f\,\mathbf v_0\times\mathbf r_0, \]

because the cross products of a vector with itself vanish.

Since

\[ \mathbf v_0\times\mathbf r_0 = -\mathbf r_0\times\mathbf v_0, \]

we obtain

\[ \mathbf r\times\mathbf v = (f\dot g-g\dot f) \mathbf r_0\times\mathbf v_0. \]

Conservation of angular momentum requires

\[ \mathbf r\times\mathbf v = \mathbf r_0\times\mathbf v_0. \]

Therefore,

\[ \boxed{ f\dot g-\dot f g=1. } \]

This relation is an important analytical and numerical consistency check.

10. Why \(f\) and \(g\) do not yet solve the conic-dependence problem

The structure

\[ \mathbf r = f\mathbf r_0+g\mathbf v_0 \]

is universal.

However, the formulas used to evaluate \(f\) and \(g\) may still be written in terms of the eccentric anomaly for elliptic motion, hyperbolic anomaly for hyperbolic motion, and a separate limiting form for parabolic motion.

So the same structural problem remains:

\[ \boxed{ \text{different anomaly variables for different conics}. } \]

The goal is therefore stronger:

We want universal expressions for \(f\), \(g\), \(\dot f\), and \(\dot g\), driven by a single anomaly-like variable.

This motivates universal variables.

11. The energy parameter \(\alpha\)

The specific orbital energy is

\[ \varepsilon = \frac{v^2}{2} - \frac{\mu}{r}. \]

For a Keplerian orbit,

\[ \varepsilon = -\frac{\mu}{2a}. \]

At the initial state,

\[ \frac{v_0^2}{2} - \frac{\mu}{r_0} = -\frac{\mu}{2a}. \]

Multiply by \(2\):

\[ v_0^2-\frac{2\mu}{r_0} = -\frac{\mu}{a}. \]

Divide by \(\mu\):

\[ \frac{v_0^2}{\mu} - \frac{2}{r_0} = -\frac1a. \]

Hence,

\[ \boxed{ \alpha = \frac1a = \frac{2}{r_0} - \frac{v_0^2}{\mu}. } \]

Equivalently,

\[ \boxed{ \alpha = -\frac{2\varepsilon}{\mu}. } \]

Thus,

\[ \alpha>0 \quad\Longrightarrow\quad \text{ellipse}, \]

\[ \alpha=0 \quad\Longrightarrow\quad \text{parabola}, \]

and

\[ \alpha<0 \quad\Longrightarrow\quad \text{hyperbola}. \]

The orbit class is therefore encoded in the sign of one scalar.

12. Motivation for the universal anomaly \(\chi\)

For elliptic motion,

\[ M=E-e\sin E. \]

Also,

\[ M=n(t-\tau), \]

where

\[ n=\sqrt{\frac{\mu}{a^3}}. \]

Differentiate the Kepler equation:

\[ n\,dt = (1-e\cos E)\,dE. \]

But

\[ r=a(1-e\cos E). \]

Therefore,

\[ 1-e\cos E = \frac{r}{a}. \]

Hence,

\[ n\,dt = \frac{r}{a}\,dE. \]

Thus,

\[ dt = \frac{r}{an}\,dE. \]

Since

\[ n=\sqrt{\frac{\mu}{a^3}}, \]

we obtain

\[ an = a\sqrt{\frac{\mu}{a^3}} = \sqrt{\frac{\mu}{a}}. \]

Therefore,

\[ dt = \frac{r\sqrt a}{\sqrt\mu}\,dE. \]

This suggests defining a new variable \(\chi\) by

\[ d\chi=\sqrt a\,dE \]

for elliptic motion.

Then

\[ dt = \frac{r}{\sqrt\mu}\,d\chi. \]

This motivates the general definition

\[ \boxed{ \frac{d\chi}{dt} = \frac{\sqrt\mu}{r}. } \]

Equivalently,

\[ \boxed{ dt = \frac{r}{\sqrt\mu}\,d\chi. } \]

For an ellipse,

\[ \boxed{ \chi=\sqrt a\,\Delta E. } \]

But the definition itself makes no reference to a conic type.

That is what makes \(\chi\) a universal anomaly.

13. Units of \(\chi\)

Since

\[ [\mu]=\frac{L^3}{T^2}, \]

we have

\[ [\sqrt\mu]=\frac{L^{3/2}}{T}. \]

Therefore,

\[ \left[ \frac{\sqrt\mu}{r} \right] = \frac{L^{1/2}}{T}. \]

From

\[ \frac{d\chi}{dt} = \frac{\sqrt\mu}{r}, \]

it follows that

\[ \boxed{ [\chi]=L^{1/2}. } \]

This agrees with

\[ \chi=\sqrt a\,\Delta E. \]

14. Derivation of the universal radial equation

The energy equation is

\[ \varepsilon = \frac12 \left( \dot r^2+\frac{h^2}{r^2} \right) - \frac{\mu}{r}. \]

Using

\[ \varepsilon = -\frac{\mu\alpha}{2}, \]

we obtain

\[ -\frac{\mu\alpha}{2} = \frac12 \left( \dot r^2+\frac{h^2}{r^2} \right) - \frac{\mu}{r}. \]

Multiply by \(2\):

\[ -\mu\alpha = \dot r^2+\frac{h^2}{r^2}-\frac{2\mu}{r}. \]

Multiply by \(r^2\):

\[ -\mu\alpha r^2 = r^2\dot r^2+h^2-2\mu r. \]

From

\[ \frac{d\chi}{dt} = \frac{\sqrt\mu}{r}, \]

we have

\[ \dot r = \frac{dr}{d\chi} \frac{d\chi}{dt}. \]

Thus,

\[ \dot r = \frac{\sqrt\mu}{r} \frac{dr}{d\chi}. \]

Hence,

\[ r^2\dot r^2 = \mu \left( \frac{dr}{d\chi} \right)^2. \]

Substitute into the energy equation:

\[ -\mu\alpha r^2 = \mu \left( \frac{dr}{d\chi} \right)^2 +h^2-2\mu r. \]

Divide by \(\mu\):

\[ -\alpha r^2 = \left( \frac{dr}{d\chi} \right)^2 + \frac{h^2}{\mu} - 2r. \]

Therefore,

\[ \left( \frac{dr}{d\chi} \right)^2 = 2r-\alpha r^2-\frac{h^2}{\mu}. \]

Differentiate with respect to \(\chi\):

\[ 2\frac{dr}{d\chi} \frac{d^2r}{d\chi^2} = 2\frac{dr}{d\chi} - 2\alpha r\frac{dr}{d\chi}. \]

Away from a turning point, divide by

\[ 2\frac{dr}{d\chi}. \]

Then

\[ \frac{d^2r}{d\chi^2} = 1-\alpha r. \]

By continuity, this differential equation is valid through radial turning points as well.

Therefore,

\[ \boxed{ \frac{d^2r}{d\chi^2} + \alpha r = 1. } \]

This is the fundamental scalar equation underlying the universal formulation.

15. Initial conditions for the radial equation

At

\[ \chi=0, \]

we are at the initial state:

\[ \boxed{ r(0)=r_0. } \]

Also,

\[ \frac{dr}{d\chi} = \frac{dr}{dt}\frac{dt}{d\chi}. \]

Since

\[ \frac{dt}{d\chi} = \frac{r}{\sqrt\mu}, \]

we get

\[ \frac{dr}{d\chi} = \dot r\frac{r}{\sqrt\mu}. \]

At the initial point,

\[ r'(0) = \frac{r_0\dot r_0}{\sqrt\mu}. \]

But

\[ \dot r_0 = \frac{\mathbf r_0\cdot\mathbf v_0}{r_0}. \]

Therefore,

\[ \boxed{ r'(0) = \frac{\mathbf r_0\cdot\mathbf v_0}{\sqrt\mu}. } \]

16. Why Stumpff functions are needed

The radial equation is

\[ r''+\alpha r=1. \]

If

\[ \alpha>0, \]

the homogeneous equation has trigonometric solutions.

If

\[ \alpha<0, \]

the homogeneous equation has hyperbolic solutions.

If

\[ \alpha=0, \]

the solution is polynomial.

Thus a naive solution would still branch into three cases.

To avoid this, introduce the Stumpff functions.

Define

\[ \boxed{ C(z) = \sum_{k=0}^{\infty} \frac{(-z)^k}{(2k+2)!} } \]

and

\[ \boxed{ S(z) = \sum_{k=0}^{\infty} \frac{(-z)^k}{(2k+3)!}. } \]

The first few terms are

\[ C(z) = \frac12 - \frac{z}{24} + \frac{z^2}{720} -\cdots, \]

and

\[ S(z) = \frac16 - \frac{z}{120} + \frac{z^2}{5040} -\cdots. \]

Therefore,

\[ \boxed{ C(0)=\frac12, \qquad S(0)=\frac16. } \]

Define

\[ \boxed{ z=\alpha\chi^2. } \]

For \(z>0\),

\[ \boxed{ C(z) = \frac{1-\cos\sqrt z}{z}, } \]

and

\[ \boxed{ S(z) = \frac{\sqrt z-\sin\sqrt z}{z^{3/2}}. } \]

For \(z<0\),

\[ \boxed{ C(z) = \frac{\cosh\sqrt{-z}-1}{-z}, } \]

and

\[ \boxed{ S(z) = \frac{\sinh\sqrt{-z}-\sqrt{-z}} {(-z)^{3/2}}. } \]

Thus the same pair of functions represents elliptic, parabolic, and hyperbolic behavior continuously.

17. Figure: Stumpff functions

Code 
import numpy as np
import matplotlib.pyplot as plt
import math

def stumpff_C(z):
    z = np.asarray(z, dtype=float)
    C = np.empty_like(z)

    positive = z > 1e-8
    negative = z < -1e-8
    small = ~(positive | negative)

    C[positive] = (1 - np.cos(np.sqrt(z[positive]))) / z[positive]
    C[negative] = (
        np.cosh(np.sqrt(-z[negative])) - 1
    ) / (-z[negative])

    zs = z[small]
    C[small] = (
        1/2
        - zs/math.factorial(4)
        + zs**2/math.factorial(6)
        - zs**3/math.factorial(8)
    )
    return C

def stumpff_S(z):
    z = np.asarray(z, dtype=float)
    S = np.empty_like(z)

    positive = z > 1e-8
    negative = z < -1e-8
    small = ~(positive | negative)

    S[positive] = (
        np.sqrt(z[positive]) - np.sin(np.sqrt(z[positive]))
    ) / z[positive]**1.5

    S[negative] = (
        np.sinh(np.sqrt(-z[negative])) - np.sqrt(-z[negative])
    ) / (-z[negative])**1.5

    zs = z[small]
    S[small] = (
        1/6
        - zs/math.factorial(5)
        + zs**2/math.factorial(7)
        - zs**3/math.factorial(9)
    )
    return S

z = np.linspace(-5, 10, 700)

plt.figure(figsize=(8, 5))
plt.plot(z, stumpff_C(z), label="C(z)")
plt.plot(z, stumpff_S(z), label="S(z)")
plt.axvline(0, linewidth=1)
plt.xlabel("z")
plt.ylabel("Function value")
plt.grid(True)
plt.legend()
plt.show()
Figure 1: The Stumpff functions C(z) and S(z), showing smooth behavior through the parabolic limit z = 0.

The important point is that the parabolic limit

\[ z=0 \]

is not singular.

Instead,

\[ C(z)\rightarrow\frac12, \]

and

\[ S(z)\rightarrow\frac16. \]

18. Universal solution for the radius

The solution of

\[ r''+\alpha r=1 \]

satisfying the initial conditions can be written as

\[ \boxed{ r = \chi^2C(z) + \frac{\mathbf r_0\cdot\mathbf v_0}{\sqrt\mu} \chi[1-zS(z)] + r_0[1-zC(z)]. } \]

Introduce

\[ U_0=1-zC(z), \]

\[ U_1=\chi[1-zS(z)], \]

and

\[ U_2=\chi^2C(z). \]

Then

\[ r = r_0U_0 + \frac{\mathbf r_0\cdot\mathbf v_0}{\sqrt\mu}U_1 + U_2. \]

The useful derivative identities are

\[ \boxed{ \frac{d}{d\chi} \left[ \chi^2C(z) \right] = \chi[1-zS(z)] } \]

and

\[ \boxed{ \frac{d}{d\chi} \left[ \chi[1-zS(z)] \right] = 1-zC(z). } \]

Therefore,

\[ U_2'=U_1, \]

and

\[ U_1'=U_0. \]

At \(\chi=0\),

\[ U_0=1, \qquad U_1=0, \qquad U_2=0. \]

Hence,

\[ r(0)=r_0. \]

Also,

\[ r'(0) = \frac{\mathbf r_0\cdot\mathbf v_0}{\sqrt\mu}, \]

as required.

19. Derivation of the universal Kepler equation

By definition,

\[ dt = \frac{r}{\sqrt\mu}\,d\chi. \]

Integrate:

\[ \boxed{ \sqrt\mu\,\Delta t = \int_0^\chi r(u)\,du. } \]

Substitute the universal expression for \(r\):

\[ r(u) = u^2C(\alpha u^2) + \frac{\mathbf r_0\cdot\mathbf v_0}{\sqrt\mu} u[1-\alpha u^2S(\alpha u^2)] + r_0[1-\alpha u^2C(\alpha u^2)]. \]

Using the identities

\[ \int_0^\chi u^2C(\alpha u^2)\,du = \chi^3S(z), \]

\[ \int_0^\chi u[1-\alpha u^2S(\alpha u^2)]\,du = \chi^2C(z), \]

and

\[ \int_0^\chi [1-\alpha u^2C(\alpha u^2)]\,du = \chi[1-zS(z)], \]

we obtain

\[ \sqrt\mu\,\Delta t = \chi^3S(z) + \frac{\mathbf r_0\cdot\mathbf v_0}{\sqrt\mu} \chi^2C(z) + r_0\chi[1-zS(z)]. \]

Since

\[ z=\alpha\chi^2, \]

the last term becomes

\[ r_0\chi[1-zS(z)] = r_0\chi - \alpha r_0\chi^3S(z). \]

Therefore,

\[ \boxed{ \sqrt\mu\,\Delta t = \frac{\mathbf r_0\cdot\mathbf v_0}{\sqrt\mu} \chi^2C(z) + (1-\alpha r_0)\chi^3S(z) + r_0\chi. } \]

Using

\[ \mathbf r_0\cdot\mathbf v_0 = r_0v_{r0}, \]

we may also write

\[ \boxed{ \sqrt\mu\,\Delta t = \frac{r_0v_{r0}}{\sqrt\mu} \chi^2C(z) + (1-\alpha r_0)\chi^3S(z) + r_0\chi. } \]

This is the universal Kepler equation.

20. Why the universal Kepler equation replaces the three classical forms

The single equation above applies when

\[ \alpha>0, \]

\[ \alpha=0, \]

and

\[ \alpha<0. \]

Thus,

\[ \boxed{ \text{elliptic Kepler} + \text{Barker} + \text{hyperbolic Kepler} } \]

are replaced computationally by

\[ \boxed{ \text{one universal Kepler equation}. } \]

The orbit type is carried by

\[ \alpha \]

and therefore by

\[ z=\alpha\chi^2. \]

21. Newton solution of the universal Kepler equation

Define

\[ F(\chi) = \frac{r_0v_{r0}}{\sqrt\mu} \chi^2C(z) + (1-\alpha r_0)\chi^3S(z) + r_0\chi - \sqrt\mu\,\Delta t. \]

We seek

\[ F(\chi)=0. \]

From

\[ \sqrt\mu\,\Delta t = \int_0^\chi r(u)\,du, \]

the fundamental theorem of calculus gives

\[ \boxed{ F'(\chi)=r(\chi). } \]

Therefore Newton iteration is

\[ \chi_{n+1} = \chi_n - \frac{F(\chi_n)}{F'(\chi_n)}. \]

Hence,

\[ \boxed{ \chi_{n+1} = \chi_n - \frac{F(\chi_n)}{r(\chi_n)}. } \]

This is one of the elegant computational features of the universal formulation.

22. Universal Lagrange \(f\) and \(g\) functions

Once \(\chi\) is known,

\[ \boxed{ f = 1-\frac{\chi^2}{r_0}C(z) } \]

and

\[ \boxed{ g = \Delta t - \frac{\chi^3}{\sqrt\mu}S(z). } \]

Thus,

\[ \boxed{ \mathbf r = \left[ 1-\frac{\chi^2}{r_0}C(z) \right]\mathbf r_0 + \left[ \Delta t-\frac{\chi^3}{\sqrt\mu}S(z) \right]\mathbf v_0. } \]

At

\[ \chi=0, \]

\[ f=1, \]

and

\[ g=0, \]

so the initial position is recovered.

23. Connection with the elliptic \(f\) and \(g\) functions

For elliptic motion,

\[ \alpha=\frac1a, \]

and

\[ \chi=\sqrt a\,\Delta E. \]

Therefore,

\[ z=\alpha\chi^2 = \Delta E^2. \]

Since

\[ C(z) = \frac{1-\cos\Delta E}{\Delta E^2}, \]

we get

\[ \chi^2C(z) = a(1-\cos\Delta E). \]

Therefore,

\[ \boxed{ f = 1-\frac{a}{r_0}(1-\cos\Delta E). } \]

Similarly,

\[ S(z) = \frac{\Delta E-\sin\Delta E}{\Delta E^3}, \]

so

\[ \chi^3S(z) = a^{3/2}(\Delta E-\sin\Delta E). \]

Thus,

\[ \boxed{ g = \Delta t - \sqrt{\frac{a^3}{\mu}} (\Delta E-\sin\Delta E). } \]

Hence the universal expressions reduce exactly to the familiar elliptic Lagrange coefficients.

24. Derivation of \(\dot f\)

Start from

\[ f = 1-\frac{\chi^2}{r_0}C(z). \]

Differentiate with respect to time:

\[ \dot f = -\frac1{r_0} \frac{d}{dt} [\chi^2C(z)]. \]

Using

\[ \frac{d}{dt} = \frac{d\chi}{dt}\frac{d}{d\chi}, \]

and

\[ \frac{d\chi}{dt} = \frac{\sqrt\mu}{r}, \]

we obtain

\[ \dot f = -\frac{\sqrt\mu}{rr_0} \frac{d}{d\chi} [\chi^2C(z)]. \]

Since

\[ \frac{d}{d\chi} [\chi^2C(z)] = \chi[1-zS(z)], \]

we have

\[ \dot f = -\frac{\sqrt\mu}{rr_0} \chi[1-zS(z)]. \]

Hence,

\[ \boxed{ \dot f = \frac{\sqrt\mu}{rr_0} [zS(z)-1]\chi. } \]

25. Derivation of \(\dot g\)

Start from

\[ g = \Delta t - \frac{\chi^3}{\sqrt\mu}S(z). \]

Differentiate:

\[ \dot g = 1 - \frac1{\sqrt\mu} \frac{d}{dt} [\chi^3S(z)]. \]

Using

\[ \frac{d}{dt} = \frac{\sqrt\mu}{r}\frac{d}{d\chi}, \]

we obtain

\[ \dot g = 1 - \frac1r \frac{d}{d\chi} [\chi^3S(z)]. \]

Since

\[ \frac{d}{d\chi} [\chi^3S(z)] = \chi^2C(z), \]

we get

\[ \boxed{ \dot g = 1-\frac{\chi^2}{r}C(z). } \]

26. Complete universal state propagation

Given

\[ \mathbf r_0, \qquad \mathbf v_0, \qquad \Delta t, \qquad \mu, \]

compute

\[ r_0=\|\mathbf r_0\|, \]

\[ v_0=\|\mathbf v_0\|, \]

\[ v_{r0} = \frac{\mathbf r_0\cdot\mathbf v_0}{r_0}, \]

and

\[ \alpha = \frac{2}{r_0} - \frac{v_0^2}{\mu}. \]

Solve

\[ \sqrt\mu\,\Delta t = \frac{r_0v_{r0}}{\sqrt\mu} \chi^2C(z) + (1-\alpha r_0)\chi^3S(z) + r_0\chi, \]

where

\[ z=\alpha\chi^2. \]

Then calculate

\[ f = 1-\frac{\chi^2}{r_0}C(z), \]

and

\[ g = \Delta t - \frac{\chi^3}{\sqrt\mu}S(z). \]

Propagate the position:

\[ \boxed{ \mathbf r = f\mathbf r_0+g\mathbf v_0. } \]

Then compute

\[ r=\|\mathbf r\|. \]

Next,

\[ \dot f = \frac{\sqrt\mu}{rr_0} [zS(z)-1]\chi, \]

and

\[ \dot g = 1-\frac{\chi^2}{r}C(z). \]

Finally,

\[ \boxed{ \mathbf v = \dot f\mathbf r_0+\dot g\mathbf v_0. } \]

27. Parabolic limit

For parabolic motion,

\[ \alpha=0. \]

Hence,

\[ z=0. \]

Using

\[ C(0)=\frac12, \]

and

\[ S(0)=\frac16, \]

the universal Kepler equation becomes

\[ \boxed{ \sqrt\mu\,\Delta t = r_0\chi + \frac{r_0v_{r0}}{2\sqrt\mu}\chi^2 + \frac16\chi^3. } \]

Thus the parabolic solution emerges smoothly from the same formulation.

28. Hyperbolic limit

For hyperbolic motion,

\[ \alpha<0, \]

so

\[ z<0. \]

The Stumpff functions automatically become hyperbolic:

\[ C(z) = \frac{\cosh\sqrt{-z}-1}{-z}, \]

and

\[ S(z) = \frac{\sinh\sqrt{-z}-\sqrt{-z}} {(-z)^{3/2}}. \]

No separate derivation of the state propagator is required.

29. Figure: the two propagation philosophies

Code 
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(11, 5))
ax.axis("off")

classical = [
    (0.07, 0.72, r"$\Delta t$"),
    (0.28, 0.72, r"$E,\ H,$ or parabolic anomaly"),
    (0.50, 0.72, r"$f$"),
    (0.70, 0.72, r"$\mathbf{r}_{PQW},\mathbf{v}_{PQW}$"),
    (0.91, 0.72, r"$\mathbf{r}_I,\mathbf{v}_I$")
]

universal = [
    (0.07, 0.28, r"$\Delta t$"),
    (0.32, 0.28, r"$\chi$"),
    (0.59, 0.28, r"$f,g,\dot f,\dot g$"),
    (0.86, 0.28, r"$\mathbf{r},\mathbf{v}$")
]

for x, y, text in classical:
    ax.text(
        x, y, text,
        ha="center", va="center",
        bbox=dict(boxstyle="round,pad=0.35", fill=False)
    )

for x, y, text in universal:
    ax.text(
        x, y, text,
        ha="center", va="center",
        bbox=dict(boxstyle="round,pad=0.35", fill=False)
    )

for (x1, y1, _), (x2, y2, _) in zip(classical[:-1], classical[1:]):
    ax.annotate(
        "", xy=(x2-0.06, y2), xytext=(x1+0.06, y1),
        arrowprops=dict(arrowstyle="->")
    )

for (x1, y1, _), (x2, y2, _) in zip(universal[:-1], universal[1:]):
    ax.annotate(
        "", xy=(x2-0.06, y2), xytext=(x1+0.06, y1),
        arrowprops=dict(arrowstyle="->")
    )

ax.text(0.02, 0.88, "Classical element-based propagation", fontsize=12)
ax.text(0.02, 0.44, "Universal state-vector propagation", fontsize=12)

plt.show()
Figure 2: Element-based and universal state-vector propagation routes.

30. What each idea contributes

The Kepler, Barker, and hyperbolic Kepler equations answer

\[ \boxed{ \text{Given time, where am I along the conic?} } \]

The perifocal frame answers

\[ \boxed{ \text{How do I write that conic location as vectors?} } \]

The orbital-angle transformation answers

\[ \boxed{ \text{How is the orbital plane oriented in inertial space?} } \]

The Lagrange functions answer

\[ \boxed{ \text{Can I propagate the Cartesian state directly?} } \]

The universal-variable formulation answers

\[ \boxed{ \text{Can that direct propagation use one set of equations for all conics?} } \]

32. Final conceptual summary

The entire development can be summarized as

\[ \boxed{ \text{time-of-flight equations} } \]

\[ \downarrow \]

\[ \boxed{ \text{location along a conic} } \]

\[ \downarrow \]

\[ \boxed{ \text{perifocal vector representation} } \]

\[ \downarrow \]

\[ \boxed{ \text{direct state propagation with } f \text{ and } g } \]

\[ \downarrow \]

\[ \boxed{ \text{unification through } \chi,\alpha,C(z),S(z) } \]

The final universal propagation chain is

\[ \boxed{ (\mathbf r_0,\mathbf v_0,\Delta t) \rightarrow \alpha \rightarrow \chi \rightarrow C(z),S(z) \rightarrow f,g,\dot f,\dot g \rightarrow (\mathbf r,\mathbf v). } \]

The central pedagogical message is:

Universal variables are not a new orbit theory. They are a unified parameterization of the same Kepler two-body dynamics that students have already encountered through elliptic, parabolic, and hyperbolic anomaly equations.