Browse Source
- Consolidate 5 TEST_CASEs into 1 SCENARIO with 10 SECTIONs - Add REL_TOL=1e-8 to tolerance table in continue.md - Replace all qualitative checks with quantitative WithinAbs assertions - Use precalculated velocity values (from precalc script) - Replace C-style arrays with std::array - Remove unused variables and dead sections - TOML 1.0 inline table syntax config - Add precalc_extreme_eccentricity.py for expected values 40 assertions, 0 warnings, 0 FIXMEstest-refactor
4 changed files with 381 additions and 0 deletions
@ -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() |
||||
@ -0,0 +1,210 @@
|
||||
#include <catch2/catch_test_macros.hpp> |
||||
#include <catch2/matchers/catch_matchers_floating_point.hpp> |
||||
#include "../src/physics.h" |
||||
#include "../src/orbital_mechanics.h" |
||||
#include "../src/simulation.h" |
||||
#include "../src/config_loader.h" |
||||
#include <cmath> |
||||
#include <array> |
||||
|
||||
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<double, 4> 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); |
||||
} |
||||
@ -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 } |
||||
Loading…
Reference in new issue