feat: Implement Barker's equation for parabolic orbit propagation

- Add solve_barker_equation() function using cubic formula: D + D³/3 = M - Integrate Barker's equation into propagate_orbital_elements() for parabolic orbits - Add comprehensive test suite (11 tests, 239 assertions) following TDD - Use cbrt() for cube root (handles negative numbers properly) - Parabolic propagation now uses exact analytical solution instead of iterative solver - All 93 tests passing (239,872 assertions)
5 months ago · 5eb3a84473
3 changed files with 205 additions and 18 deletions
--- a/src/orbital_mechanics.cpp
+++ b/src/orbital_mechanics.cpp
@ -130,6 +130,20 @@ double mean_anomaly_to_true_anomaly(double mean_anomaly, double eccentricity) {
    }
 }

+double solve_barker_equation(double mean_anomaly) {
+    if (fabs(mean_anomaly) < 1e-15) {
+        return 0.0;
+    }
+
+    double c = 1.5 * mean_anomaly;
+    double discriminant = c * c + 1.0;
+    double sqrt_discriminant = sqrt(discriminant);
+    double D = cbrt(c + sqrt_discriminant) + cbrt(c - sqrt_discriminant);
+    double nu = 2.0 * atan(D);
+
+    return nu;
+}
+
 // TODO: refactor for readability
 OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, double parent_mass) {
    double mu = G * parent_mass;
@ -251,30 +265,45 @@ OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, doub
    double e = elements.eccentricity;
    double nu = elements.true_anomaly;
    double mu = G * parent_mass;
-    double n = sqrt(mu / pow(fabs(a), 3.0));

-    double E = 2.0 * atan(sqrt((1.0 - e) / (1.0 + e)) * tan(nu / 2.0));
+    if (fabs(e - 1.0) < PARABOLIC_TOLERANCE) {
+        double p = elements.semi_latus_rectum;
+        double D = tan(nu / 2.0);
+        double M = D + (D * D * D) / 3.0;
+        double n = sqrt(mu / pow(p, 3.0));
+        M = M + n * dt;
+        double nu_new = solve_barker_equation(M);
+
+        OrbitalElements result = elements;
+        result.true_anomaly = nu_new;
+        return result;
+    } else {
+        double p = a * (1.0 - e * e);
+        double n = sqrt(mu / pow(fabs(a), 3.0));

-    double M = E - e * sin(E);
+        double E = 2.0 * atan(sqrt((1.0 - e) / (1.0 + e)) * tan(nu / 2.0));

-    M = M + n * dt;
+        double M = E - e * sin(E);

-    double E_new = get_initial_trial_value(M, e);
+        M = M + n * dt;

-    const double CONVERGENCE_TOLERANCE = 1.0e-10;
-    const int MAX_ITERATIONS = 50;
+        double E_new = get_initial_trial_value(M, e);

-    int iterations = 0;
-    double E_prev = E_new + 2.0 * CONVERGENCE_TOLERANCE;
-    while (fabs(E_new - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) {
-        E_prev = E_new;
-        double sin_E = sin(E_new);
-        E_new = E_new - (E_new - e * sin_E - M) / (1.0 - e * cos(E_new));
-        iterations++;
-    }
+        const double CONVERGENCE_TOLERANCE = 1.0e-10;
+        const int MAX_ITERATIONS = 50;

-    OrbitalElements result = elements;
-    result.true_anomaly = 2.0 * atan(sqrt((1.0 + e) / (1.0 - e)) * tan(E_new / 2.0));
+        int iterations = 0;
+        double E_prev = E_new + 2.0 * CONVERGENCE_TOLERANCE;
+        while (fabs(E_new - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) {
+            E_prev = E_new;
+            double sin_E = sin(E_new);
+            E_new = E_new - (E_new - e * sin_E - M) / (1.0 - e * cos(E_new));
+            iterations++;
+        }
+
+        OrbitalElements result = elements;
+        result.true_anomaly = 2.0 * atan(sqrt((1.0 + e) / (1.0 - e)) * tan(E_new / 2.0));

-    return result;
+        return result;
+    }
 }
--- a/src/orbital_mechanics.h
+++ b/src/orbital_mechanics.h
@ -39,6 +39,9 @@ double hyperbolic_to_true_anomaly(double hyperbolic_anomaly, double eccentricity
 // Automatically dispatches to elliptical or hyperbolic based on eccentricity
 double mean_anomaly_to_true_anomaly(double mean_anomaly, double eccentricity);

+// Barker's equation: D + D³/3 = M, where D = tan(ν/2)
+double solve_barker_equation(double mean_anomaly);
+
 OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, double dt, double parent_mass);

 #endif
--- a/tests/test_barkers_equation.cpp
+++ b/tests/test_barkers_equation.cpp
@ -0,0 +1,155 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include "../src/orbital_mechanics.h"
+#include <cmath>
+
+TEST_CASE("Barker's equation - zero mean anomaly", "[barker][analytical]") {
+    double M = 0.0;
+    double nu = solve_barker_equation(M);
+    double nu_expected = 0.0;
+    REQUIRE_THAT(nu, Catch::Matchers::WithinAbs(nu_expected, 1e-15));
+}
+
+TEST_CASE("Barker's equation - small positive mean anomaly", "[barker][analytical]") {
+    double M = 0.1;
+    double nu = solve_barker_equation(M);
+    REQUIRE(nu > 0.0);
+    REQUIRE(nu < M_PI);
+    double D = tan(nu / 2.0);
+    double M_recovered = D + (D * D * D) / 3.0;
+    REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14));
+}
+
+TEST_CASE("Barker's equation - moderate positive mean anomaly", "[barker][analytical]") {
+    double M = 1.0;
+    double nu = solve_barker_equation(M);
+    REQUIRE(nu > 0.0);
+    REQUIRE(nu < M_PI);
+    double D = tan(nu / 2.0);
+    double M_recovered = D + (D * D * D) / 3.0;
+    REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14));
+}
+
+TEST_CASE("Barker's equation - large positive mean anomaly", "[barker][analytical]") {
+    double M = 5.0;
+    double nu = solve_barker_equation(M);
+    REQUIRE(nu > 0.0);
+    REQUIRE(nu < M_PI);
+    double D = tan(nu / 2.0);
+    double M_recovered = D + (D * D * D) / 3.0;
+    REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14));
+}
+
+TEST_CASE("Barker's equation - very large mean anomaly", "[barker][analytical]") {
+    double M = 20.0;
+    double nu = solve_barker_equation(M);
+    REQUIRE(nu > 0.0);
+    REQUIRE(nu < M_PI);
+    double D = tan(nu / 2.0);
+    double M_recovered = D + (D * D * D) / 3.0;
+    REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-13));
+}
+
+TEST_CASE("Barker's equation - small negative mean anomaly", "[barker][analytical]") {
+    double M = -0.1;
+    double nu = solve_barker_equation(M);
+    REQUIRE(nu < 0.0);
+    REQUIRE(nu > -M_PI);
+    double D = tan(nu / 2.0);
+    double M_recovered = D + (D * D * D) / 3.0;
+    REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14));
+}
+
+TEST_CASE("Barker's equation - moderate negative mean anomaly", "[barker][analytical]") {
+    double M = -1.0;
+    double nu = solve_barker_equation(M);
+    REQUIRE(nu < 0.0);
+    REQUIRE(nu > -M_PI);
+    double D = tan(nu / 2.0);
+    double M_recovered = D + (D * D * D) / 3.0;
+    REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14));
+}
+
+TEST_CASE("Barker's equation - large negative mean anomaly", "[barker][analytical]") {
+    double M = -5.0;
+    double nu = solve_barker_equation(M);
+    REQUIRE(nu < 0.0);
+    REQUIRE(nu > -M_PI);
+    double D = tan(nu / 2.0);
+    double M_recovered = D + (D * D * D) / 3.0;
+    REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M, 1e-14));
+}
+
+TEST_CASE("Barker's equation - round-trip conversion", "[barker][analytical]") {
+    std::vector<double> test_values = {-10.0, -5.0, -1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0, 5.0, 10.0};
+
+    for (double M_original : test_values) {
+        double nu = solve_barker_equation(M_original);
+        double D = tan(nu / 2.0);
+        double M_recovered = D + (D * D * D) / 3.0;
+        REQUIRE_THAT(M_recovered, Catch::Matchers::WithinAbs(M_original, 1e-13));
+    }
+}
+
+TEST_CASE("Barker's equation - true anomaly range", "[barker][analytical]") {
+    for (double M = -50.0; M <= 50.0; M += 1.0) {
+        double nu = solve_barker_equation(M);
+        REQUIRE(nu > -M_PI * 0.99);
+        REQUIRE(nu < M_PI * 0.99);
+    }
+}
+
+TEST_CASE("Parabolic orbit propagation using Barker's equation", "[barker][propagation]") {
+    const double PARENT_MASS = 1.989e30;
+    const double TIME_STEP = 3600.0;
+    const int NUM_STEPS = 24;
+
+    OrbitalElements initial;
+    initial.semi_latus_rectum = 2.992e11;
+    initial.eccentricity = 1.0;
+    initial.true_anomaly = 0.0;
+    initial.inclination = 0.0;
+    initial.longitude_of_ascending_node = 0.0;
+    initial.argument_of_periapsis = 0.0;
+
+    Vec3 pos, vel;
+    orbital_elements_to_cartesian(initial, PARENT_MASS, &pos, &vel);
+
+    double initial_distance = vec3_magnitude(pos);
+    double initial_velocity = vec3_magnitude(vel);
+    double escape_velocity = sqrt(2.0 * G * PARENT_MASS / initial_distance);
+
+    INFO("Initial distance: " << initial_distance / 1.496e11 << " AU");
+    INFO("Initial velocity: " << initial_velocity / 1000.0 << " km/s");
+    INFO("Escape velocity: " << escape_velocity / 1000.0 << " km/s");
+
+    REQUIRE_THAT(initial_velocity, Catch::Matchers::WithinAbs(escape_velocity, 1.0));
+
+    OrbitalElements current = initial;
+    double total_time = 0.0;
+
+    for (int step = 0; step < NUM_STEPS; step++) {
+        OrbitalElements next = propagate_orbital_elements(current, TIME_STEP, PARENT_MASS);
+        current = next;
+        total_time += TIME_STEP;
+    }
+
+    Vec3 pos_final, vel_final;
+    orbital_elements_to_cartesian(current, PARENT_MASS, &pos_final, &vel_final);
+
+    double final_distance = vec3_magnitude(pos_final);
+    double final_velocity = vec3_magnitude(vel_final);
+
+    INFO("Final true anomaly: " << current.true_anomaly << " rad");
+    INFO("Final distance: " << final_distance / 1.496e11 << " AU");
+    INFO("Final velocity: " << final_velocity / 1000.0 << " km/s");
+
+    REQUIRE(final_distance > initial_distance);
+
+    REQUIRE(final_velocity < initial_velocity);
+
+    double final_escape_velocity = sqrt(2.0 * G * PARENT_MASS / final_distance);
+    INFO("Final escape velocity: " << final_escape_velocity / 1000.0 << " km/s");
+
+    REQUIRE_THAT(final_velocity, Catch::Matchers::WithinAbs(final_escape_velocity, 1.0));
+}