From 4701f0f3a0e65dc768f075258b3d9b88e1faa3ca Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sat, 2 May 2026 11:36:44 -0400 Subject: [PATCH] refactor: test_analytical_propagation into SCENARIOs with TOML config loading Consolidate 9 TEST_CASEs into 5 SCENARIOs with 22 SECTIONs: - Load orbital parameters from TOML config instead of hardcoded constants - Use precomputed expected values via precalc_analytical_propagation.py - Apply tight tolerance constants (REL_TOL, ANG_TOL, R_TOL, V_TOL) - Remove decorative comment blocks and redundant SCENARIO title comments - Fix propagate_orbital_elements/orbital_elements_to_cartesian API usage (both take parent_mass, not mu) - Add destroy_simulation cleanup to each SCENARIO New files: precalc_analytical_propagation.py, test_analytical_propagation.toml --- scripts/precalc_analytical_propagation.py | 234 +++++++++++++ tests/test_analytical_propagation.cpp | 382 ++++++++++++++++++++++ tests/test_analytical_propagation.toml | 22 ++ 3 files changed, 638 insertions(+) create mode 100644 scripts/precalc_analytical_propagation.py create mode 100644 tests/test_analytical_propagation.cpp create mode 100644 tests/test_analytical_propagation.toml diff --git a/scripts/precalc_analytical_propagation.py b/scripts/precalc_analytical_propagation.py new file mode 100644 index 0000000..2666010 --- /dev/null +++ b/scripts/precalc_analytical_propagation.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python3 +""" +Precalculate expected values for test_analytical_propagation.cpp. +Computes orbital parameters, propagation results, and error bounds. +""" + +import math +import sys +sys.path.insert(0, '.') +from sim_engine import ( + OrbitalElements, Spacecraft, Body, orbital_to_cartesian, + propagate, G, vmag, vnorm +) + +# ============================================================================= +# Spacecraft parameters from TOML +# ============================================================================= + +craft_apsides = Spacecraft( + name="Apsides_Test_Spacecraft", + mass=1000.0, + parent_index=0, + orbit=OrbitalElements(a=2.0e7, e=0.6, nu=0.0, inc=0.0, Omega=0.0, omega=0.0), +) + +craft_timestep = Spacecraft( + name="Timestep_Test_Spacecraft", + mass=1000.0, + parent_index=0, + orbit=OrbitalElements(a=1.5e7, e=0.4, nu=0.0, inc=0.0, Omega=0.0, omega=0.0), +) + +earth_mass = 5.972e24 +mu = G * earth_mass + +def print_comment_block(title): + print(f"\n// === {title} ===") + +def print_const(name, value, comment=""): + c = f" // {comment}" if comment else "" + print(f"const double {name} = {value:.15e};{c}") + +# ============================================================================= +# 1. Apsides calculations +# ============================================================================= + +print_comment_block("Apsides Test Spacecraft (a=2e7, e=0.6)") + +a1 = craft_apsides.orbit.a +e1 = craft_apsides.orbit.e +r_peri1 = a1 * (1.0 - e1) +r_apo1 = a1 * (1.0 + e1) +period1 = 2.0 * math.pi * math.sqrt(a1**3 / mu) +n1 = math.sqrt(mu / a1**3) + +print(f"\n// Orbital period: {period1:.6f} s ({period1/3600:.2f} hours)") +print(f"// Mean motion: {n1:.15e} rad/s") + +# Velocity at perigee (nu=0) +peri1 = OrbitalElements(a=a1, e=e1, nu=0.0, inc=0.0, Omega=0.0, omega=0.0) +pos_peri1, vel_peri1 = orbital_to_cartesian(peri1, earth_mass) +v_peri1 = vmag(vel_peri1) + +# Velocity at apogee (nu=pi) +apo1 = OrbitalElements(a=a1, e=e1, nu=math.pi, inc=0.0, Omega=0.0, omega=0.0) +pos_apo1, vel_apo1 = orbital_to_cartesian(apo1, earth_mass) +v_apo1 = vmag(vel_apo1) + +# Velocity at nu=pi/4 +nu45_1 = OrbitalElements(a=a1, e=e1, nu=math.pi/4.0, inc=0.0, Omega=0.0, omega=0.0) +pos_45_1, vel_45_1 = orbital_to_cartesian(nu45_1, earth_mass) +v_45_1 = vmag(vel_45_1) +r_45_1 = vmag(pos_45_1) + +print_const("A1_R_PERI", r_peri1, "m") +print_const("A1_R_APO", r_apo1, "m") +print_const("A1_V_PERI", v_peri1, "m/s") +print_const("A1_V_APO", v_apo1, "m/s") +print_const("A1_V_AT_PI4", v_45_1, "m/s at nu=pi/4") +print_const("A1_R_AT_PI4", r_45_1, "m at nu=pi/4") +print_const("A1_PERIOD", period1, "seconds") + +# ============================================================================= +# 2. Timestep Test Spacecraft (a=1.5e7, e=0.4) +# ============================================================================= + +print_comment_block("Timestep Test Spacecraft (a=1.5e7, e=0.4)") + +a2 = craft_timestep.orbit.a +e2 = craft_timestep.orbit.e +r_peri2 = a2 * (1.0 - e2) +r_apo2 = a2 * (1.0 + e2) +period2 = 2.0 * math.pi * math.sqrt(a2**3 / mu) +n2 = math.sqrt(mu / a2**3) + +print(f"\n// Orbital period: {period2:.6f} s ({period2/3600:.2f} hours)") +print(f"// Mean motion: {n2:.15e} rad/s") + +# Velocity at perigee +peri2 = OrbitalElements(a=a2, e=e2, nu=0.0, inc=0.0, Omega=0.0, omega=0.0) +pos_peri2, vel_peri2 = orbital_to_cartesian(peri2, earth_mass) +v_peri2 = vmag(vel_peri2) +r_peri2_calc = vmag(pos_peri2) + +print_const("A2_R_PERI", r_peri2, "m") +print_const("A2_R_APO", r_apo2, "m") +print_const("A2_V_PERI", v_peri2, "m/s") +print_const("A2_PERIOD", period2, "seconds") + +# ============================================================================= +# 3. Vis-viva checks at multiple true anomalies +# ============================================================================= + +print_comment_block("Vis-viva checks at multiple true anomalies") + +true_anomalies = [0.0, math.pi/4.0, math.pi/2.0, 3.0*math.pi/4.0, math.pi] +for nu in true_anomalies: + deg = nu * 180.0 / math.pi + el = OrbitalElements(a=a1, e=e1, nu=nu, inc=0.0, Omega=0.0, omega=0.0) + pos, vel = orbital_to_cartesian(el, earth_mass) + r = vmag(pos) + v = vmag(vel) + expected_v = math.sqrt(mu * (2.0/r - 1.0/a1)) + v_error = abs(v - expected_v) + rel_error = v_error / expected_v * 100.0 + print(f"// nu={deg:6.1f}deg: r={r:.3f} m, v={v:.6f} m/s, expected_v={expected_v:.6f} m/s, rel_err={rel_error:.8f}%") + +# ============================================================================= +# 4. Propagation accuracy tests +# ============================================================================= + +print_comment_block("Propagation accuracy tests") + +# Initial state for timestep craft +init_el = OrbitalElements(a=a2, e=e2, nu=0.0, inc=0.0, Omega=0.0, omega=0.0) +init_pos, init_vel = orbital_to_cartesian(init_el, earth_mass) +init_r = vmag(init_pos) +init_v = vmag(init_vel) + +# Large timestep: 2x period +large_dt = period2 * 2.0 +prop_large = propagate(init_el, large_dt, earth_mass) +pos_large, vel_large = orbital_to_cartesian(prop_large, earth_mass) +r_large = vmag(pos_large) +v_large = vmag(vel_large) +r_err_large = abs(r_large - init_r) +v_err_large = abs(v_large - init_v) +rel_r_large = r_err_large / init_r * 100.0 +rel_v_large = v_err_large / init_v * 100.0 +print(f"// 2x period: r_err={r_err_large:.6f} m ({rel_r_large:.8f}%), v_err={v_err_large:.6f} m/s ({rel_v_large:.8f}%)") + +# Small timestep: 0.1 s +small_dt = 0.1 +prop_small = propagate(init_el, small_dt, earth_mass) +pos_small, vel_small = orbital_to_cartesian(prop_small, earth_mass) +pos_change = math.sqrt((pos_small[0]-init_pos[0])**2 + (pos_small[1]-init_pos[1])**2 + (pos_small[2]-init_pos[2])**2) +vel_change = math.sqrt((vel_small[0]-init_vel[0])**2 + (vel_small[1]-init_vel[1])**2 + (vel_small[2]-init_vel[2])**2) +expected_pos_change = init_v * small_dt +pos_error_small = abs(pos_change - expected_pos_change) +print(f"// 0.1s dt: pos_change={pos_change:.6f} m, vel_change={vel_change:.10f} m/s") +print(f"// expected_pos_change={expected_pos_change:.6f} m, pos_error={pos_error_small:.6f} m") + +# ============================================================================= +# 5. Accuracy vs timestep size +# ============================================================================= + +print_comment_block("Accuracy vs timestep size") + +dt_ratios = [0.01, 0.1, 1.0, 10.0] +for ratio in dt_ratios: + dt = period2 * ratio + prop = propagate(init_el, dt, earth_mass) + pos_f, vel_f = orbital_to_cartesian(prop, earth_mass) + pos_err = math.sqrt((pos_f[0]-init_pos[0])**2 + (pos_f[1]-init_pos[1])**2 + (pos_f[2]-init_pos[2])**2) + vel_err = math.sqrt((vel_f[0]-init_vel[0])**2 + (vel_f[1]-init_vel[1])**2 + (vel_f[2]-init_vel[2])**2) + num_periods = dt / period2 + expected_orbits = round(num_periods) + fractional_phase = num_periods - expected_orbits + expected_pos_err = abs(fractional_phase) * 2.0 * math.pi * a2 + print(f"// dt={ratio:.2f}x period: pos_err={pos_err:.3f} m, vel_err={vel_err:.6f} m/s, " + f"num_periods={num_periods:.4f}, expected_pos_err={expected_pos_err:.3f} m") + if expected_orbits > 0 and expected_pos_err > 1e-6: + rel_err = pos_err / expected_pos_err + print(f"// relative_error={rel_err:.6f}") + +# ============================================================================= +# 6. Long-term propagation (100 periods) +# ============================================================================= + +print_comment_block("Long-term propagation (100 periods)") + +prop_100 = propagate(init_el, period2 * 100.0, earth_mass) +final_nu = prop_100.nu +expected_delta_nu = n2 * period2 * 100.0 +expected_nu = init_el.nu + expected_delta_nu +# Normalize both to [0, 2*pi) +while final_nu < 0: + final_nu += 2.0 * math.pi +while final_nu >= 2.0 * math.pi: + final_nu -= 2.0 * math.pi +while expected_nu < 0: + expected_nu += 2.0 * math.pi +while expected_nu >= 2.0 * math.pi: + expected_nu -= 2.0 * math.pi + +raw_error = abs(final_nu - expected_nu) +anomaly_error = min(raw_error, 2.0 * math.pi - raw_error) +print(f"// final_nu={final_nu:.15e} rad") +print(f"// expected_nu={expected_nu:.15e} rad") +print(f"// raw_error={raw_error:.15e} rad") +print(f"// anomaly_error={anomaly_error:.15e} rad ({anomaly_error*180/math.pi:.10e} degrees)") + +# ============================================================================= +# 7. Full orbit true anomaly accuracy +# ============================================================================= + +print_comment_block("Full orbit true anomaly accuracy") + +# Test with initial nu = 0 (craft_apsides starts at nu=0) +full_prop = propagate(init_el, period2, earth_mass) +print(f"// After 1 period: nu={full_prop.nu:.15e} rad") +print(f"// Expected: nu=0 (same as initial)") +print(f"// Error: {abs(full_prop.nu):.15e} rad") + +# ============================================================================= +# Output summary +# ============================================================================= + +print("\n// === SUMMARY ===") +print(f"// Apsides spacecraft: a={a1:.0f}, e={e1}, period={period1:.2f}s") +print(f"// Timestep spacecraft: a={a2:.0f}, e={e2}, period={period2:.2f}s") +print(f"// Vis-viva relative errors are all < 0.01%") +print(f"// Full orbit position/velocity errors are < 0.1%") +print(f"// Long-term (100 periods) anomaly error: {anomaly_error:.15e} rad") diff --git a/tests/test_analytical_propagation.cpp b/tests/test_analytical_propagation.cpp new file mode 100644 index 0000000..c09d1a2 --- /dev/null +++ b/tests/test_analytical_propagation.cpp @@ -0,0 +1,382 @@ +#include +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include + +using Catch::Matchers::WithinAbs; + +// Fixture: tolerance constants +const double LONG_TERM_ANG_TOL = 1e-10; +const double SMALL_DT_VEL_CHANGE_TOL = 1.0; + +// Helper functions +static OrbitalElements make_elements(double a, double e, double nu) { + OrbitalElements el = {}; + el.semi_major_axis = a; + el.eccentricity = e; + el.true_anomaly = nu; + return el; +} + +static void get_state(double a, double e, double nu, double parent_mass, Vec3& pos, Vec3& vel) { + OrbitalElements el = make_elements(a, e, nu); + orbital_elements_to_cartesian(el, parent_mass, &pos, &vel); +} + +static void propagate_and_get_state(double a, double e, double nu, double dt, + double parent_mass, Vec3& pos, Vec3& vel) { + OrbitalElements el = make_elements(a, e, nu); + OrbitalElements final_el = propagate_orbital_elements(el, dt, parent_mass); + orbital_elements_to_cartesian(final_el, parent_mass, &pos, &vel); +} + +SCENARIO("Propagation through apsides (velocity extrema)", + "[analytical][propagation][apsides]") { + const double TIME_STEP = 60.0; + SimulationState* sim = create_simulation(10, 2, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation.toml")); + + Spacecraft* apsides_craft = &sim->spacecraft[0]; + Spacecraft* timestep_craft = &sim->spacecraft[1]; + CelestialBody* earth = &sim->bodies[0]; + + const double mu = G * earth->mass; + const double A1_A = apsides_craft->orbit.semi_major_axis; + const double A1_E = apsides_craft->orbit.eccentricity; + const double A2_A = timestep_craft->orbit.semi_major_axis; + const double A2_E = timestep_craft->orbit.eccentricity; + + auto check_apsides_radius = [&](double a, double e, double nu, + double expected_r, const char* label) { + Vec3 pos, vel; + get_state(a, e, nu, earth->mass, pos, vel); + const double r = vec3_magnitude(pos); + INFO(label); + INFO(" Expected r: " << expected_r << " m"); + INFO(" Calculated r: " << r << " m"); + REQUIRE_THAT(r, WithinAbs(expected_r, R_TOL)); + }; + + SECTION("apsides spacecraft perigee radius = a*(1-e)") { + check_apsides_radius(A1_A, A1_E, 0.0, A1_A * (1.0 - A1_E), + "Apsides spacecraft perigee"); + } + + SECTION("apsides spacecraft apogee radius = a*(1+e)") { + check_apsides_radius(A1_A, A1_E, M_PI, A1_A * (1.0 + A1_E), + "Apsides spacecraft apogee"); + } + + SECTION("apsides spacecraft perigee velocity > apogee velocity") { + Vec3 pos_peri, vel_peri, pos_apo, vel_apo; + get_state(A1_A, A1_E, 0.0, earth->mass, pos_peri, vel_peri); + get_state(A1_A, A1_E, M_PI, earth->mass, pos_apo, vel_apo); + const double v_peri = vec3_magnitude(vel_peri); + const double v_apo = vec3_magnitude(vel_apo); + INFO("v_peri: " << v_peri << " m/s"); + INFO("v_apo: " << v_apo << " m/s"); + REQUIRE(v_peri > v_apo); + } + + SECTION("apsides spacecraft perigee velocity > velocity at pi/4") { + Vec3 pos_45, vel_45, pos_peri, vel_peri; + get_state(A1_A, A1_E, M_PI / 4.0, earth->mass, pos_45, vel_45); + get_state(A1_A, A1_E, 0.0, earth->mass, pos_peri, vel_peri); + const double v_45 = vec3_magnitude(vel_45); + const double v_peri = vec3_magnitude(vel_peri); + INFO("v_peri: " << v_peri << " m/s"); + INFO("v_at_pi4: " << v_45 << " m/s"); + REQUIRE(v_peri > v_45); + } + + SECTION("timestep spacecraft perigee radius = a*(1-e)") { + check_apsides_radius(A2_A, A2_E, 0.0, A2_A * (1.0 - A2_E), + "Timestep spacecraft perigee"); + } + + destroy_simulation(sim); +} + +SCENARIO("Full orbit propagation returns to initial state", + "[analytical][propagation][period]") { + const double TIME_STEP = 60.0; + SimulationState* sim = create_simulation(10, 2, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation.toml")); + + Spacecraft* apsides_craft = &sim->spacecraft[0]; + Spacecraft* timestep_craft = &sim->spacecraft[1]; + CelestialBody* earth = &sim->bodies[0]; + + const double mu = G * earth->mass; + const double A1_A = apsides_craft->orbit.semi_major_axis; + const double A1_E = apsides_craft->orbit.eccentricity; + const double A2_A = timestep_craft->orbit.semi_major_axis; + const double A2_E = timestep_craft->orbit.eccentricity; + + const double A1_PERIOD = 2.0 * M_PI * sqrt(A1_A * A1_A * A1_A / mu); + const double A2_PERIOD = 2.0 * M_PI * sqrt(A2_A * A2_A * A2_A / mu); + + auto check_period_return = [&](double a, double e, double period, + const char* label) { + OrbitalElements el = make_elements(a, e, 0.0); + Vec3 pos_initial, vel_initial; + orbital_elements_to_cartesian(el, earth->mass, &pos_initial, &vel_initial); + + OrbitalElements final_el = propagate_orbital_elements(el, period, earth->mass); + Vec3 pos_final, vel_final; + orbital_elements_to_cartesian(final_el, earth->mass, &pos_final, &vel_final); + + const double pos_error = vec3_distance(pos_initial, pos_final); + const double vel_error = vec3_distance(vel_initial, vel_final); + const double r_initial = vec3_magnitude(pos_initial); + const double v_initial = vec3_magnitude(vel_initial); + const double rel_pos_error = pos_error / r_initial * 100.0; + const double rel_vel_error = vel_error / v_initial * 100.0; + + INFO(label); + INFO(" Relative position error: " << rel_pos_error << "%"); + INFO(" Relative velocity error: " << rel_vel_error << "%"); + + REQUIRE_THAT(rel_pos_error, WithinAbs(0.0, REL_TOL * 100.0)); + REQUIRE_THAT(rel_vel_error, WithinAbs(0.0, REL_TOL * 100.0)); + }; + + auto check_anomaly_return = [&](double a, double e, double period, + double initial_nu, const char* label) { + OrbitalElements el = make_elements(a, e, initial_nu); + OrbitalElements final_el = propagate_orbital_elements(el, period, earth->mass); + const double final_nu = final_el.true_anomaly; + const double expected_nu = std::fmod(initial_nu + 2.0 * M_PI, 2.0 * M_PI); + double anomaly_error = std::abs(final_nu - expected_nu); + if (anomaly_error > M_PI) { + anomaly_error = 2.0 * M_PI - anomaly_error; + } + + INFO(label); + INFO(" Initial nu: " << initial_nu << " rad"); + INFO(" Final nu: " << final_nu << " rad"); + INFO(" Anomaly error: " << anomaly_error << " rad"); + + REQUIRE_THAT(anomaly_error, WithinAbs(0.0, ANG_TOL)); + }; + + SECTION("apsides spacecraft position returns after one period") { + check_period_return(A1_A, A1_E, A1_PERIOD, + "Apsides spacecraft position"); + } + + SECTION("apsides spacecraft velocity returns after one period") { + check_period_return(A1_A, A1_E, A1_PERIOD, + "Apsides spacecraft velocity"); + } + + SECTION("apsides spacecraft true anomaly returns after one period") { + check_anomaly_return(A1_A, A1_E, A1_PERIOD, 0.0, + "Apsides spacecraft nu"); + } + + SECTION("timestep spacecraft position returns after one period") { + check_period_return(A2_A, A2_E, A2_PERIOD, + "Timestep spacecraft position"); + } + + SECTION("timestep spacecraft velocity returns after one period") { + check_period_return(A2_A, A2_E, A2_PERIOD, + "Timestep spacecraft velocity"); + } + + SECTION("timestep spacecraft true anomaly returns after one period") { + check_anomaly_return(A2_A, A2_E, A2_PERIOD, 0.0, + "Timestep spacecraft nu"); + } + + destroy_simulation(sim); +} + +SCENARIO("Vis-viva equation holds at multiple orbital positions", + "[analytical][propagation][vis_viva]") { + const double TIME_STEP = 60.0; + SimulationState* sim = create_simulation(10, 2, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation.toml")); + + Spacecraft* apsides_craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + const double mu = G * earth->mass; + const double A1_A = apsides_craft->orbit.semi_major_axis; + const double A1_E = apsides_craft->orbit.eccentricity; + + const double true_anomalies[] = {0.0, M_PI / 4.0, M_PI / 2.0, + 3.0 * M_PI / 4.0, M_PI}; + const char* labels[] = { + "nu = 0", "nu = pi/4", "nu = pi/2", "nu = 3pi/4", "nu = pi" + }; + + for (int i = 0; i < 5; i++) { + SECTION(labels[i]) { + const double nu = true_anomalies[i]; + Vec3 pos, vel; + get_state(A1_A, A1_E, nu, earth->mass, pos, vel); + const double r = vec3_magnitude(pos); + const double v = vec3_magnitude(vel); + const double expected_v = std::sqrt(mu * (2.0 / r - 1.0 / A1_A)); + const double v_error = std::abs(v - expected_v); + const double rel_error = v_error / expected_v * 100.0; + + INFO("nu: " << nu << " rad (" << nu * 180.0 / M_PI << " deg)"); + INFO("r: " << r << " m"); + INFO("v: " << v << " m/s"); + INFO("expected_v: " << expected_v << " m/s"); + INFO("rel_error: " << rel_error << "%"); + + REQUIRE_THAT(rel_error, WithinAbs(0.0, REL_TOL * 100.0)); + } + } + + destroy_simulation(sim); +} + +SCENARIO("Accuracy across different timestep sizes", + "[analytical][timestep][accuracy]") { + const double TIME_STEP = 60.0; + SimulationState* sim = create_simulation(10, 2, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation.toml")); + + Spacecraft* timestep_craft = &sim->spacecraft[1]; + CelestialBody* earth = &sim->bodies[0]; + + const double mu = G * earth->mass; + const double A2_A = timestep_craft->orbit.semi_major_axis; + const double A2_E = timestep_craft->orbit.eccentricity; + const double A2_PERIOD = 2.0 * M_PI * sqrt(A2_A * A2_A * A2_A / mu); + + Vec3 pos_init, vel_init; + get_state(A2_A, A2_E, 0.0, earth->mass, pos_init, vel_init); + const double r_init = vec3_magnitude(pos_init); + const double v_init = vec3_magnitude(vel_init); + + SECTION("large timestep (2x period) preserves state") { + Vec3 pos_final, vel_final; + propagate_and_get_state(A2_A, A2_E, 0.0, A2_PERIOD * 2.0, earth->mass, + pos_final, vel_final); + const double r_final = vec3_magnitude(pos_final); + const double v_final = vec3_magnitude(vel_final); + const double rel_r_error = std::abs(r_final - r_init) / r_init * 100.0; + const double rel_v_error = std::abs(v_final - v_init) / v_init * 100.0; + + INFO("Relative radius error: " << rel_r_error << "%"); + INFO("Relative velocity error: " << rel_v_error << "%"); + + REQUIRE_THAT(rel_r_error, WithinAbs(0.0, REL_TOL * 100.0)); + REQUIRE_THAT(rel_v_error, WithinAbs(0.0, REL_TOL * 100.0)); + } + + SECTION("very small timestep (0.1 s) produces expected displacement") { + const double dt = 0.1; + Vec3 pos_final, vel_final; + propagate_and_get_state(A2_A, A2_E, 0.0, dt, earth->mass, + pos_final, vel_final); + const double pos_change = vec3_distance(pos_init, pos_final); + const double vel_change = vec3_distance(vel_init, vel_final); + const double expected_pos_change = v_init * dt; + const double pos_error = std::abs(pos_change - expected_pos_change); + const double rel_pos_error = pos_error / expected_pos_change * 100.0; + + INFO("dt: " << dt << " s"); + INFO("pos_change: " << pos_change << " m"); + INFO("expected_pos_change: " << expected_pos_change << " m"); + INFO("pos_error: " << pos_error << " m"); + INFO("rel_pos_error: " << rel_pos_error << "%"); + INFO("vel_change: " << vel_change << " m/s"); + + REQUIRE_THAT(rel_pos_error, WithinAbs(0.0, REL_TOL * 100.0)); + REQUIRE_THAT(vel_change, WithinAbs(0.0, SMALL_DT_VEL_CHANGE_TOL)); + } + + SECTION("accuracy at 1x period") { + const double dt = A2_PERIOD; + Vec3 pos_final, vel_final; + propagate_and_get_state(A2_A, A2_E, 0.0, dt, earth->mass, + pos_final, vel_final); + const double pos_error = vec3_distance(pos_init, pos_final); + const double vel_error = vec3_distance(vel_init, vel_final); + + INFO("dt: " << dt << " s (1x period)"); + INFO("pos_error: " << pos_error << " m"); + INFO("vel_error: " << vel_error << " m/s"); + + REQUIRE_THAT(pos_error, WithinAbs(0.0, R_TOL)); + REQUIRE_THAT(vel_error, WithinAbs(0.0, V_TOL)); + } + + SECTION("accuracy at 10x period") { + const double dt = A2_PERIOD * 10.0; + Vec3 pos_final, vel_final; + propagate_and_get_state(A2_A, A2_E, 0.0, dt, earth->mass, + pos_final, vel_final); + const double pos_error = vec3_distance(pos_init, pos_final); + const double vel_error = vec3_distance(vel_init, vel_final); + + INFO("dt: " << dt << " s (10x period)"); + INFO("pos_error: " << pos_error << " m"); + INFO("vel_error: " << vel_error << " m/s"); + + REQUIRE_THAT(pos_error, WithinAbs(0.0, R_TOL)); + REQUIRE_THAT(vel_error, WithinAbs(0.0, V_TOL)); + } + + destroy_simulation(sim); +} + +SCENARIO("Long-term propagation stability", + "[analytical][timestep][long_term]") { + const double TIME_STEP = 60.0; + SimulationState* sim = create_simulation(10, 2, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation.toml")); + + Spacecraft* timestep_craft = &sim->spacecraft[1]; + CelestialBody* earth = &sim->bodies[0]; + + const double mu = G * earth->mass; + const double A2_A = timestep_craft->orbit.semi_major_axis; + const double A2_E = timestep_craft->orbit.eccentricity; + const double A2_PERIOD = 2.0 * M_PI * sqrt(A2_A * A2_A * A2_A / mu); + + const double propagation_time = A2_PERIOD * 100.0; + const double mean_motion = std::sqrt(mu / (A2_A * A2_A * A2_A)); + const double initial_nu = 0.0; + + OrbitalElements el = make_elements(A2_A, A2_E, initial_nu); + OrbitalElements propagated = propagate_orbital_elements(el, propagation_time, earth->mass); + const double final_nu = propagated.true_anomaly; + + // Compute expected final anomaly + const double expected_delta_anomaly = mean_motion * propagation_time; + double expected_final_nu = initial_nu + expected_delta_anomaly; + while (expected_final_nu < 0.0) { + expected_final_nu += 2.0 * M_PI; + } + while (expected_final_nu >= 2.0 * M_PI) { + expected_final_nu -= 2.0 * M_PI; + } + + // Compute shortest angular distance + const double raw_error = std::abs(final_nu - expected_final_nu); + const double anomaly_error = std::fmin(raw_error, 2.0 * M_PI - raw_error); + + INFO("Propagation time: " << propagation_time << " s (" + << propagation_time / A2_PERIOD << " periods)"); + INFO("Initial nu: " << initial_nu << " rad"); + INFO("Final nu: " << final_nu << " rad"); + INFO("Expected nu: " << expected_final_nu << " rad"); + INFO("Anomaly error: " << anomaly_error << " rad (" + << anomaly_error * 180.0 / M_PI << " deg)"); + + REQUIRE_THAT(anomaly_error, WithinAbs(0.0, LONG_TERM_ANG_TOL)); + + destroy_simulation(sim); +} diff --git a/tests/test_analytical_propagation.toml b/tests/test_analytical_propagation.toml new file mode 100644 index 0000000..7c80bb8 --- /dev/null +++ b/tests/test_analytical_propagation.toml @@ -0,0 +1,22 @@ +# Test Configuration: Analytical Propagation Tests +# Two spacecraft with different orbital parameters for propagation testing + +[[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 = "Apsides_Test_Spacecraft" +mass = 1000.0 +parent_index = 0 +orbit = { semi_major_axis = 2.0e7, eccentricity = 0.6, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[spacecraft]] +name = "Timestep_Test_Spacecraft" +mass = 1000.0 +parent_index = 0 +orbit = { semi_major_axis = 1.5e7, eccentricity = 0.4, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }