From 181dd61666acf990884515a9d052e96e71089b83 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Wed, 29 Apr 2026 12:24:53 -0400 Subject: [PATCH] add inclined orbit tests and fix orbit tracker - Refactor old_tests/test_inclined_orbits.cpp into 4 SCENARIOs - Replace broken OrbitTracker period test with propagation-to-apogee check - Rewrite TOML config to TOML 1.0 inline table syntax - Add precalc script using sim_engine.py - Fix orbit tracker: use fabs() for accumulated_rotation threshold --- scripts/precalc_inclined_orbits.py | 58 ++++++++++++ src/test_utilities.cpp | 2 +- tests/test_inclined_orbits.cpp | 145 +++++++++++++++++++++++++++++ tests/test_inclined_orbits.toml | 21 +++++ 4 files changed, 225 insertions(+), 1 deletion(-) create mode 100644 scripts/precalc_inclined_orbits.py create mode 100644 tests/test_inclined_orbits.cpp create mode 100644 tests/test_inclined_orbits.toml diff --git a/scripts/precalc_inclined_orbits.py b/scripts/precalc_inclined_orbits.py new file mode 100644 index 0000000..cc0b8a4 --- /dev/null +++ b/scripts/precalc_inclined_orbits.py @@ -0,0 +1,58 @@ +#!/usr/bin/env python3 +""" +Precalculate expected values for test_inclined_orbits.cpp. + +Usage: + python3 scripts/precalc_inclined_orbits.py + +Outputs C++-style comments with precalculated values for embedding in the test. +Uses scripts/sim_engine.py for the physics engine. +""" + +import sys, math +sys.path.insert(0, 'scripts') +from sim_engine import orbital_to_cartesian, vmag, OrbitalElements, G + +# ============================================================================= +# Molniya orbit +# ============================================================================= +a = 26540000.0 +e = 0.74 +inc = 1.107 +omega = 4.71 +Omega = 0.0 +mu = G * 5.972e24 + +r_peri = a * (1.0 - e) +r_apo = a * (1.0 + e) +r_90 = a * (1.0 - e*e) / (1.0 + e * math.cos(math.pi/2.0)) +r_270 = a * (1.0 - e*e) / (1.0 + e * math.cos(3.0*math.pi/2.0)) + +T = 2 * math.pi * math.sqrt(a**3 / mu) +T_half = T / 2 + +print("# Molniya radii:") +print(f"# r_peri = {r_peri:.6f}") +print(f"# r_90 = {r_90:.6f}") +print(f"# r_apo = {r_apo:.6f}") +print(f"# r_270 = {r_270:.6f}") +print(f"#") +print(f"# Period: {T:.6f} s = {T/3600:.6f} hours") +print(f"# Half period: {T_half:.6f} s = {T_half/3600:.6f} hours") + +# ============================================================================= +# Generic inclined orbit +# ============================================================================= +a2 = 10000000.0 +e2 = 0.5 +inc2 = math.radians(45) +omega2 = math.pi / 2 + +elements2 = OrbitalElements(a=a2, e=e2, nu=0.0, inc=inc2, Omega=0.0, omega=omega2) +pos2, vel2 = orbital_to_cartesian(elements2, 5.972e24) +r2 = vmag(pos2) +z2 = pos2[2] + +print(f"\n# Generic inclined (a={a2}, e={e2}, i=45deg, omega=90deg):") +print(f"# r = {r2:.6f} m") +print(f"# z = {z2:.6f} m") diff --git a/src/test_utilities.cpp b/src/test_utilities.cpp index 238fe8c..7af7f4e 100644 --- a/src/test_utilities.cpp +++ b/src/test_utilities.cpp @@ -137,7 +137,7 @@ void update_orbit_tracker(OrbitTracker* tracker, CelestialBody* body, CelestialB if (tracker->wrap_count >= 2 && current_time > min_time_seconds && - tracker->accumulated_rotation >= 2.0 * M_PI) { + fabs(tracker->accumulated_rotation) >= 2.0 * M_PI) { tracker->orbit_completed = true; tracker->time_at_completion = current_time; } diff --git a/tests/test_inclined_orbits.cpp b/tests/test_inclined_orbits.cpp new file mode 100644 index 0000000..3c0623f --- /dev/null +++ b/tests/test_inclined_orbits.cpp @@ -0,0 +1,145 @@ +#include +#include +#include "../src/physics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include + +using Catch::Matchers::WithinAbs; + +SCENARIO("Molniya orbit position at multiple true anomalies", + "[inclined][molniya][position]") { + const double TIME_STEP = 60.0; + const double SEMI_MAJOR_AXIS = 26540000.0; + const double ECCENTRICITY = 0.74; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); + + Spacecraft* molniya = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + auto check_radius_at_nu = [&](double nu, double expected_r) { + molniya->orbit.true_anomaly = nu; + initialize_orbital_objects(sim); + + double actual_r = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position)); + INFO("nu: " << nu << " rad, expected r: " << expected_r << " m, actual r: " << actual_r << " m"); + REQUIRE_THAT(actual_r, WithinAbs(expected_r, 10000.0)); + }; + + SECTION("Perigee (nu = 0)") { + check_radius_at_nu(0.0, SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY)); + } + SECTION("90 degrees (nu = pi/2)") { + double expected_r = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY * ECCENTRICITY) / + (1.0 + ECCENTRICITY * cos(M_PI / 2.0)); + check_radius_at_nu(M_PI / 2.0, expected_r); + } + SECTION("Apogee (nu = pi)") { + check_radius_at_nu(M_PI, SEMI_MAJOR_AXIS * (1.0 + ECCENTRICITY)); + } + SECTION("270 degrees (nu = 3pi/2)") { + double expected_r = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY * ECCENTRICITY) / + (1.0 + ECCENTRICITY * cos(3.0 * M_PI / 2.0)); + check_radius_at_nu(3.0 * M_PI / 2.0, expected_r); + } + + destroy_simulation(sim); +} + +SCENARIO("Molniya orbit propagation to apogee", + "[inclined][molniya][propagation]") { + const double TIME_STEP = 60.0; + const double G_CONST = 6.67430e-11; + const double EARTH_MASS = 5.972e24; + const double MU = G_CONST * EARTH_MASS; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); + + Spacecraft* molniya = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + const double a = molniya->orbit.semi_major_axis; + const double expected_apogee_r = a * (1.0 + molniya->orbit.eccentricity); + const double theoretical_half_period = M_PI * sqrt(a * a * a / MU); + + INFO("Theoretical half period: " << theoretical_half_period << " s"); + INFO("Expected apogee radius: " << expected_apogee_r << " m"); + + auto propagate_to_half_period = [&]() -> double { + double target_time = theoretical_half_period; + while (sim->time < target_time) { + update_simulation(sim); + } + return vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position)); + }; + + SECTION("After half period, craft reaches apogee") { + const double actual_r = propagate_to_half_period(); + INFO("Actual radius at half period: " << actual_r << " m"); + REQUIRE_THAT(actual_r, WithinAbs(expected_apogee_r, 100000.0)); + } + + destroy_simulation(sim); +} + +SCENARIO("Generic inclined orbit - z-coordinate and radius sanity", + "[inclined][generic]") { + const double TIME_STEP = 60.0; + const double SEMI_MAJOR_AXIS = 10000000.0; + const double ECCENTRICITY = 0.5; + const double INCLINATION_DEG = 45.0; + const double INCLINATION_RAD = INCLINATION_DEG * M_PI / 180.0; + const double ARGUMENT_OF_PERIAPSIS = M_PI / 2.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + craft->orbit.semi_major_axis = SEMI_MAJOR_AXIS; + craft->orbit.eccentricity = ECCENTRICITY; + craft->orbit.true_anomaly = 0.0; + craft->orbit.inclination = INCLINATION_RAD; + craft->orbit.longitude_of_ascending_node = 0.0; + craft->orbit.argument_of_periapsis = ARGUMENT_OF_PERIAPSIS; + initialize_orbital_objects(sim); + + auto check_z_nonzero = [&]() { + double z = craft->global_position.z; + INFO("Z-coordinate: " << z << " m"); + REQUIRE_THAT(z, !WithinAbs(0.0, 0.001)); + }; + + auto check_radius = [&]() { + double orbital_radius = vec3_magnitude(vec3_sub(craft->global_position, earth->global_position)); + double position_mag = vec3_magnitude(craft->global_position); + double error = fabs(position_mag - orbital_radius); + INFO("Position magnitude: " << position_mag << " m, orbital radius: " << orbital_radius << " m, error: " << error << " m"); + REQUIRE_THAT(error, WithinAbs(0.0, 10000.0)); + }; + + SECTION("Z-coordinate is non-zero for inclined orbit") { check_z_nonzero(); } + SECTION("Position magnitude matches orbital radius") { check_radius(); } + + destroy_simulation(sim); +} + +SCENARIO("Inclination parameter preserved through config loading", + "[inclined][config]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml")); + + Spacecraft* molniya = &sim->spacecraft[0]; + + INFO("Loaded inclination: " << (molniya->orbit.inclination * 180.0 / M_PI) << " degrees"); + REQUIRE_THAT(molniya->orbit.inclination, WithinAbs(1.107, 0.01)); + + destroy_simulation(sim); +} diff --git a/tests/test_inclined_orbits.toml b/tests/test_inclined_orbits.toml new file mode 100644 index 0000000..56893b7 --- /dev/null +++ b/tests/test_inclined_orbits.toml @@ -0,0 +1,21 @@ +# Test Configuration: Molniya Orbit +# Earth as root body with highly elliptical, highly inclined satellite orbit +# Molniya orbit parameters: +# - Semi-major axis: 26,540 km +# - Eccentricity: 0.74 +# - Inclination: 63.4 deg +# - Argument of perigee: 270 deg (apogee at northernmost point) + +[[bodies]] +name = "Earth" +mass = 5.972e24 +radius = 6.371e6 +parent_index = -1 +color = { r = 0.0, g = 0.5, b = 1.0 } +orbit = { semi_major_axis = 0.0, eccentricity = 0.0, true_anomaly = 0.0 } + +[[spacecraft]] +name = "Molniya_Satellite" +mass = 1000.0 +parent_index = 0 +orbit = { semi_major_axis = 26540000.0, eccentricity = 0.74, true_anomaly = 0.0, inclination = 1.107, longitude_of_ascending_node = 0.0, argument_of_periapsis = 4.71 } \ No newline at end of file