Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# j-faria/kepler.cpp

Created Jan 16, 2017
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); }
to join this conversation on GitHub. Already have an account? Sign in to comment