Solving Kepler’s Equation

Warning: The routines below, that has B2 perturbations have a bug. For some values the orbit is reflected through the origin. I will put a fix here as soon as I find one. For the time being, please do not use drift_shep_pert() routines!

Warning2: This is simply a copy of old documentation. It may not be up to date if I make changes to the code (note added on 2017-06-30).

Solving Kepler’s equation for arbitrary epoch and eccentricity is a common problem in celestial mechanics. Below, I present my implementation of such a solver, written in C. It is an extension of Shepperd’s (1985) work, and uses ideas from Danby (1992) and Mikkola (1999).

Given the initial position and velocity, and the time interval; the code calculates the new position and the velocity of a particle for the Hamiltonian

H = -p2 - μ-- B2-.
    2m    r   r2
(1)

Why a new routine

There are of course codes already available for this purpose; perhaps most notably the SWIFT routines. They are heavily tested and very reliable. However, I wanted to write my own routines, for a number of reasons.

Please note that my experience with SWIFT routines is very limited. It is possible that one can make a few modifications to increase their accuracy while maintaining their reliability. I’d appreciate any feedback on this matter.

Solving Kepler’s equation

For simplicity let us initially assume B2 = 0. Then the orbit is a Keplerian conic section. We start by calculating the quantities

                                    2μ
r0 = |r0|, v0 = |v0|, σ0 = r0⋅v0,  β = --- v20,
                                    r0
(2)

where r0 and v0 are the initial position and velocity vectors, respectively. If β > 0, we have an elliptical orbit and can calculate the period

         -3∕2
P = 2πμβ     ,
(3)

and the quantity

                   ⌊               ⌋
         -5∕2       δt-+P-∕2--2σ0∕β
δU = n2πβ    , n =         P         ,
(4)

where δt is the given time interval. This will be useful if δt > P, a situation that I do not encounter since I always choose timesteps much shorter than the period.

For our independent variable, we choose the initial value u = 0. The main loop consists of the following steps (prime denotes differentiation with respect to u):

    -βu2---
q = 1+ βu2 ,
(5)

 ′  ---2βu---
q = (1+ βu2)2 ,
(6)

 ′′  ----2β---   --8β2u2--
q = (1+ βu2)2 - (1+ βu2)3 ,
(7)

U˜0 = 1- 2q,
(8)

U˜1 = 2(1 - q)u ,
(9)

    16  5
U = 15U˜1 G5 (q)+ δU ,
(10)

U0 = 2 ˜U02 - 1,
(11)

U1 = 2 ˜U0U˜1 ,
(12)

      ˜ 2
U2 = 2U1 ,
(13)

U3 = βU + U1U2∕3,
(14)

r = r0U0 + σ0U1 + μU2,
(15)

 ′
r = 4(1- q)(σ0U0 + (μ - βr0)U1),
(16)

r′′ = - 4q′(σ0U0 + (μ- βr0)U1)+ 16(1- q)2(- βσ0U1 + (μ - βr0)U0),
(17)

Δt = r0U1 + σ0U2 + μU3 ,
(18)

f = Δt - δt,
(19)

f′ = 4(1 - q)r,
(20)

f′′ = 4(r′(1 - q) - rq′),
(21)

 ′′′     ′′     ′′         ′′
f  = - 8r q- 4rq + 4(1- q)r ,
(22)

δu1 = - f∕f′,
(23)

δu2 = - f ∕(f′ + δu1f′′∕2),
(24)

           ′     ′′      2 ′′′
δu3 = - f∕(f +δu2f ∕2 + δu2f ∕6),
(25)

u = u+ δu3.
(26)

Once the increment δu3 or f is smaller than a preset tolerance, we exit the loop. The function G5(x) is a special case of Gauss’s hypergeometric function: G5(x) = 2F1(5,1;72;x) (Abramowitz and Stegun1972, Ch. 15). The position and velocities are calculated by

f = 1- (μ∕r0)U2 ,
(27)

g = r0U1 + σ0U2,
(28)

F = - μU1∕(rr0),
(29)

G = 1 - (μ∕r)U2 ,
(30)

r = fr0 + gv0,
(31)

v = Fr0 + Gv0 .
(32)

This method of solution is almost identical to Shepperd’s. The only thing I did is to increase the order of the iteration by using Danby’s technique. In the actual implementation, if q > 12 during iteration, or if convergence is not achieved after 12 iterations, we stop, and try to cover δt in two steps of δt∕2.

When B20, we need two modifications. In the following, I adopt the approach of Mikkola (1999); see that paper for generating functions etc. First note that in the plane of motion, we can transform to polar coordinates (r,θ) and write the Hamiltonian as

      (      2)
H =  1  p2r + pθ2 - μ-- B22-
     2(     r     r ) r
  =  1  p2+ p2θ --2B2 -  μ-
     2   r     r2       r
     1(     p2)   μ
  =  -  p2r +-ψ2- - --,
     2      r      r
(33)

which is in Keplerian form. The first modification is

β = 2μ-- v2+ 2B2-.
    r0   0    r20
(34)

Furthermore, we cannot use f-g formulation directly, so at the end of the loop, we revert to a more general method, which is applicable for any motion with a radial force:

pθ = r0 × v0,
(35)

p2ψ = p2θ - 2B2 ,
(36)

f = 1- (μ∕r0)U2 ,
(37)

g = r0U1 + σ0U2,
(38)

    ---g-----
ξ = r20f + σ0g ,
(39)

yψ = ξat1(p2ψξ2),
(40)

˙r = (σ0U0 + (μ - βr0)U1)∕r,
(41)

 2   2 2
θ = pθyψ ,
(42)

r =-r(c0(θ2)r0 + yψc1(θ2)pθ × r0),
   r0
(43)

v = r˙r + 12pθ × r.
    r    r
(44)

In this formulation c0(x2) and c1(x2) are Stumpff functions, and at1(x2) is arctan(x)∕x. All these functions, and the function G5(x) mentioned earlier, have simple continued fraction expansions. However, in my implementation, I only calculate G5(x) by continued fraction expansion to keep things simple. An efficient way to do this is given by Shepperd (1985). For calculating c0 and c1 by continued fractions, probably the best way is to use the approach of Flanders and Frame (1987). The continued fraction for at1(x) goes back to Lambert (1770) and Lagrange (1776) according to Olds (1963, Appendix II, formula 14):

              √ --
at1(x) = arcta√n(--x)= -----------1---------- .
             x      1 + -------1-⋅x--------
                           ------4-⋅x-----
                        3+         9⋅x
                           5 + -----16-⋅x---
                               7+ --------
                                  9 + 25⋅x-
                                       ...
(45)

The code

You can download the code from GitHub. The repository contains a few files. The most useful one is CF_kep_drift.c, which includes two routines. void drift_shep(double mu, double *r0, double *v0, double dt) solves Kepler’s equation for a particle with initial position r0, initial velocity v0, moving around a mass mu, over a time interval dt. The final positions and velocities are returned in r0 and v0. void drift_shep_pert(double mu, double B2, double *r0, double *v0, double dt) does the same but with perturbation B2. There is also a non-recursive version: drift_shep_nonrecurs(). This may be useful if you want to port the code to GPUs.

Kep_Drift_test.c is a small piece of code that I wrote to test my routines. It also has hooks for SWIFT routines, but I do not distribute them. If you download SWIFT, you can copy over the necessary files (listed in the Makefile), uncomment a few lines in the Makefile, change the definition of the macro DRIFT(mu, r, v) and compare the routines.

If there is some interest, I may post my comparison and other tests here.

Terms for using

This software is provided to help you get the job done. Hopefully to solve a problem, but possibly as a nice pastime. There are no restrictions whatsoever on the usage of this software. You can modify and adapt it to any use you see fit. However, I cannot be held responsible for any damage arising from using this software, including but not limited to wasting your time. Also, distributing the software, even parts of it as part of another software, is not usage. Such actions are governed by "Terms for copying and distribution", as indicated below.

I appreciate, but do not require, attribution. If you are unsure about what reference to use for attribution, feel free to ask. Note that, if you simply use this software to obtain some results to write a paper, you are not required to offer me co-authorship or such. Only if I provide extensive help to you for using the software, and hence I contribute significant amount of my time specifically for your project, I would like to be offered co-authorship. Do not let this to discourage you from asking questions, but please read the documentation (at least this page and possibly the references) before doing so.

In particular, feel free to port this code to GPUs. I also plan to do this at some point, but I think the code is already useful on its own.

Terms for copying and distribution

This software is copyright (c) 2008 by Atakan G├╝rkan, except where noted in the source files.

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License version 3 as published by the Free Software Foundation.

Note that I leave out the phrase "... or any later version", that is, you cannot merge this software with software under different versions of GPL. If you need to do this, please contact me and if the version you propose is reasonable for me, I will make an addition or exception to these terms.

This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

Parts of this software may be provided with licenses different from but compatible with GNU General Public License version 3. All such parts are clearly marked in the source code files. For those parts, you are free to use GNU General Public License version 3 or the particular license indicated in the source code.

Improvements and Bugfixes

If you make an improvement to this software, I would appreciate if you let me know about it. I also hope that you would allow me to further distribute your enhancements as part of this software, under the terms of this license. Note that, you do not need to transfer your copyright to me to make your improvements part of this software.

If you find an error in the formulae above, I really would like to hear about it. If there is an inconsistency between this text and the code, rely on the code. That is what I tested.

References

   M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions. Dover, New York, 1972.

   Richard H. Battin. An Introduction to the Mathematics and Methods of Astrodynamics. American Institute of Aeronautics and Astronautics, 1987.

   J. M. A. Danby. Fundamentals of celestial mechanics, 2nd ed. Willman-Bell, Richmond, 1992.

   Harley Flanders and J. Sutherland Frame. Algorithm of the bi-month: Elementary transcendental functions. The College Mathematics Journal, 18(5):417–421, 1987. ISSN 07468342. URL http://www.jstor.org/stable/2686970.

   R. W. Gosper. Continued Fractions. MIT HAKMEM Artificial Intelligence Memo, no 239: item 101, 1972. http://www.inwap.com/pdp10/hbaker/hakmem/cf.html and http://keithbriggs.info/cfup.htm.

   S. Mikkola. Efficient Symplectic Integration of Satellite Orbits. Celestial Mechanics and Dynamical Astronomy, 74:275–285, August 1999. doi: 10.1023/A:1008398121638.

   Carl Douglas Olds. Continued Fractions. Random House, 1963.

   S. W. Shepperd. Universal Keplerian state transition matrix. Celestial Mechanics, 35:129–144, February 1985.

You can download a better typeset version of first few sections here (LaTeX source).