Refactor orbital_mechanics: separate elliptical and hyperbolic Kepler solvers

New modular API:
- solve_kepler_elliptical(M, e): Newton-Raphson for E - e·sin(E) = M
- solve_kepler_hyperbolic(M, e): Solver for H - e·sinh(H) = M
- eccentric_to_true_anomaly(E, e): Convert eccentric to true anomaly
- hyperbolic_to_true_anomaly(H, e): Convert hyperbolic to true anomaly
- mean_anomaly_to_true_anomaly(M, e): Unified wrapper (dispatches based on e)

Changes:
- Renamed solve_kepler_equation() → solve_kepler_elliptical() for clarity
- Extracted KEPLER_TOLERANCE and KEPLER_MAX_ITERATIONS constants
- Separated hyperbolic solver logic from combined function
- Fixed test_newton_raphson_convergence to verify Kepler's equation
  (instead of incorrectly expecting E ≈ M for small e)
- Added TODO comment for future cartesian_to_orbital_elements refactoring

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;
 

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

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);
+
+        // 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);
 
-        double error = fabs(eccentric_anomaly - expected_eccentric_anomaly);
         INFO("Eccentric anomaly: " << eccentric_anomaly << " rad");
-        INFO("Expected: " << expected_eccentric_anomaly << " rad");
-        INFO("Error: " << error);
+        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);
 
-        REQUIRE(error < 1.0e-6);
+        // 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;