diff --git a/continue.md b/continue.md index 90d7d11..d875884 100644 --- a/continue.md +++ b/continue.md @@ -27,6 +27,7 @@ | `R_TOL` | `1e-6` | Radius / distance magnitudes | | `V_TOL` | `1e-6` | Velocity magnitudes | | `M_TOL` | `1e-6` | Time / period values | +| `REL_TOL` | `1e-8` | Relative / percentage errors (dimensionless) | - Declare tolerance constants in the fixture (between `SCENARIO` opening and first `SECTION`) - Tighten aggressively: if observed error is `1e-8`, use `1e-6` (two orders of margin) diff --git a/scripts/precalc_extreme_eccentricity.py b/scripts/precalc_extreme_eccentricity.py new file mode 100644 index 0000000..1b5712a --- /dev/null +++ b/scripts/precalc_extreme_eccentricity.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python3 +""" +Precalculate expected values for test_extreme_eccentricity.cpp. + +Usage: + python3 scripts/precalc_extreme_eccentricity.py + +Outputs C++-style comments with precalculated values for embedding in the test. +""" + +import sys, math +sys.path.insert(0, 'scripts') +from sim_engine import orbital_to_cartesian, cartesian_to_orbital_elements, vmag, OrbitalElements, G + +# ============================================================================= +# Spacecraft 0: Highly_Elliptical (e=0.99, a=6.5e8) +# ============================================================================= +mu = G * 5.972e24 +a0 = 6.5e8 +e0 = 0.99 +nu0 = 0.0 + +elements0 = OrbitalElements(a=a0, e=e0, nu=nu0, inc=0.0, Omega=0.0, omega=0.0) +pos0, vel0 = orbital_to_cartesian(elements0, 5.972e24) + +r0 = vmag(pos0) +v0 = vmag(vel0) +expected_r_peri0 = a0 * (1.0 - e0) +expected_r_apo0 = a0 * (1.0 + e0) + +# Round-trip +elements0_rt = cartesian_to_orbital_elements(pos0, vel0, 5.972e24) + +print("# Spacecraft 0: Highly_Elliptical (e=0.99, a=6.5e8)") +print(f"# r_peri = {expected_r_peri0:.6f} m") +print(f"# r_apo = {expected_r_apo0:.6f} m") +print(f"# r = {r0:.6f} m") +print(f"# v = {v0:.6f} m/s") +print(f"# dr = {abs(r0 - expected_r_peri0):.2e} m") +print(f"# dr_apo = {abs(r0 - expected_r_apo0):.2e} m") +print(f"# e_rt = {elements0_rt.e:.15f} (error: {abs(elements0_rt.e - e0):.2e})") +print(f"# a_rt = {elements0_rt.a:.6f} m") +print() + +# Nu = pi (apoapsis) +elements0_pi = OrbitalElements(a=a0, e=e0, nu=math.pi, inc=0.0, Omega=0.0, omega=0.0) +pos0_pi, vel0_pi = orbital_to_cartesian(elements0_pi, 5.972e24) +r0_pi = vmag(pos0_pi) +v0_pi = vmag(vel0_pi) + +print(f"# At apoapsis (nu=pi):") +print(f"# r = {r0_pi:.6f} m (expected: {expected_r_apo0:.6f} m)") +print(f"# v = {v0_pi:.6f} m/s") +print(f"# dr = {abs(r0_pi - expected_r_apo0):.2e} m") +print() + +# ============================================================================= +# Spacecraft 1: Near_Parabolic (e=0.99, a=7.0e8) +# ============================================================================= +a1 = 7.0e8 +e1 = 0.99 +nu1 = 0.0 + +elements1 = OrbitalElements(a=a1, e=e1, nu=nu1, inc=0.0, Omega=0.0, omega=0.0) +pos1, vel1 = orbital_to_cartesian(elements1, 5.972e24) + +r1 = vmag(pos1) +v1 = vmag(vel1) +expected_r_peri1 = a1 * (1.0 - e1) +expected_r_apo1 = a1 * (1.0 + e1) + +# Apoapsis +elements1_pi = OrbitalElements(a=a1, e=e1, nu=math.pi, inc=0.0, Omega=0.0, omega=0.0) +pos1_pi, vel1_pi = orbital_to_cartesian(elements1_pi, 5.972e24) +r1_pi = vmag(pos1_pi) +v1_pi = vmag(vel1_pi) + +print("# Spacecraft 1: Near_Parabolic (e=0.99, a=7.0e8)") +print(f"# r_peri = {expected_r_peri1:.6f} m") +print(f"# r_apo = {expected_r_apo1:.6f} m") +print(f"# r_peri_actual = {r1:.6f} m") +print(f"# v_peri = {v1:.6f} m/s") +print(f"# r_apo_actual = {r1_pi:.6f} m") +print(f"# v_apo = {v1_pi:.6f} m/s") +print(f"# dr_peri = {abs(r1 - expected_r_peri1):.2e} m") +print(f"# dr_apo = {abs(r1_pi - expected_r_apo1):.2e} m") +print(f"# v_peri > v_apo: {v1 > v1_pi}") +print() + +# ============================================================================= +# Spacecraft 2: Slightly_Hyperbolic (e=1.05, a=-1.3e8) +# ============================================================================= +a2 = -1.3e8 +e2 = 1.05 +nu2 = 0.0 + +elements2 = OrbitalElements(a=a2, e=e2, nu=nu2, inc=0.0, Omega=0.0, omega=0.0) +pos2, vel2 = orbital_to_cartesian(elements2, 5.972e24) + +r2 = vmag(pos2) +v2 = vmag(vel2) + +escape_vel = math.sqrt(2.0 * mu / r2) +circular_vel = math.sqrt(mu / r2) +expected_v_sq = mu * (2.0 / r2 - 1.0 / a2) +expected_v = math.sqrt(expected_v_sq) + +print("# Spacecraft 2: Slightly_Hyperbolic (e=1.05, a=-1.3e8)") +print(f"# r = {r2:.6f} m") +print(f"# v = {v2:.6f} m/s") +print(f"# v_exp = {expected_v:.6f} m/s") +print(f"# v_err = {abs(v2 - expected_v):.2e} m/s") +print(f"# rel_err = {abs(v2 - expected_v) / expected_v:.2e}") +print(f"# escape_vel = {escape_vel:.6f} m/s") +print(f"# circular_vel = {circular_vel:.6f} m/s") +print(f"# a < 0: {a2 < 0}") +print() + +# ============================================================================= +# Velocity at different true anomalies for each spacecraft +# ============================================================================= +print("# Velocity magnitudes at different true anomalies:") +print("# (vis-viva: v = sqrt(mu * (2/r - 1/a)))") +print() + +for idx, (a_val, e_val, name) in enumerate([(a0, e0, "Highly_Elliptical"), + (a1, e1, "Near_Parabolic"), + (a2, e2, "Slightly_Hyperbolic")]): + print(f"# {name} (a={a_val:.2e}, e={e_val:.2f}):") + for nu in [0.0, math.pi/2.0, math.pi, 3.0*math.pi/2.0]: + if e_val > 1.0: + max_nu = math.acos(-1.0 / e_val) + if abs(nu) >= max_nu: + print(f"# nu={nu:.4f} rad: SKIPPED (hyperbolic limit +/- {max_nu:.4f})") + continue + elem = OrbitalElements(a=a_val, e=e_val, nu=nu, inc=0.0, Omega=0.0, omega=0.0) + p, v = orbital_to_cartesian(elem, 5.972e24) + r = vmag(p) + v_mag = vmag(v) + v_exp = math.sqrt(mu * (2.0/r - 1.0/a_val)) + rel_err = abs(v_mag - v_exp) / v_exp + print(f"# nu={nu:.4f} rad: v={v_mag:.6f} m/s, v_exp={v_exp:.6f} m/s, rel_err={rel_err:.2e}") + print() diff --git a/tests/test_extreme_eccentricity.cpp b/tests/test_extreme_eccentricity.cpp new file mode 100644 index 0000000..e5ace25 --- /dev/null +++ b/tests/test_extreme_eccentricity.cpp @@ -0,0 +1,210 @@ +#include +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include +#include + +using Catch::Matchers::WithinAbs; + +SCENARIO("Extreme eccentricity orbital conversions and vis-viva accuracy", + "[extreme][eccentricity][high]") { + const double TIME_STEP = 60.0; + const double parent_mass = 5.972e24; + const double mu = G * parent_mass; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_extreme_eccentricity.toml")); + + Spacecraft* high_e = &sim->spacecraft[0]; + Spacecraft* near_parabolic = &sim->spacecraft[1]; + Spacecraft* hyperbolic = &sim->spacecraft[2]; + + // Tolerances + const double R_TOL = 1e-6; + const double V_TOL = 1e-6; + const double E_TOL = 1e-12; + const double REL_TOL = 1e-8; + + // Precomputed analytical values for spacecraft 0 (a=6.5e8, e=0.99) + const double a0 = high_e->orbit.semi_major_axis; + const double e0 = high_e->orbit.eccentricity; + const double expected_r_peri0 = a0 * (1.0 - e0); // 6.5e6 + const double expected_r_apo0 = a0 * (1.0 + e0); // 1.2935e9 + + // Precomputed analytical values for spacecraft 1 (a=7.0e8, e=0.99) + const double a1 = near_parabolic->orbit.semi_major_axis; + const double e1 = near_parabolic->orbit.eccentricity; + const double expected_r_peri1 = a1 * (1.0 - e1); // 7.0e6 + const double expected_r_apo1 = a1 * (1.0 + e1); // 1.393e9 + + // Precomputed analytical values for spacecraft 2 (e=1.05) + const double e2 = hyperbolic->orbit.eccentricity; + const double max_nu_hyperbolic = acos(-1.0 / e2); // ~2.8317 rad + + // Helper: convert elements to cartesian and check vis-viva consistency + auto check_visviva = [&](const OrbitalElements& orbit, double r, double v) { + double expected_v_sq = mu * (2.0 / r - 1.0 / orbit.semi_major_axis); + REQUIRE(expected_v_sq > 0.0); + const double expected_v = sqrt(expected_v_sq); + const double rel_err = fabs(v - expected_v) / expected_v; + INFO("v=" << v << " m/s, v_exp=" << expected_v << " m/s, rel_err=" << rel_err); + REQUIRE_THAT(rel_err, WithinAbs(0.0, REL_TOL)); + }; + + // Helper: convert elements to cartesian at given true anomaly + auto convert_at_nu = [&](Spacecraft* craft, double nu) { + craft->orbit.true_anomaly = nu; + orbital_elements_to_cartesian(craft->orbit, parent_mass, &craft->local_position, &craft->local_velocity); + }; + + // Helper: round-trip check + auto roundtrip = [&](double a, double e, double nu) { + OrbitalElements elements = {}; + elements.semi_major_axis = a; + elements.eccentricity = e; + elements.true_anomaly = nu; + Vec3 pos, vel; + orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); + OrbitalElements recovered = cartesian_to_orbital_elements(pos, vel, parent_mass); + return recovered; + }; + + SECTION("highly elliptical: periapsis radius = a*(1-e)") { + convert_at_nu(high_e, 0.0); + const double r = vec3_magnitude(high_e->local_position); + const double v = vec3_magnitude(high_e->local_velocity); + + INFO("r=" << r << " m, expected=" << expected_r_peri0 << " m"); + INFO("v=" << v << " m/s"); + + REQUIRE_THAT(r, WithinAbs(expected_r_peri0, R_TOL)); + REQUIRE_THAT(v, WithinAbs(11046.701562, V_TOL)); + check_visviva(high_e->orbit, r, v); + + // Round-trip eccentricity accuracy + const OrbitalElements recovered = roundtrip(a0, e0, 0.0); + INFO("e_recovered=" << recovered.eccentricity << ", error=" << fabs(recovered.eccentricity - e0)); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(e0, E_TOL)); + } + + SECTION("highly elliptical: apoapsis radius = a*(1+e)") { + convert_at_nu(high_e, M_PI); + const double r = vec3_magnitude(high_e->local_position); + const double v = vec3_magnitude(high_e->local_velocity); + + INFO("r=" << r << " m, expected=" << expected_r_apo0 << " m"); + INFO("v=" << v << " m/s"); + + REQUIRE_THAT(r, WithinAbs(expected_r_apo0, R_TOL)); + REQUIRE_THAT(v, WithinAbs(55.511063, V_TOL)); + check_visviva(high_e->orbit, r, v); + } + + SECTION("near-parabolic: periapsis radius") { + convert_at_nu(near_parabolic, 0.0); + const double r = vec3_magnitude(near_parabolic->local_position); + const double v = vec3_magnitude(near_parabolic->local_velocity); + + INFO("r=" << r << " m, expected=" << expected_r_peri1 << " m"); + INFO("v=" << v << " m/s"); + + REQUIRE_THAT(r, WithinAbs(expected_r_peri1, R_TOL)); + REQUIRE_THAT(v, WithinAbs(10644.867979, V_TOL)); + check_visviva(near_parabolic->orbit, r, v); + } + + SECTION("near-parabolic: apoapsis radius") { + convert_at_nu(near_parabolic, M_PI); + const double r = vec3_magnitude(near_parabolic->local_position); + const double v = vec3_magnitude(near_parabolic->local_velocity); + + INFO("r=" << r << " m, expected=" << expected_r_apo1 << " m"); + INFO("v=" << v << " m/s"); + + REQUIRE_THAT(r, WithinAbs(expected_r_apo1, R_TOL)); + REQUIRE_THAT(v, WithinAbs(53.491799, V_TOL)); + check_visviva(near_parabolic->orbit, r, v); + } + + SECTION("near-parabolic: velocity at periapsis and apoapsis") { + near_parabolic->orbit.true_anomaly = 0.0; + Vec3 dummy, vel_peri; + orbital_elements_to_cartesian(near_parabolic->orbit, parent_mass, &dummy, &vel_peri); + const double v_peri = vec3_magnitude(vel_peri); + + near_parabolic->orbit.true_anomaly = M_PI; + Vec3 vel_apo; + orbital_elements_to_cartesian(near_parabolic->orbit, parent_mass, &dummy, &vel_apo); + const double v_apo = vec3_magnitude(vel_apo); + + INFO("v_peri=" << v_peri << " m/s, v_apo=" << v_apo << " m/s"); + REQUIRE_THAT(v_peri, WithinAbs(10644.867979, V_TOL)); + REQUIRE_THAT(v_apo, WithinAbs(53.491799, V_TOL)); + } + + SECTION("hyperbolic: velocity matches vis-viva") { + convert_at_nu(hyperbolic, 0.0); + const double r = vec3_magnitude(hyperbolic->local_position); + const double v = vec3_magnitude(hyperbolic->local_velocity); + + INFO("r=" << r << " m"); + INFO("v=" << v << " m/s"); + + REQUIRE_THAT(v, WithinAbs(11211.998050, V_TOL)); + } + + SECTION("hyperbolic: true anomaly limits") { + INFO("max_nu=" << max_nu_hyperbolic << " rad (±" << max_nu_hyperbolic * 180.0 / M_PI << "°)"); + + // pi and 3pi/2 should be outside hyperbolic range + const double pi = M_PI; + const double three_pi_half = 3.0 * M_PI / 2.0; + INFO("pi=" << pi << " rad, exceeds limit: " << (fabs(pi) >= max_nu_hyperbolic)); + INFO("3pi/2=" << three_pi_half << " rad, exceeds limit: " << (fabs(three_pi_half) >= max_nu_hyperbolic)); + REQUIRE(fabs(pi) >= max_nu_hyperbolic); + REQUIRE(fabs(three_pi_half) >= max_nu_hyperbolic); + } + + SECTION("vis-viva accuracy at multiple true anomalies") { + const std::array true_anomalies = {0.0, M_PI / 2.0, M_PI, 3.0 * M_PI / 2.0}; + + for (int i = 0; i < sim->craft_count; i++) { + Spacecraft* craft = &sim->spacecraft[i]; + const double a = craft->orbit.semi_major_axis; + const double e = craft->orbit.eccentricity; + + INFO("Spacecraft " << i << ": e=" << e << ", a=" << a); + + for (int j = 0; j < 4; j++) { + double nu = true_anomalies[j]; + + if (e > 1.0) { + if (fabs(nu) >= max_nu_hyperbolic) { + INFO(" nu=" << nu << " rad: SKIPPED (exceeds hyperbolic limit)"); + continue; + } + } + + craft->orbit.true_anomaly = nu; + Vec3 pos, vel; + orbital_elements_to_cartesian(craft->orbit, parent_mass, &pos, &vel); + + const double r = vec3_magnitude(pos); + const double v = vec3_magnitude(vel); + + double expected_v_sq = mu * (2.0 / r - 1.0 / a); + if (expected_v_sq > 0.0) { + const double expected_v = sqrt(expected_v_sq); + const double rel_err = fabs(v - expected_v) / expected_v; + INFO(" nu=" << nu << " rad: v=" << v << " m/s, rel_err=" << rel_err); + REQUIRE_THAT(rel_err, WithinAbs(0.0, REL_TOL)); + } + } + } + } + + destroy_simulation(sim); +} diff --git a/tests/test_extreme_eccentricity.toml b/tests/test_extreme_eccentricity.toml new file mode 100644 index 0000000..1e8a58a --- /dev/null +++ b/tests/test_extreme_eccentricity.toml @@ -0,0 +1,27 @@ +# Test Configuration: Extreme Eccentricity Orbits + +[[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 = "Highly_Elliptical" +mass = 1000.0 +parent_index = 0 +orbit = { semi_major_axis = 6.5e8, eccentricity = 0.99, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[spacecraft]] +name = "Near_Parabolic" +mass = 1000.0 +parent_index = 0 +orbit = { semi_major_axis = 7.0e8, eccentricity = 0.99, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 } + +[[spacecraft]] +name = "Slightly_Hyperbolic" +mass = 1000.0 +parent_index = 0 +orbit = { semi_major_axis = -1.3e8, eccentricity = 1.05, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }