diff --git a/src/orbital_mechanics.cpp b/src/orbital_mechanics.cpp index 10798ad..880fa5c 100644 --- a/src/orbital_mechanics.cpp +++ b/src/orbital_mechanics.cpp @@ -66,20 +66,21 @@ void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass, *out_velocity = mat3_multiply_vec3(rotation, velocity); } +// Shared solver constants +static const double KEPLER_TOLERANCE = 1.0e-10; +static const int KEPLER_MAX_ITERATIONS = 50; + double get_initial_trial_value(double mean_anomaly, double eccentricity) { return mean_anomaly + eccentricity * sin(mean_anomaly) + ((pow(eccentricity, 2) / 2.0) * sin(2.0 * mean_anomaly)); } -double solve_kepler_equation(double mean_anomaly, double eccentricity) { - const double CONVERGENCE_TOLERANCE = 1.0e-10; - const int MAX_ITERATIONS = 50; - +double solve_kepler_elliptical(double mean_anomaly, double eccentricity) { double E = get_initial_trial_value(mean_anomaly, eccentricity); - double E_prev = E + 2.0 * CONVERGENCE_TOLERANCE; + double E_prev = E + 2.0 * KEPLER_TOLERANCE; int iterations = 0; - while (fabs(E - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) { + while (fabs(E - E_prev) > KEPLER_TOLERANCE && iterations < KEPLER_MAX_ITERATIONS) { E_prev = E; double sin_E = sin(E); E = E - (E - eccentricity * sin_E - mean_anomaly) / (1.0 - eccentricity * cos(E)); @@ -89,6 +90,52 @@ double solve_kepler_equation(double mean_anomaly, double eccentricity) { return E; } +double solve_kepler_hyperbolic(double mean_anomaly, double eccentricity) { + // Initial guess for hyperbolic anomaly + double H = mean_anomaly; + if (eccentricity * sinh(mean_anomaly) > mean_anomaly) { + H = log(2.0 * mean_anomaly / eccentricity); + } + + double H_prev = H + 2.0 * KEPLER_TOLERANCE; + int iterations = 0; + + while (fabs(H - H_prev) > KEPLER_TOLERANCE && iterations < KEPLER_MAX_ITERATIONS) { + H_prev = H; + double sinh_H = sinh(H); + double cosh_H = cosh(H); + H = H - (H - eccentricity * sinh_H - mean_anomaly) / (1.0 - eccentricity * cosh_H); + iterations++; + } + + return H; +} + +double eccentric_to_true_anomaly(double eccentric_anomaly, double eccentricity) { + // E to true anomaly conversion: tan(ν/2) = √((1+e)/(1-e)) · tan(E/2) + double tan_half_E = tan(eccentric_anomaly / 2.0); + double tan_half_nu = sqrt((1.0 + eccentricity) / (1.0 - eccentricity)) * tan_half_E; + return 2.0 * atan(tan_half_nu); +} + +double hyperbolic_to_true_anomaly(double hyperbolic_anomaly, double eccentricity) { + // Hyperbolic E to true anomaly: tanh(ν/2) = √((e-1)/(e+1)) · tanh(H/2) + double tanh_half_H = tanh(hyperbolic_anomaly / 2.0); + double tanh_half_nu = sqrt((eccentricity - 1.0) / (eccentricity + 1.0)) * tanh_half_H; + return 2.0 * atanh(tanh_half_nu); +} + +double mean_anomaly_to_true_anomaly(double mean_anomaly, double eccentricity) { + if (eccentricity < 1.0) { + double E = solve_kepler_elliptical(mean_anomaly, eccentricity); + return eccentric_to_true_anomaly(E, eccentricity); + } else { + double H = solve_kepler_hyperbolic(mean_anomaly, eccentricity); + return hyperbolic_to_true_anomaly(H, eccentricity); + } +} + +// TODO: refactor for readability OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, double parent_mass) { double mu = G * parent_mass; diff --git a/src/orbital_mechanics.h b/src/orbital_mechanics.h index 9ea9bca..f3ffbf9 100644 --- a/src/orbital_mechanics.h +++ b/src/orbital_mechanics.h @@ -20,10 +20,23 @@ void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass, OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, double parent_mass); -double solve_kepler_equation(double mean_anomaly, double eccentricity); - +// Initial guess for Newton-Raphson: M + e·sin(M) + (e²/2)·sin(2M) double get_initial_trial_value(double mean_anomaly, double eccentricity); +// Elliptical Kepler equation solver: E - e·sin(E) = M +double solve_kepler_elliptical(double mean_anomaly, double eccentricity); + +// Hyperbolic Kepler equation solver: H - e·sinh(H) = M +double solve_kepler_hyperbolic(double mean_anomaly, double eccentricity); + +// Conversions between anomaly types +double eccentric_to_true_anomaly(double eccentric_anomaly, double eccentricity); +double hyperbolic_to_true_anomaly(double hyperbolic_anomaly, double eccentricity); + +// Unified mean anomaly to true anomaly conversion +// Automatically dispatches to elliptical or hyperbolic based on eccentricity +double mean_anomaly_to_true_anomaly(double mean_anomaly, double eccentricity); + OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, double dt, double parent_mass); #endif diff --git a/tests/test_newton_raphson_convergence.cpp b/tests/test_newton_raphson_convergence.cpp index e0e4eaf..80b6b4c 100644 --- a/tests/test_newton_raphson_convergence.cpp +++ b/tests/test_newton_raphson_convergence.cpp @@ -18,15 +18,24 @@ TEST_CASE("Newton-Raphson solver - very low eccentricity (e < 0.01)", "[newton][ double mean_anomaly = M_PI / 2.0; - double eccentric_anomaly = solve_kepler_equation(mean_anomaly, e); - double expected_eccentric_anomaly = mean_anomaly; + double eccentric_anomaly = solve_kepler_elliptical(mean_anomaly, e); - double error = fabs(eccentric_anomaly - expected_eccentric_anomaly); - INFO("Eccentric anomaly: " << eccentric_anomaly << " rad"); - INFO("Expected: " << expected_eccentric_anomaly << " rad"); - INFO("Error: " << error); + // Verify Kepler's equation is satisfied: E - e*sin(E) = M + // Note: E ≈ M only when e=0 exactly. For small e, E ≈ M + e*sin(M) + double kepler_residual = eccentric_anomaly - e * sin(eccentric_anomaly) - mean_anomaly; + double first_order_approx = mean_anomaly + e * sin(mean_anomaly); + double approximation_error = fabs(eccentric_anomaly - first_order_approx); - REQUIRE(error < 1.0e-6); + INFO("Eccentric anomaly: " << eccentric_anomaly << " rad"); + INFO("Mean anomaly: " << mean_anomaly << " rad"); + INFO("Kepler's equation residual |E - e·sin(E) - M|: " << kepler_residual); + INFO("First-order approx E ≈ M + e·sin(M): " << first_order_approx << " rad"); + INFO("|E - approx|: " << approximation_error); + + // Verify solver correctly solves Kepler's equation + REQUIRE(fabs(kepler_residual) < CONVERGENCE_TOLERANCE); + // Verify result matches first-order approximation (valid for small e) + REQUIRE(approximation_error < 0.01); } } @@ -39,7 +48,7 @@ TEST_CASE("Newton-Raphson solver - moderate eccentricity (0.1 < e < 0.5)", "[new double mean_anomaly = M_PI / 4.0; - double eccentric_anomaly = solve_kepler_equation(mean_anomaly, e); + double eccentric_anomaly = solve_kepler_elliptical(mean_anomaly, e); double rhs = mean_anomaly + e * sin(eccentric_anomaly); double residual = eccentric_anomaly - rhs;