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

Last updated 12 Feb 98

A recent 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 software, in hopes that it will prove useful to someone. (Please don't use this in commercial software without notifying me. We might work out a trade or something, but I'm not too interested in supporting my competition to this extent! But for use in your own personal projects... use it freely.)

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:

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. 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 Pentium 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);
}