I know it is a controversial opinion, but you might need to know where things are when running a simulation. As a bonus, it helps you in knowing what to draw on the screen.

Of course, you can just use a position vector (x, y, z) where x, y and z are the coordinates of the object. But things have the inconvenient tendency not to stay in the same place for all eternity. So we’ll want to know where that object is at the given time.

Now, the more generic approach is to run a physics simulation. We compute the acceleration of the object, update its speed depending on the acceleration and update its position depending on its speed. By the way, use Runge-Kutta 4 for this, not Euler’s method.

That would work for things in orbit too. But we can do much better in this case. This one guy, Johannes Kepler found simpler formulas to describe the movement of two orbiting objects (e.g. the Sun and the Earth, or the Earth and the Moon). Basically, the trajectories will be ellipses, and we can deduce the position of an object in its orbit at any given time.

We can describe the position of the object in its orbit as the angle between the direction of the periapsis (where the planet is closest to its sun) and the direction to the object, taken at the relevant focus. This is the true anomaly. Everything is an illusion, so we also like to pretend the ellipse is actually a circle and define the eccentric anomaly. But these angles are not linear with time. So we also define the mean anomaly, which works like an angle (goes from 0 to 2π) but is linear with time. Once we know the direction of the object, we will only need to know its distance, r.

And here are the formulas to get from time to the true anomaly, and how we then obtain the distance from focus (e is the eccentricity of the ellipse):

• mean anomaly (M): M = M₀ + n × (t – t₀)
• eccentric anomaly (E): M = E – e × sin E
• true anomaly (f): $latex \displaystyle f=2\arctan\sqrt{\frac{1+e}{1-e}} \tan\frac E2$
• distance from focus (r): $latex \displaystyle r=\frac p{1 + e \cos f}$

You won’t be calculating this in your head, but the expressions are pretty easy to translate to working code. Well, except for one. Notice how the eccentric anomaly is not expressed as a function of the mean anomaly? Well, that’s the hard part. And it’s called Kepler’s equation.

We do not know a simple way to express E as a formula depending only on e and M. Instead, we must proceed by successive approximations until we decide we are close enough. For this, we will need two things:

• a starting value E₀ for E
• a way to improve the current approximation of E

Most introductory books and online articles I have found will tell you to take E₀ = M. This has been obsolete since 1978. It works well enough, but fails in corner cases (near-parabolic orbits around periapsis). And modern starting values are close enough to limit the number of improvements you then need to make. I reviewed 21 papers, selected the 3 most relevant ones, and implemented 21 methods for picking starting values (3 papers for the hyperbolic case with 5 methods for the hyperbolic case).

To improve the starting value, old references will make you do E := M + e sin E, but this is pretty naive. Most modern references will use Newton’s method, which helps you solve an equation faster as long as you know its first derivative. But you can go further by using the second derivative using Halley’s method. You can go arbitrarily far using Householder’s method.

Great, now I have 21 × 4 = 84 approaches to try (5 × 3 = 15 for the hyperbolic case). Which one is the most accurate? Which one is the fastest? And, more importantly, can I find a fast and accurate approach?

You can spend years studying the theoretical properties of these functions. Or, we just try them for real. It is the most efficient way to know how they actually behave on real hardware with real implementations of numbers.

I wrote a benchmark tool to test both the accuracy and the speed of each approach. It turns out that only a handful of approaches are accurate enough for my taste:

# Elliptical case
## Naive method
## Newton's method
## Halley's method
     580 c  3.580e-16  gooding_10
     764 c  3.926e-16  gooding_11
     577 c  3.604e-16  mikkola_1
     742 c  3.653e-16  mikkola_2
## Householder's third order method
     606 c  3.645e-16  gooding_10
     809 c  3.623e-16  gooding_11
     641 c  3.935e-16  mikkola_1
     764 c  3.904e-16  mikkola_2

# Hyperbolic case
## Newton's method
## Halley's method
    1454 c  3.731e-15  mikkola_1
    1578 c  3.725e-15  mikkola_2
## Householder's third order method
    1704 c  3.618e-15  mikkola_1
    1757 c  3.618e-15  mikkola_2

The first column is the number of CPU cycles per solve of Kepler’s equation. The second column is the worst relative error over the test set. The third column is a short name for the method for picking a starting value for E.

As we can see above, no approach seems to work well enough with the naive improvement, or with Newton’s method. And using third derivatives does not do much. gooding_11 and mikkola_2 are relatively slow. This leaves gooding_10 and mikkola_1. I picked mikkola_1 because the paper also treated the hyperbolic case. It comes from A cubic approximation for Kepler’s equation by Seppo Mikkola.

And after all this, I have a very robust way to know where stuff is at any time! Of course, you should not forget to take care of floating point precision issues when implementing all of this…

Appendix: References

• smith_* are from A simple, efficient starting value for the iterative solution of Kepler’s equation (1979), by Gary R Smith
• mikkola_* are from A cubic approximation for Kepler’s equation (1987), by Seppo Mikkola
• gooding_* are from Procedures for Solving Kepler’s Equation (1986) for the elliptic case, and from The hyperbolic Kepler equation (and the elliptic equation revisited) (1988) for the hyperbolic case, by A. W. Odell and R. H. Gooding
• pfelger is from Astronomy on the Personal Computer (2000) by Oliver Montenbruck and Thomas Pfleger

Appendix: Full benchmark

# Elliptical case
## Naive method
     264 c  2.500e-01  smith_1
     276 c  5.303e-01  smith_2
     360 c  2.000e-01  smith_3
     444 c  9.701e-02  smith_4
     385 c  1.667e-01  smith_5
     480 c  6.495e+01  smith_6
     246 c  2.500e-01  gooding_1
     299 c  2.000e-01  gooding_2
     356 c  1.667e-01  gooding_3
     246 c  5.303e-01  gooding_4
     460 c  9.701e-02  gooding_5
     253 c  2.417e-01  gooding_6
     325 c  7.497e-02  gooding_7
     334 c  7.497e-02  gooding_7b
     378 c  1.277e-01  gooding_8
     423 c  7.499e-02  gooding_9
     455 c  2.376e-01  gooding_10
     608 c  8.020e-02  gooding_11
     286 c  1.439e-01  gooding_12
     448 c  4.342e-01  mikkola_1
     629 c  4.426e-01  mikkola_2
## Newton's method
     522 c  8.517e-03  smith_1
     539 c  1.379e-02  smith_2
     578 c  8.111e-03  smith_3
     697 c  8.111e-03  smith_4
     706 c  8.111e-03  smith_5
     900 c  1.028e+02  smith_6
     481 c  8.517e-03  gooding_1
     583 c  8.111e-03  gooding_2
     561 c  8.111e-03  gooding_3
     476 c  1.379e-02  gooding_4
     686 c  8.111e-03  gooding_5
     537 c  3.528e-02  gooding_6
     536 c  1.378e-02  gooding_7
     585 c  1.378e-02  gooding_7b
     702 c  2.524e-04  gooding_8
     680 c  1.379e-02  gooding_9
     731 c  1.117e-10  gooding_10
     940 c  2.072e-08  gooding_11
     528 c  5.673e-03  gooding_12
     697 c  1.825e-06  mikkola_1
     815 c  1.825e-06  mikkola_2
## Halley's method
     395 c  1.250e-01  smith_1
     419 c  2.553e-03  smith_2
     488 c  6.250e-02  smith_3
     596 c  1.710e-02  smith_4
     557 c  4.167e-02  smith_5
     989 c  5.130e+01  smith_6
     429 c  1.250e-01  gooding_1
     503 c  6.250e-02  gooding_2
     524 c  4.167e-02  gooding_3
     426 c  2.553e-03  gooding_4
     587 c  1.710e-02  gooding_5
     410 c  6.127e-03  gooding_6
     467 c  2.548e-03  gooding_7
     440 c  2.548e-03  gooding_7b
     508 c  9.153e-05  gooding_8
     552 c  2.550e-03  gooding_9
     598 c  3.580e-16  gooding_10
     821 c  3.926e-16  gooding_11
     430 c  1.028e-03  gooding_12
     594 c  3.604e-16  mikkola_1
     742 c  3.653e-16  mikkola_2
## Householder's third order method
     438 c  1.562e-02  smith_1
     439 c  6.779e-04  smith_2
     501 c  7.812e-03  smith_3
     667 c  2.138e-03  smith_4
     533 c  5.208e-03  smith_5
    1053 c  6.535e+01  smith_6
     453 c  1.562e-02  gooding_1
     501 c  7.812e-03  gooding_2
     594 c  5.208e-03  gooding_3
     429 c  6.779e-04  gooding_4
     595 c  2.138e-03  gooding_5
     454 c  1.638e-03  gooding_6
     526 c  6.736e-04  gooding_7
     585 c  6.736e-04  gooding_7b
     583 c  9.493e-09  gooding_8
     652 c  6.741e-04  gooding_9
     681 c  3.645e-16  gooding_10
     846 c  3.623e-16  gooding_11
     486 c  2.724e-04  gooding_12
     651 c  3.935e-16  mikkola_1
     752 c  3.904e-16  mikkola_2

# Hyperbolic case
## Newton's method
     780 c  7.245e+25  basic
    1064 c  1.421e-03  mikkola_1
    1201 c  1.970e-04  mikkola_2
    1014 c  8.559e-01  gooding
     892 c  1.515e+19  pfelger
## Halley's method
    1134 c  1.729e+17  basic
    1396 c  3.731e-15  mikkola_1
    1526 c  3.725e-15  mikkola_2
    1265 c  9.160e+10  gooding
    1237 c  1.791e+13  pfelger
## Householder's third order method
    1260 c  1.323e+07  basic
    1631 c  3.618e-15  mikkola_1
    1595 c  3.618e-15  mikkola_2
    1467 c  2.276e+17  gooding
    1365 c  9.997e+01  pfelger