From 432f8db7bb6be9225ac7dff1581d327e267ef13c Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sat, 2 May 2026 19:35:58 -0400 Subject: [PATCH] refactor: consolidate test_analytical_propagation into single SCENARIO with precalculated values MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Merge 5 SCENARIOs into 1 with shared fixture - Reorganize 23 SECTIONs by test type (geometry → vis-viva → period → timestep → long-term) - Add 2 missing timestep tests (apogee radius, velocity comparison) - Replace all qualitative checks (a > b) with WithinAbs against precalculated values - Use named tolerance constants everywhere — zero hardcoded numbers in WithinAbs - Add VEL_CHANGE_SMALL_DT fixture constant - Refactor precalc script to load from TOML via Simulator class - Output matches new test ordering --- scripts/precalc_analytical_propagation.py | 437 +++++++++++----------- tests/test_analytical_propagation.cpp | 402 ++++++++++---------- 2 files changed, 400 insertions(+), 439 deletions(-) diff --git a/scripts/precalc_analytical_propagation.py b/scripts/precalc_analytical_propagation.py index 2666010..e8ba745 100644 --- a/scripts/precalc_analytical_propagation.py +++ b/scripts/precalc_analytical_propagation.py @@ -1,234 +1,217 @@ #!/usr/bin/env python3 """ Precalculate expected values for test_analytical_propagation.cpp. -Computes orbital parameters, propagation results, and error bounds. +Loads config from tests/test_analytical_propagation.toml, then 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") +sys.path.insert(0, "scripts") +from sim_engine import Simulator, propagate, orbital_to_cartesian, vmag, G + +def main(): + sim = Simulator("tests/test_analytical_propagation.toml", dt=60.0) + + earth = sim.get_body("Earth") + craft_apsides = sim.get_craft("Apsides_Test_Spacecraft") + craft_timestep = sim.get_craft("Timestep_Test_Spacecraft") + + earth_mass = earth.mass + mu = G * earth_mass + + a1 = craft_apsides.orbit.a + e1 = craft_apsides.orbit.e + period1 = 2.0 * math.pi * math.sqrt(a1**3 / mu) + n1 = math.sqrt(mu / a1**3) + + a2 = craft_timestep.orbit.a + e2 = craft_timestep.orbit.e + period2 = 2.0 * math.pi * math.sqrt(a2**3 / mu) + n2 = math.sqrt(mu / a2**3) + + 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 geometry — both spacecraft + # ============================================================================= + + print_comment_block("Apsides geometry (both spacecraft)") + + # Apsides spacecraft + r_peri1 = a1 * (1.0 - e1) + r_apo1 = a1 * (1.0 + e1) + peri1 = craft_apsides.orbit + _, vel_peri1 = orbital_to_cartesian(peri1, earth_mass) + v_peri1 = vmag(vel_peri1) + apo1_el = type(peri1)(a=a1, e=e1, nu=math.pi, inc=0.0, Omega=0.0, omega=0.0) + _, vel_apo1 = orbital_to_cartesian(apo1_el, earth_mass) + v_apo1 = vmag(vel_apo1) + nu45_1_el = type(peri1)(a=a1, e=e1, nu=math.pi/4.0, inc=0.0, Omega=0.0, omega=0.0) + _, vel_45_1 = orbital_to_cartesian(nu45_1_el, earth_mass) + v_45_1 = vmag(vel_45_1) + + print(f"// Apsides spacecraft: a={a1:.0f}, e={e1}, period={period1:.2f}s") + print(f"// Mean motion: {n1:.15e} rad/s") + print(f"// r_peri={r_peri1:.3f} m, r_apo={r_apo1:.3f} m") + print(f"// v_peri={v_peri1:.6f} m/s, v_apo={v_apo1:.6f} m/s") + print(f"// v_at_pi4={v_45_1:.6f} m/s") + + # Timestep spacecraft + r_peri2 = a2 * (1.0 - e2) + r_apo2 = a2 * (1.0 + e2) + peri2 = craft_timestep.orbit + _, vel_peri2 = orbital_to_cartesian(peri2, earth_mass) + v_peri2 = vmag(vel_peri2) + apo2_el = type(peri2)(a=a2, e=e2, nu=math.pi, inc=0.0, Omega=0.0, omega=0.0) + _, vel_apo2 = orbital_to_cartesian(apo2_el, earth_mass) + v_apo2 = vmag(vel_apo2) + + print(f"// Timestep spacecraft: a={a2:.0f}, e={e2}, period={period2:.2f}s") + print(f"// Mean motion: {n2:.15e} rad/s") + print(f"// r_peri={r_peri2:.3f} m, r_apo={r_apo2:.3f} m") + print(f"// v_peri={v_peri2:.6f} m/s, v_apo={v_apo2:.6f} m/s") + + 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_PERIOD", period1, "seconds") + 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_V_APO", v_apo2, "m/s") + print_const("A2_PERIOD", period2, "seconds") + + # ============================================================================= + # 2. 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 = type(craft_apsides.orbit)(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}%") + + # ============================================================================= + # 3. Period return — full orbit closure for both spacecraft + # ============================================================================= + + print_comment_block("Period return — full orbit closure") + + # Apsides spacecraft: propagate 1 period from nu=0 + el1 = type(craft_apsides.orbit)(a=a1, e=e1, nu=0.0, inc=0.0, Omega=0.0, omega=0.0) + _, vel1_init = orbital_to_cartesian(el1, earth_mass) + prop1 = propagate(el1, period1, earth_mass) + _, vel1_final = orbital_to_cartesian(prop1, earth_mass) + vel_change1 = math.sqrt((vel1_final[0]-vel1_init[0])**2 + (vel1_final[1]-vel1_init[1])**2 + (vel1_final[2]-vel1_init[2])**2) + print(f"// Apsides after 1 period: vel_change={vel_change1:.15e} m/s, final_nu={prop1.nu:.15e} rad") + + # Timestep spacecraft: propagate 1 period from nu=0 + el2 = type(craft_timestep.orbit)(a=a2, e=e2, nu=0.0, inc=0.0, Omega=0.0, omega=0.0) + _, vel2_init = orbital_to_cartesian(el2, earth_mass) + prop2 = propagate(el2, period2, earth_mass) + _, vel2_final = orbital_to_cartesian(prop2, earth_mass) + vel_change2 = math.sqrt((vel2_final[0]-vel2_init[0])**2 + (vel2_final[1]-vel2_init[1])**2 + (vel2_final[2]-vel2_init[2])**2) + print(f"// Timestep after 1 period: vel_change={vel_change2:.15e} m/s, final_nu={prop2.nu:.15e} rad") + + # ============================================================================= + # 4. Timestep accuracy + # ============================================================================= + + print_comment_block("Timestep accuracy") + + # Initial state for timestep craft + init_el = type(craft_timestep.orbit)(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") + + # Accuracy at various multiples of period + dt_ratios = [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) + print(f"// {ratio:.0f}x period: pos_err={pos_err:.6f} m, vel_err={vel_err:.10f} m/s") + + # ============================================================================= + # 5. Long-term stability (100 periods) + # ============================================================================= + + print_comment_block("Long-term stability (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"// Propagation time: {period2*100.0:.2f} s ({100.0} periods)") + 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)") + + # ============================================================================= + # 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") + +if __name__ == "__main__": + main() diff --git a/tests/test_analytical_propagation.cpp b/tests/test_analytical_propagation.cpp index c09d1a2..229e9b1 100644 --- a/tests/test_analytical_propagation.cpp +++ b/tests/test_analytical_propagation.cpp @@ -9,33 +9,9 @@ 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]") { +SCENARIO("Analytical propagation: apsides, periods, vis-viva, timestep accuracy, long-term stability", + "[analytical][propagation][apsides][period][vis_viva][timestep][accuracy][long_term]") { + // === Fixture === 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")); @@ -50,10 +26,46 @@ SCENARIO("Propagation through apsides (velocity extrema)", 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); + + // Precalculated velocity magnitudes (from scripts/precalc_analytical_propagation.py) + const double A1_V_PERI = 8928.484709064580e+00; // m/s at perigee + const double A1_V_APO = 2232.121177266145e+00; // m/s at apogee + const double A1_V_AT_PI4 = 8292.953779787020e+00; // m/s at nu=pi/4 + const double A2_V_PERI = 7874.183374942587e+00; // m/s at perigee + const double A2_V_APO = 3374.650017832537e+00; // m/s at apogee + + // Precalculated displacement for 0.1s timestep (from precalc script) + const double VEL_CHANGE_SMALL_DT = 0.4920854266; + + const double LONG_TERM_ANG_TOL = 1e-10; + + // === Helper lambdas === + auto make_elements = [&](double a, double e, double nu) { + OrbitalElements el = {}; + el.semi_major_axis = a; + el.eccentricity = e; + el.true_anomaly = nu; + return el; + }; + + auto get_state = [&](double a, double e, double nu, Vec3& pos, Vec3& vel) { + OrbitalElements el = make_elements(a, e, nu); + orbital_elements_to_cartesian(el, earth->mass, &pos, &vel); + }; + + auto propagate_and_get_state = [&](double a, double e, double nu, double dt, + Vec3& pos, Vec3& vel) { + OrbitalElements el = make_elements(a, e, nu); + OrbitalElements final_el = propagate_orbital_elements(el, dt, earth->mass); + orbital_elements_to_cartesian(final_el, earth->mass, &pos, &vel); + }; + 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); + get_state(a, e, nu, pos, vel); const double r = vec3_magnitude(pos); INFO(label); INFO(" Expected r: " << expected_r << " m"); @@ -61,65 +73,6 @@ SCENARIO("Propagation through apsides (velocity extrema)", 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); @@ -164,107 +117,146 @@ SCENARIO("Full orbit propagation returns to initial state", 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"); + auto check_vis_viva = [&](double a, double e, double nu, const char* label) { + Vec3 pos, vel; + get_state(a, e, nu, 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 / a)); + const double rel_error = std::abs(v - expected_v) / expected_v * 100.0; + + INFO(label); + 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)); + }; + + // === SECTIONs === + + // --- 1. Apsides geometry (both spacecraft) --- + SECTION("apsides 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 velocity returns after one period") { - check_period_return(A1_A, A1_E, A1_PERIOD, - "Apsides spacecraft velocity"); + SECTION("apsides 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 true anomaly returns after one period") { - check_anomaly_return(A1_A, A1_E, A1_PERIOD, 0.0, - "Apsides spacecraft nu"); + SECTION("apsides perigee velocity > apogee velocity") { + Vec3 pos_peri, vel_peri, pos_apo, vel_apo; + get_state(A1_A, A1_E, 0.0, pos_peri, vel_peri); + get_state(A1_A, A1_E, M_PI, 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_THAT(v_peri, WithinAbs(A1_V_PERI, V_TOL)); + REQUIRE_THAT(v_apo, WithinAbs(A1_V_APO, V_TOL)); } - SECTION("timestep spacecraft position returns after one period") { - check_period_return(A2_A, A2_E, A2_PERIOD, - "Timestep spacecraft position"); + SECTION("apsides 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, pos_45, vel_45); + get_state(A1_A, A1_E, 0.0, 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_THAT(v_peri, WithinAbs(A1_V_PERI, V_TOL)); + REQUIRE_THAT(v_45, WithinAbs(A1_V_AT_PI4, V_TOL)); } - SECTION("timestep spacecraft velocity returns after one period") { - check_period_return(A2_A, A2_E, A2_PERIOD, - "Timestep spacecraft velocity"); + SECTION("timestep perigee radius = a*(1-e)") { + check_apsides_radius(A2_A, A2_E, 0.0, A2_A * (1.0 - A2_E), + "Timestep spacecraft perigee"); } - SECTION("timestep spacecraft true anomaly returns after one period") { - check_anomaly_return(A2_A, A2_E, A2_PERIOD, 0.0, - "Timestep spacecraft nu"); + SECTION("timestep apogee radius = a*(1+e)") { + check_apsides_radius(A2_A, A2_E, M_PI, A2_A * (1.0 + A2_E), + "Timestep spacecraft apogee"); } - destroy_simulation(sim); -} + SECTION("timestep perigee velocity > apogee velocity") { + Vec3 pos_peri, vel_peri, pos_apo, vel_apo; + get_state(A2_A, A2_E, 0.0, pos_peri, vel_peri); + get_state(A2_A, A2_E, M_PI, 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_THAT(v_peri, WithinAbs(A2_V_PERI, V_TOL)); + REQUIRE_THAT(v_apo, WithinAbs(A2_V_APO, V_TOL)); + } -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")); + // --- 2. Vis-viva --- + SECTION("vis-viva at nu = 0") { + check_vis_viva(A1_A, A1_E, 0.0, "vis-viva nu=0"); + } - Spacecraft* apsides_craft = &sim->spacecraft[0]; - CelestialBody* earth = &sim->bodies[0]; + SECTION("vis-viva at nu = pi/4") { + check_vis_viva(A1_A, A1_E, M_PI / 4.0, "vis-viva nu=pi/4"); + } - const double mu = G * earth->mass; - const double A1_A = apsides_craft->orbit.semi_major_axis; - const double A1_E = apsides_craft->orbit.eccentricity; + SECTION("vis-viva at nu = pi/2") { + check_vis_viva(A1_A, A1_E, M_PI / 2.0, "vis-viva nu=pi/2"); + } - 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" - }; + SECTION("vis-viva at nu = 3pi/4") { + check_vis_viva(A1_A, A1_E, 3.0 * M_PI / 4.0, "vis-viva nu=3pi/4"); + } - 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)); - } + SECTION("vis-viva at nu = pi") { + check_vis_viva(A1_A, A1_E, M_PI, "vis-viva nu=pi"); } - destroy_simulation(sim); -} + // --- 3. Period return (both spacecraft) --- + SECTION("apsides position returns after one period") { + check_period_return(A1_A, A1_E, A1_PERIOD, + "Apsides spacecraft position"); + } -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")); + SECTION("apsides velocity returns after one period") { + check_period_return(A1_A, A1_E, A1_PERIOD, + "Apsides spacecraft velocity"); + } - Spacecraft* timestep_craft = &sim->spacecraft[1]; - CelestialBody* earth = &sim->bodies[0]; + SECTION("apsides true anomaly returns after one period") { + check_anomaly_return(A1_A, A1_E, A1_PERIOD, 0.0, + "Apsides spacecraft nu"); + } - 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); + SECTION("timestep position returns after one period") { + check_period_return(A2_A, A2_E, A2_PERIOD, + "Timestep spacecraft position"); + } + + SECTION("timestep velocity returns after one period") { + check_period_return(A2_A, A2_E, A2_PERIOD, + "Timestep spacecraft velocity"); + } - 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("timestep true anomaly returns after one period") { + check_anomaly_return(A2_A, A2_E, A2_PERIOD, 0.0, + "Timestep spacecraft nu"); + } + // --- 4. Timestep accuracy --- 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, + propagate_and_get_state(A2_A, A2_E, 0.0, A2_PERIOD * 2.0, pos_final, vel_final); + Vec3 pos_init, vel_init; + get_state(A2_A, A2_E, 0.0, pos_init, vel_init); const double r_final = vec3_magnitude(pos_final); const double v_final = vec3_magnitude(vel_final); + const double r_init = vec3_magnitude(pos_init); + const double v_init = vec3_magnitude(vel_init); 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; @@ -278,11 +270,12 @@ SCENARIO("Accuracy across different timestep sizes", 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); + propagate_and_get_state(A2_A, A2_E, 0.0, dt, pos_final, vel_final); + Vec3 pos_init, vel_init; + get_state(A2_A, A2_E, 0.0, pos_init, vel_init); 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 expected_pos_change = vec3_magnitude(vel_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; @@ -294,18 +287,18 @@ SCENARIO("Accuracy across different timestep sizes", 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)); + REQUIRE_THAT(vel_change, WithinAbs(VEL_CHANGE_SMALL_DT, V_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); + propagate_and_get_state(A2_A, A2_E, 0.0, A2_PERIOD, pos_final, vel_final); + Vec3 pos_init, vel_init; + get_state(A2_A, A2_E, 0.0, pos_init, vel_init); 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("dt: " << A2_PERIOD << " s (1x period)"); INFO("pos_error: " << pos_error << " m"); INFO("vel_error: " << vel_error << " m/s"); @@ -314,14 +307,15 @@ SCENARIO("Accuracy across different timestep sizes", } 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, + propagate_and_get_state(A2_A, A2_E, 0.0, A2_PERIOD * 10.0, pos_final, vel_final); + Vec3 pos_init, vel_init; + get_state(A2_A, A2_E, 0.0, pos_init, vel_init); 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("dt: " << A2_PERIOD * 10.0 << " s (10x period)"); INFO("pos_error: " << pos_error << " m"); INFO("vel_error: " << vel_error << " m/s"); @@ -329,54 +323,38 @@ SCENARIO("Accuracy across different timestep sizes", REQUIRE_THAT(vel_error, WithinAbs(0.0, V_TOL)); } - destroy_simulation(sim); -} + // --- 5. Long-term stability --- + SECTION("100 periods propagation stability") { + 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; -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]; + 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; - 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 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; + } - 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; + 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); - 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; + 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)"); - // 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; + REQUIRE_THAT(anomaly_error, WithinAbs(0.0, LONG_TERM_ANG_TOL)); } - // 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); }