# Solving Kepler's equation: discussion of theory, plus C code

Last updated 2013 April 19

In 1998, a post on the sci.astro.amateur newsgroup asked about a solution to Kepler's equation for the hyperbolic case. The following code is extracted from that used in my Guide and Find_Orb software, in hopes that it will prove useful to someone. You can find the current version in astfuncs.cpp in the basic astronomical functions library also on this site. I recommend use of the updated version; it includes a couple of fixes for high-eccentricity and near-parabolic cases.

There have been an amazing number of solutions offered to Kepler's equation over the years, most of which try to optimize the following factors:

Simplicity. The solution offered below is somewhat complex; for my particular purpose, that was entirely acceptable, since it allowed me to do quite well in some other areas. If I were teaching an introductory course in astronomical computing, I might start with something simpler. And, of course, complexity is the enemy of bug-free software.

Completeness. A lot of solutions converge for elliptical orbits only, or for low-eccentricity cases. (For example, the Kepler solver used in most artificial satellite software is in this category. You don't deal with hyperbolic or near-parabolic cases in that instance.) Again, for my particular purpose, all cases must converge to a solution to be acceptable. For example, during a planetary fly-by, it's possible to have eccentricities in the thousands.

Worst-case behavior. Some solutions converge quite rapidly in most cases, but may end up with really slow convergence (usually for nearly-parabolic orbits). In my case, having poor behavior in such cases would not be acceptable.

Average-case behavior. The definition of an "average" case will obviously vary! In my software, there are an enormous number of cases with eccentricity < .3 (asteroids), and a lot of cases that are nearly parabolic, but not quite (comets). So I naturally invested some thought in making sure those converge as quickly as possible.

Generalness of the solution. You can, without much trouble, write code that will find the solution to a generalized continuous function, given that you know a solution exists between two points. The benefit of this is that your "equation solver" can then be used for pretty much anything you want. In my case, this very particular equation was so important that I felt it worthwhile to sacrifice "generalness" to get a much faster solution.

For my purposes, it was necessary that the code converge for all cases, elliptical, hyperbolic, and parabolic; and that it do so as briskly as possible. (I've omitted the parabolic case below, by the way. In that case, an analytical, non-iterative solution is used.) In some cases, Guide must compute positions for tens of thousands of asteroids, and do so quickly enough to update the chart right away. Find_Orb may, in the course of numerically integrating an object's orbit, call this function millions of times.

Also, the first version of this code was created in 1991, when math chips were still scarce; I was willing to do almost anything to make the code run faster. I probably wouldn't go to such extremes in the modern era... but since I've already got the code, and it's tested and ready... why not?

And now, a blow-by-blow description of how the kepler( ) function works. The first step is just to check for a zero mean anomaly; in this case, the return value is trivially zero, and the function returns right away. (This also evades divide-by-zero cases later on.)

Next, if the eccentricity is less than .3, we use a "low-eccentricity" formula from Meeus' Astronomical Algorithms. This gets us a very good solution; follow up with one iteration (pardon the oxymoron!) of Newton's method, and we have our solution. This was important for that "zillions of asteroids" case. Almost all asteroids will stop here; few have high eccentricities.

For what follows, it's necessary to be on the positive half of the orbit. If the mean anomaly is negative, we flip the sign and record that fact by setting is_negative = 1.

Now, the step following this will require most of the explanation. In the end, the only thing we'll be doing is a simple iteration of Newton's method. But to make sure it converges, we'll need a good starting value. For elliptical orbits that either have an eccentricity less than .8 or a mean anomaly greater than 60 degrees, setting our initial guess for the eccentric anomaly equal to the mean anomaly is good enough. But if we've got a hyperbolic orbit, or if the eccentricity is greater than .8 and the mean anomaly is below 60 degrees, a smarter initial guess may be absolutely essential to guarantee convergence. (In other cases, we'd get convergence, but it would be embarrassingly slow.)

The basic idea is to perform a series expansion of Kepler's equation, and to keep only the first few terms. That equation can be solved or approximated analytically, and the result is a good enough solution to Kepler's equation that we can be sure of getting convergence. Here's the math, for both the elliptical and hyperbolic cases.

Elliptical case:

```M = E - ecc * sin( E)
M = E - ecc * (E - E^3 / 6 + higher terms)
M = (1 - ecc) * E + (ecc / 6) * E^3 + higher terms
```

Hyperbolic case:

```M = ecc * sinh( E) - E
M = ecc * (E + E^3 / 6 + higher terms) - E
M = (ecc - 1) * E + (ecc / 6) * E^3 + higher terms
```

(1 - ecc) is always positive in the elliptical case (ecc < 1) and (ecc - 1) is also always positive in the hyperbolic case (ecc > 1). So really, we have only one equation to consider:

```M = fabs( 1 - ecc) * E + (ecc / 6) * E^3
```

(Note to non-C programmers: fabs = floating-point absolute value.)

For a starter, we guess that E^3 is not going to be all that important. If so, then a good guess for E would be E = M / fabs( 1 - ecc). So we compute that value, and then test our guess that this first term will dominate the second term. If, instead, our guess for E leads to

```fabs( 1 - ecc) * E < (ecc / 6) * E^3
6 * fabs( 1 - ecc) / ecc < E^2
(assume ecc is close to 1)
6 * fabs( 1 - ecc) < E^2
```

then it's really the cubic, or the higher-order terms that dominate. So we switch instead to dropping the term in E, and we get

```M = (ecc / 6) * E^3
E = cube_root( 6 * M / ecc)
```

Now, in truth, if you have a hyperbolic case with (roughly) M > pi, those "higher-order" terms in the sinh expansion start to dominate. So in that case, we switch to:

```M = ecc * sinh( E) - E
(E is a lot smaller than ecc * sinh( E))
M = ecc * sinh( E)
E = inverse_sinh( M / ecc)
```

For each of these cases, we've got a starting value of E that is guaranteed to converge to a solution. The approximations made above are admittedly crude; there is always the tradeoff between "better starting approximation allowing fewer Newton steps, at the cost of more math up-front" versus "lousy approximation that can be done with little math, but that then requires seemingly endless iterations." The above is the result of creating a lot of contour plots, showing the number of iterations as a function of M and ecc. In particular, most comets have eccentricities very close to, but not equal to, 1, with small M during their time near perihelion; I regarded that as a crucial case, and invested much skull sweat in getting that to converge briskly.

An interesting point concerning the hyperbolic case has been made by Chris Marriott, author of SkyMap . He noted that the two "guesses" for an initial anomaly, E1=cube_root(6*M) and E2=asinh(M/ecc), define bounds for the actual value of E; you can be certain that E2 < E < E1. This makes implementing a secant or binary-search method much easier.

```#include <math.h>

#define PI 3.14159265358979323
#define THRESH 1.e-8
#define CUBE_ROOT( X)  (exp( log( X) / 3.))

static double asinh( const double z);
double kepler( const double ecc, double mean_anom);

static double asinh( const double z)
{
return( log( z + sqrt( z * z + 1.)));
}

static double kepler( const double ecc, double mean_anom)
{
double curr, err, thresh;
int is_negative = 0, n_iter = 0;

if( !mean_anom)
return( 0.);

if( ecc < .3)     /* low-eccentricity formula from Meeus,  p. 195 */
{
curr = atan2( sin( mean_anom), cos( mean_anom) - ecc);
/* one correction step,  and we're done */
err = curr - ecc * sin( curr) - mean_anom;
curr -= err / (1. - ecc * cos( curr));
return( curr);
}

if( mean_anom < 0.)
{
mean_anom = -mean_anom;
is_negative = 1;
}

curr = mean_anom;
thresh = THRESH * fabs( 1. - ecc);
if( ecc > .8 && mean_anom < PI / 3. || ecc > 1.)    /* up to 60 degrees */
{
double trial = mean_anom / fabs( 1. - ecc);

if( trial * trial > 6. * fabs(1. - ecc))   /* cubic term is dominant */
{
if( mean_anom < PI)
trial = CUBE_ROOT( 6. * mean_anom);
else        /* hyperbolic w/ 5th & higher-order terms predominant */
trial = asinh( mean_anom / ecc);
}
curr = trial;
}

if( ecc < 1.)
{
err = curr - ecc * sin( curr) - mean_anom;
while( fabs( err) > thresh)
{
n_iter++;
curr -= err / (1. - ecc * cos( curr));
err = curr - ecc * sin( curr) - mean_anom;
}
}
else
{
err = ecc * sinh( curr) - curr - mean_anom;
while( fabs( err) > thresh)
{
n_iter++;
curr -= err / (ecc * cosh( curr) - 1.);
err = ecc * sinh( curr) - curr - mean_anom;
}
}
return( is_negative ? -curr : curr);
}
```