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

Last updated 2016 September 15

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. The latest version can always be found on GitHub. I recommend using that version. It's happened that I made changes to the code that this page didn't quite keep up with.

The first thing I should mention is that solving Kepler's equation is almost trivial if you're only dealing with elliptical orbits. There are about ten lines of code below starting with if (ecc < .9) which pretty much do everything you'd need if your eccentricity is below that level. (Aside from some "housekeeping" problems of making sure the mean anomaly is between -180 and +180 degrees.) You can actually raise that limit to 0.9999 and use just that part of the code. It'll be slightly less efficient than the full code, but not by much, and will suffice for almost all work with asteroids.

But for comets and spacecraft, you have to consider near-parabolic, parabolic, slightly hyperbolic, and extremely hyperbolic orbits. Here, things become more difficult. In the course of reviewing the code just now, eighteen years after writing the first version, I found an unnecessary inefficiency and that the threshhold was poorly set for extremely hyperbolic orbits. There are all sorts of ways to go wrong. In hindsight, I might have done well to follow the advice of those who suggested using a universal variable formulation instead. (At least for eccentricities greater than .9.) UVF would add some overhead, but results in one code path instead of three cases. I think it would also remove the loss of precision issues one encounters for nearly parabolic orbits.

So here's a blow-by-blow description of how the kepler( ) function below 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 avoids divide-by-zero cases later on.

For elliptical cases, convergence is faster for -180 < mean_anom < 180. We do store the number of full revolutions added/subtracted, and subtract/add it to the final result. Thus, if you feed the code mean_anom = 1457 degrees and a low eccentricity, the return value will be somewhere near 1457 degrees, rather than near 1457 - 360*4 = 17 degrees. (Note that the actual math is all done in radians.)

Next, if the eccentricity is less than 0.9, we use a "low-eccentricity" formula from Meeus' Astronomical Algorithms. This gets us a very good solution; we follow up with Newton-Raphson iterations to get our solution. For an eccentricity of 0.2, three iterations will suffice, which covers the vast majority of asteroids and satellites. For e=0.9, six iterations will sometimes be required. (As notetd above, this is actually so quick and easy that I'd probably stick with it rather than using universal variables for these e<0.9 cases.)

Also as mentioned above, one can extend use of this loop for eccentricities up to about 0.9999. At that point, you may sometimes do eleven iterations.

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 = true.

We set a lower threshhold of convergence tolerance for near-parabolic orbits, where even a tiny change in eccentric anomaly can result in a large change of position.

Now, the step following this will require most of the explanation. In the end, we'll be doing the same Newton-Raphson iterations as we did for the e < 0.9 case. But we have to be a little more careful in picking the starting value for the iterations. For elliptical orbits with mean anomalies greater than 60 degrees, just setting curr = mean_anom will work well. But if we've got an hyperbolic orbit, or if the 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   ... multiply both sides by 6 / E :
6 * fabs( 1 - ecc) / ecc < E^2         ...now 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.

There is one remaining complication regarding nearly-parabolic orbits (which are common among comets, so it's an important case). As the eccentricity approaches 1, computing M=E-ecc*sin(E) can involve the subtraction of two almost identical quantities, a recipe for loss of precision. The fix is described in the near_parabolic function below.

#include <math.h>
#include <stdbool.h>

/* If the eccentricity is very close to parabolic,  and the eccentric
anomaly is quite low,  you can get an unfortunate situation where
roundoff error keeps you from converging.  Consider the just-barely-
elliptical case,  where in Kepler's equation,

M = E - e sin( E)

   E and e sin( E) can be almost identical quantities.  To
around this,  near_parabolic( ) computes E - e sin( E) by expanding
the sine function as a power series:

E - e sin( E) = E - e( E - E^3/3! + E^5/5! - ...)
= (1-e)E + e( -E^3/3! + E^5/5! - ...)

   It's a little bit expensive to do this,  and you only need do it
quite rarely.  (I only encountered the problem because I had orbits
that were supposed to be 'pure parabolic',  but due to roundoff,
they had e = 1+/- epsilon,  with epsilon _very_ small.)  So 'near_parabolic'
is only called if we've gone seven iterations without converging. */

static double near_parabolic( const double ecc_anom, const double e)
{
   const double anom2 = (e > 1. ? ecc_anom * ecc_anom : -ecc_anom * ecc_anom);
   double term = e * anom2 * ecc_anom / 6.;
   double rval = (1. - e) * ecc_anom - term;
   unsigned n = 4;

   while( fabs( term) > 1e-15)
      {
      term *= anom2 / (double)(n * (n + 1));
      rval -= term;
      n += 2;
      }
   return( rval);
}

#define MAX_ITERATIONS 7
#define THRESH 1.e-12
#define MIN_THRESH 1.e-15
#define CUBE_ROOT( X)  (exp( log( X) / 3.))
#define PI 3.1415926535897932384626433832795028841971693993751058209749445923

static double kepler( const double ecc, double mean_anom)
{
   double curr, err, thresh, offset = 0.;
   double delta_curr = 1.;
   bool is_negative = false;
   unsigned n_iter = 0;

   if( !mean_anom)
      return( 0.);

   if( ecc < 1.)
      {
      if( mean_anom < -PI || mean_anom > PI)
         {
         double tmod = fmod( mean_anom, PI * 2.);

         if( tmod > PI)             /* bring mean anom within -pi to +pi */
            tmod -= 2. * PI;
         else if( tmod < -PI)
            tmod += 2. * PI;
         offset = mean_anom - tmod;
         mean_anom = tmod;
         }

      if( ecc < .99999)     /* low-eccentricity formula from Meeus,  p. 195 */
         {
         curr = atan2( sin( mean_anom), cos( mean_anom) - ecc);
         do
            {
            err = (curr - ecc * sin( curr) - mean_anom) / (1. - ecc * cos( curr));
            curr -= err;
            }
            while( fabs( err) > THRESH);
         return( curr + offset);
         }
      }


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

   curr = mean_anom;
   thresh = THRESH * fabs( 1. - ecc);
               /* Due to roundoff error,  there's no way we can hope to */
               /* get below a certain minimum threshhold anyway:        */
   if( thresh < MIN_THRESH)
      thresh = MIN_THRESH;
   if( thresh > THRESH)       /* i.e.,  ecc > 2. */
      thresh = THRESH;
   if( 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.)
      while( fabs( delta_curr) > thresh)
         {
         if( n_iter++ > MAX_ITERATIONS)
            err = near_parabolic( curr, ecc) - mean_anom;
         else
            err = curr - ecc * sin( curr) - mean_anom;
         delta_curr = -err / (1. - ecc * cos( curr));
         curr += delta_curr;
         }
   else
      while( fabs( delta_curr) > thresh)
         {
         if( n_iter++ > MAX_ITERATIONS)
            err = -near_parabolic( curr, ecc) - mean_anom;
         else
            err = ecc * sinh( curr) - curr - mean_anom;
         delta_curr = -err / (ecc * cosh( curr) - 1.);
         curr += delta_curr;
         }
   return( is_negative ? offset - curr : offset + curr);
}