Orbit Propagation

Listing 3

// This function solves Kepler's equation by the Newton
// method.
// Inputs:  m --  mean anomaly
//          e0 -- starting estimate for eccen anom
//          e  -- orbit eccentricity
//          mi -- maximum number of iterations
//          t  -- stopping tolerance
// Outputs: sk returns the final estimate for eccen anom
//   * Input m and e0 and output sk are in degrees. *

double sk(double m,double e0,double e, int mi,double t)
{
  /* Solve Kepler's equation by the Newton method. */
  double en,ep;
  int i;

  e0 = e0 * DTR;
  m  = m  * DTR;

  i  = 0;
  en = e0;
  do {
    i++;
    ep = en;

    // The next line of code is the Newton iteration:
    //   x(n+1) = x(n) - (1/(dy/dx)) * y(n).

    en = en - (en - e*sin(en) - m)/(1 - e*cos(en))  ;
    en = fmod(en,TWOPI);
    if (fabs(en-ep) < t) break;
  } while (i <= mi );
  return(en*RTD);
}





<B>Sample Input/Output</B>

Initial eccentric anom = 0
Input mean anom        = 286.879298
Input maximum its      = 15
Input stop tolerance   = 0.000001
Input eccentricity     = 0.100000
Final eccentric anom   = 281.260008
Mean anom computed from solution = 286.879298

This solution took five iterations; the solution E was
checked by plugging it back into Kepler's equation
and computing M.