C++ implementation of "A Practical Method for Solving the Kepler Equation"

kepler.cpp

/**
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);
}