C++ implementation of "A Practical Method for Solving the Kepler Equation"
/** | |
Calculates the eccentric anomaly at time t by solving Kepler's equation. | |
See "A Practical Method for Solving the Kepler Equation", Marc A. Murison, 2006 | |
@param t the time at which to calculate the eccentric anomaly. | |
@param period the orbital period of the planet | |
@param ecc the eccentricity of the orbit | |
@param t_peri time of periastron passage | |
@return eccentric anomaly. | |
*/ | |
double ecc_anomaly(double t, double period, double ecc, double time_peri) | |
{ | |
double tol; | |
if (ecc < 0.8) tol = 1e-14; | |
else tol = 1e-13; | |
double n = 2.*M_PI/period; // mean motion | |
double M = n*(t - time_peri); // mean anomaly | |
double Mnorm = fmod(M, 2.*M_PI); | |
double E0 = keplerstart3(ecc, Mnorm); | |
double dE = tol + 1; | |
double E; | |
int count = 0; | |
while (dE > tol) | |
{ | |
E = E0 - eps3(ecc, Mnorm, E0); | |
dE = abs(E-E0); | |
E0 = E; | |
count++; | |
// failed to converge, this only happens for nearly parabolic orbits | |
if (count == 100) break; | |
} | |
return E; | |
} | |
/** | |
Provides a starting value to solve Kepler's equation. | |
See "A Practical Method for Solving the Kepler Equation", Marc A. Murison, 2006 | |
@param e the eccentricity of the orbit | |
@param M mean anomaly (in radians) | |
@return starting value for the eccentric anomaly. | |
*/ | |
double keplerstart3(double e, double M) | |
{ | |
double t34 = e*e; | |
double t35 = e*t34; | |
double t33 = cos(M); | |
return M + (-0.5*t35 + e + (t34 + 1.5*t33*t35)*t33)*sin(M); | |
} | |
/** | |
An iteration (correction) method to solve Kepler's equation. | |
See "A Practical Method for Solving the Kepler Equation", Marc A. Murison, 2006 | |
@param e the eccentricity of the orbit | |
@param M mean anomaly (in radians) | |
@param x starting value for the eccentric anomaly | |
@return corrected value for the eccentric anomaly | |
*/ | |
double eps3(double e, double M, double x) | |
{ | |
double t1 = cos(x); | |
double t2 = -1 + e*t1; | |
double t3 = sin(x); | |
double t4 = e*t3; | |
double t5 = -x + t4 + M; | |
double t6 = t5/(0.5*t5*t4/t2+t2); | |
return t5/((0.5*t3 - 1/6*t1*t6)*e*t6+t2); | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment