the particle’s position in time, we make use of several special functions. Let us begin by defining the $c$-functions [@Stumpff1962] as a series expansion: $$\begin{aligned}
c_n(z) \equiv \sum_{j=0}^\infty \frac{(-z)^j}{(n+2j)!} \label{eq:c},\end{aligned}$$ which satisfy the recursion relation $$\begin{aligned}
c_n(z) &=& \frac{1}{n!} - z \,c_{n+2}. \label{eq:recur}\end{aligned}$$ The $c$-functions are related to trigonometric functions, for example $$\begin{aligned}
c_0(z) = \cos \sqrt{z} \quad\quad \text{and}\quad\quad c_1(z) = \frac{\sin\sqrt{z}}{\sqrt{z}},\end{aligned}$$ and thus satisfy the following relationships [@Mikkola1997], which are related to the half-angle formula for trigonometric functions: $$\begin{aligned}
c_5(z) &=& \frac1{16} \left[ c_5(z/4) + c_4(z/4) + c_3(z/4)c_2(z/4)\label{eq:crel5}\right]\\
c_4(z) &=& \frac1{8}  c_3(z/4) \left[1 + c_1(z/4).\label{eq:crel4}\right]\end{aligned}$$ Values for $c_0$ through $c_3$ are then readily computed from Eq. \[eq:recur\]. Next, we introduce the so called $G$-functions [@StiefelScheifele1971] which in turn depend on the $c$-functions: $$\begin{aligned}
G_n(\beta,X) \equiv X^n c_n(\beta X^2). \label{eq:G}\end{aligned}$$ The $G$-functions also satisfy recursion relationships similar to those mentioned above for the $c$-functions. We can easily calculate derivatives of $G_n$ by looking at the series expansion of $c_n$ [see @MikkolaInnanen1999 for details]. With this framework, we can now write down the steps needed to find the solution to the Kepler Hamiltonian in compact form.

First, we need to calculate the following three quantities from the initial conditions $\mathbf{r}_0, \mathbf{v}_0$: 