// Workers AI · dad joke modeWhat did the Stumpff function say? "I'm rooted in math.

In celestial mechanics, the Stumpff functions $\ c_{k}(x)\ ,$ were developed by Karl Stumpff for analyzing trajectories and orbits using the universal variable formulation. They are defined by the alternating series:

$\ c_{k}(x)~\equiv ~{\frac {1}{\ k!\ }}-{\frac {x}{~\left(k+2\right)!\ }}+{\frac {\ x^{2}}{~\left(k+4\right)!\ }}-\cdots ~=~\sum _{n=0}^{\infty }\ {\frac {\ (-1)^{n}\ x^{n}\ }{\ \left(k+2n\right)!\ }}~~$ for $~~k=0,1,2,3,\ \ldots ~~.$

Like the sine, cosine, and exponential functions, Stumpff functions are well-behaved entire functions : Their series converge absolutely for any finite argument $\ x~.$

Stumpff functions are useful for working with surface launch trajectories, and boosts from closed orbits to escape trajectories, since formulas for spacecraft trajectories using them smoothly meld from conventional closed orbits (circles and ellipses, eccentricity e : 0 ≤ e < 1 ) to open orbits (parabolas and hyperbolas, ( e ≥ 1 ), with no singularities and no imaginary numbers arising in the expressions as the launch vehicle gains speed to escape velocity and beyond. (The same advantage occurs in reverse, as a spacecraft decelerates from an arrival trajectory to go into a closed orbit around its destination, or descends to a planet's surface from a stable orbit.)

Relations to circular and hyperbolic trigononometric functions

By comparing the Taylor series expansion of the trigonometric functions sin and cos with $\ c_{0}(x)\$ and $\ c_{1}(x)\ ,$ a relationship can be found. For ${\displaystyle ~x>0\$ :}

${\begin{aligned}c_{0}(x)~&=~~\cos {\sqrt {x\ }}\ ,\\[1ex]c_{1}(x)~&=~{\frac {\ \sin {\sqrt {x\ }}\ }{\sqrt {x\ }}}~.\end{aligned}}$

Similarly, by comparing with the expansion of the hyperbolic functions sinh and cosh we find for ${\displaystyle ~x<0\$ :}

${\begin{aligned}c_{0}(x)~&=~~\cosh {\sqrt {-x\ }}\ ,\\[1ex]c_{1}(x)~&=~{\frac {\ \sinh {\sqrt {-x\ }}\ }{\sqrt {-x\ }}}~.\end{aligned}}$

Circular orbits and elliptical orbits use sine and cosine relations, and hyperbolic orbits use the sinh and cosh relations. Parabolic orbits (marginal escape orbits) formulas are a special in-between case.

Note that $\ \sin(x)\,/\,x\ =\operatorname {sinc} (x)$, the sine cardinal function; and similarly for sinh as the hyperbolic sine cardinal, sinhc.
By definition $\ \operatorname {sinc} (0)=\operatorname {sinhc} (0)=1\$.

Therefore: $${\begin{aligned}c_{0}(x)&={\begin{cases}\cos {\sqrt {x\ }}\quad &{\text{if }}x\geq 0,\\\cosh {\sqrt {-x\ }}\quad &{\text{if }}x<0.\end{cases}}\\[2ex]c_{1}(x)&={\begin{cases}\operatorname {sinc} {\sqrt {x\ }}\quad &{\text{if }}x\geq 0,\\\operatorname {sinhc} {\sqrt {-x\ }}\quad &{\text{if }}x<0.\end{cases}}\end{aligned}}$$

Recursion

For higher-order Stumpff functions needed for both ordinary trajectories and for perturbation theory, one can use the recurrence relation:

$\ x\ c_{k+2}(x)={\frac {1}{k!}}-c_{k}(x)~$ for $~k=0,1,2,\ \ldots \ ~,$

or when $\ x\neq 0\$

$\ c_{k+2}(x)~=~{\frac {\ 1\ }{x}}\left(\ {\frac {1}{k!}}-c_{k}(x)\ \right)~$ for $~k=0,1,2,\ \ldots ~~.$

Using this recursion, the two further Stumpff functions needed for the universal variable formulation are: $${\begin{aligned}c_{2}(x)&={\begin{cases}{\dfrac {1~-~\cos {\sqrt {x\ }}\ }{x}}\quad &{\text{if }}x>0,\\~{\tfrac {1}{2}}\quad &{\text{if }}x=0,\\{\dfrac {1~-~\cosh {\sqrt {-x\ }}\ }{x}}\quad &{\text{if }}x<0.\end{cases}}\\[2ex]c_{3}(x)&={\begin{cases}{\dfrac {1~-~\operatorname {sinc} {\sqrt {x\ }}\ }{x}}\quad &{\text{if }}x>0,\\~{\tfrac {1}{6}}\quad &{\text{if }}x=0,\\{\dfrac {1~-~\operatorname {sinhc} {\sqrt {-x\ }}\ }{x}}\quad &{\text{if }}x<0.\end{cases}}\end{aligned}}$$

Software implementation

By the small-angle approximation, both $\cos {\sqrt {x}}$ and $\operatorname {sinc} {\sqrt {x}}$ are close to 1 for small $x$, as are their hyperbolic counterparts $\cosh {\sqrt {-x}}$ and $\operatorname {sinhc} {\sqrt {-x}}$; $\cos {\sqrt {x}}$ is also close to 1 for $x=(4\pi ^{2})n^{2}$. Subtracting these from 1 is vulnerable to catastrophic cancellation.

For $x\leq (k+1)(k+2)/2$ the partial sum of the power series of $c_{k}(x)$ can be evaluated as a polynomial without severe loss of precision. However the larger $|x|$ the more terms are required to achieve a given accuracy.
If $|x|\leq \epsilon \,(k+1)(k+2)/2$ then only the one constant term of the polynomial is required, where $\epsilon$ is the machine epsilon.

The subtraction $(1-\cos \theta )$ can be calculated without cancellation as the versine function: $$\operatorname {versin} \theta =1-\cos \theta =2\sin ^{2}{\frac {\theta }{2}}$$ although this only gains significant precision for $\theta \approx (2\pi )n$.
Similarly: $$1-\cosh \theta =-2\sinh ^{2}{\frac {\theta }{2}}$$

A precise software implementation will blend these techniques in appropriate ranges based on a numerical analysis of the function and the machine's arithmetic precision.

Relations to other functions

The Stumpff functions can be expressed in terms of the Mittag-Leffler function:

$\ c_{k}(x)~=~E_{2,k+1}(-x)~~.$