Skip to content

Instantly share code, notes, and snippets.

@j-faria j-faria/kepler.cpp
Created Jan 16, 2017

Embed
What would you like to do?
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
You can’t perform that action at this time.