From 5eb3a844730aa6429952180a9de86baccf07ce28 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sun, 1 Feb 2026 15:10:03 -0500 Subject: [PATCH] feat: Implement Barker's equation for parabolic orbit propagation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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) --- src/orbital_mechanics.cpp | 65 ++++++++++---- src/orbital_mechanics.h | 3 + tests/test_barkers_equation.cpp | 155 ++++++++++++++++++++++++++++++++ 3 files changed, 205 insertions(+), 18 deletions(-) create mode 100644 tests/test_barkers_equation.cpp diff --git a/src/orbital_mechanics.cpp b/src/orbital_mechanics.cpp index 4ebeeb3..94be15b 100644 --- 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; + } } diff --git a/src/orbital_mechanics.h b/src/orbital_mechanics.h index 4720d52..b17c1ab 100644 --- 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 diff --git a/tests/test_barkers_equation.cpp b/tests/test_barkers_equation.cpp new file mode 100644 index 0000000..fcee8cc --- /dev/null +++ b/tests/test_barkers_equation.cpp @@ -0,0 +1,155 @@ +#include +#include +#include "../src/orbital_mechanics.h" +#include + +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 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)); +}