/// Kepler's equation: M = E - e * sin(E)
/// M is the Mean Anomaly (angle to where body would be if its orbit was actually circular)
/// E is the Eccentric Anomaly (angle to where the body is on the ellipse)
/// e is the eccentricity of the orbit (0 = perfect circle, and up to 1 is increasingly elliptical)
pub fn kepler_equation(eccentric_anomaly: f64, mean_anomaly: f64, eccentricity: f64) -> f64 {
    eccentricity.mul_add(eccentric_anomaly.sin(), mean_anomaly - eccentric_anomaly)
}