diff --git a/continue.md b/continue.md index 56b0802..d29d45e 100644 --- a/continue.md +++ b/continue.md @@ -42,7 +42,21 @@ ### 3. Assertions - `using Catch::Matchers::WithinAbs;` after includes - `REQUIRE_THAT(value, WithinAbs(expected, tolerance))` — never `Approx()` -- Tolerances based on actual observed errors, tightened aggressively (1e-12 for angles, 1e-6 for meters, etc.) +- **Always use named tolerance constants** — never hardcode raw numbers in `WithinAbs()`. + +#### Tolerance Reference +| Constant | Value | Use for | +|----------|-------|---------| +| `A_TOL` | `1e-6` | Semi-major axis (meters) | +| `E_TOL` | `1e-12` | Eccentricity | +| `ANG_TOL` | `1e-12` | Angles (true anomaly, inclination, Ω, ω) — round-trip precision | +| `ANG_TOL_COARSE` | `1e-4` | Angles — degenerate cases (polar/equatorial orbits, near-180° inclination) | +| `R_TOL` | `1e-6` | Radius / distance magnitudes | +| `V_TOL` | `1e-6` | Velocity magnitudes | +| `M_TOL` | `1e-6` | Time / period values | + +- 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) - Replace qualitative checks (`a > b`) with quantitative (`WithinAbs(expected, tol)`) - `INFO("label: " << value)` for debugging context @@ -72,9 +86,9 @@ - Build and verify: `make test-build` then `./build/orbit_test '[tag]' -s`. - Run full suite: `make test`. - Review every tolerance against actual observed errors from `-s` output. -- Tighten aggressively: round-trip conversions should use `1e-2` or better for meters, `1e-6` for angles. -- Replace hardcoded loose tolerances with named constants (`A_TOL`, `E_TOL`, `ANG_TOL`, `ANG_TOL_COARSE`). -- Ensure no `1e3`, `1e5`, `1e6` for position comparisons in pure conversion tests. +- If a constant's margin is too loose (e.g., `1e-6` when error is `1e-10`), tighten to `1e-8`. +- **If a test fails due to a tolerance being too tight, report the observed error to the user and ask whether to loosen the constant or investigate the root cause. Never silently widen a tolerance.** +- Refer to the tolerance reference table in Section 3 for constant names. ### Step 3: Code Review - Remove unused includes (only include what's actually used). @@ -93,9 +107,9 @@ ### Completed - `test_barkers_equation` ✅ — Barker's equation unit tests + parabolic propagation - `test_cartesian_to_elements_advanced` ✅ — Advanced conversion tests (eccentricity spectrum, inclination, true anomaly, 3D orientation) +- `test_cartesian_to_elements_basic` ✅ — Element round-trip conversion (semi-major axis, eccentricity, true anomaly, inclination, radius, velocity) ### Can Refactor Now (sim_engine.py supports all features needed) -- `test_cartesian_to_elements_basic` — element conversion round-trip - `test_parabolic_orbit` — parabolic propagation via Barker's - `test_extreme_eccentricity` — high-eccentricity orbits - `test_extreme_orientation_mixed` — extreme inclinations/eccentricities diff --git a/scripts/precalc_cartesian_to_elements_basic.py b/scripts/precalc_cartesian_to_elements_basic.py new file mode 100644 index 0000000..b821fd8 --- /dev/null +++ b/scripts/precalc_cartesian_to_elements_basic.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 +""" +Precalculate expected values for test_cartesian_to_elements_basic.cpp. + +Usage: + python3 scripts/precalc_cartesian_to_elements_basic.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 + +# ============================================================================= +# Test configuration: moderate eccentricity, zero inclination +# ============================================================================= +mu = G * 5.972e24 +a = 1.5e7 +e = 0.5 +nu = 0.0 +inc = 0.0 +Omega = 0.0 +omega = 0.0 + +elements = OrbitalElements(a=a, e=e, nu=nu, inc=inc, Omega=Omega, omega=omega) +pos, vel = orbital_to_cartesian(elements, 5.972e24) + +r = vmag(pos) +v = vmag(vel) + +# Round-trip: convert back +elements_rt = cartesian_to_orbital_elements(pos, vel, 5.972e24) + +print("# Test: test_cartesian_to_elements_basic") +print(f"#") +print(f"# Original elements:") +print(f"# a = {a:.6f}") +print(f"# e = {e:.6f}") +print(f"# nu = {nu:.6f}") +print(f"# inc = {inc:.6f}") +print(f"# Omega = {Omega:.6f}") +print(f"# omega = {omega:.6f}") +print(f"#") +print(f"# State vectors from elements:") +print(f"# pos = ({pos[0]:.6f}, {pos[1]:.6f}, {pos[2]:.6f}) m") +print(f"# vel = ({vel[0]:.6f}, {vel[1]:.6f}, {vel[2]:.6f}) m/s") +print(f"# r = {r:.6f} m") +print(f"# v = {v:.6f} m/s") +print(f"#") +print(f"# Round-trip recovered elements:") +print(f"# a = {elements_rt.a:.15f}") +print(f"# e = {elements_rt.e:.15f}") +print(f"# nu = {elements_rt.nu:.15f}") +print(f"# inc = {elements_rt.inc:.15f}") +print(f"# Omega = {elements_rt.Omega:.15f}") +print(f"# omega = {elements_rt.omega:.15f}") +print(f"#") +print(f"# Errors:") +print(f"# da = {abs(elements_rt.a - a):.2e}") +print(f"# de = {abs(elements_rt.e - e):.2e}") +print(f"# dnu = {abs(elements_rt.nu - nu):.2e}") +print(f"# dinc = {abs(elements_rt.inc - inc):.2e}") +print(f"# dOmega = {abs(elements_rt.Omega - Omega):.2e}") +print(f"# domega = {abs(elements_rt.omega - omega):.2e}") +print(f"# dr = {abs(r - r):.2e} (trivial)") +print(f"# dv = {abs(v - v):.2e} (trivial)") + +# Re-convert recovered elements back to state vectors +pos2, vel2 = orbital_to_cartesian(elements_rt, 5.972e24) +r2 = vmag(pos2) +v2 = vmag(vel2) + +print(f"#") +print(f"# Reconstructed from recovered elements:") +print(f"# pos = ({pos2[0]:.6f}, {pos2[1]:.6f}, {pos2[2]:.6f}) m") +print(f"# vel = ({vel2[0]:.6f}, {vel2[1]:.6f}, {vel2[2]:.6f}) m/s") +print(f"# r = {r2:.6f} m") +print(f"# v = {v2:.6f} m/s") +print(f"# dr = {abs(r2 - r):.2e} m") +print(f"# dv = {abs(v2 - v):.2e} m/s") diff --git a/tests/test_cartesian_to_elements_basic.cpp b/tests/test_cartesian_to_elements_basic.cpp new file mode 100644 index 0000000..d42b47a --- /dev/null +++ b/tests/test_cartesian_to_elements_basic.cpp @@ -0,0 +1,89 @@ +#include +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include + +using Catch::Matchers::WithinAbs; + +SCENARIO("Cartesian ↔ orbital elements round-trip conversion", + "[cartesian][elements][roundtrip]") { + const double TIME_STEP = 60.0; + const double parent_mass = 5.972e24; + + const double A_TOL = 1e-6; + const double E_TOL = 1e-12; + const double ANG_TOL = 1e-12; + const double R_TOL = 1e-6; + const double V_TOL = 1e-6; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_cartesian_to_elements_basic.toml")); + Spacecraft* craft = &sim->spacecraft[0]; + + const OrbitalElements& orig = craft->orbit; + + // Convert elements → state vectors + Vec3 pos, vel; + orbital_elements_to_cartesian(orig, parent_mass, &pos, &vel); + + const double expected_r = vec3_magnitude(pos); // 7500000.0 + const double expected_v = vec3_magnitude(vel); // 8928.484709... + + // Round-trip: state vectors → elements + const OrbitalElements recovered = cartesian_to_orbital_elements(pos, vel, parent_mass); + + // Re-convert recovered elements → state vectors + Vec3 pos2, vel2; + orbital_elements_to_cartesian(recovered, parent_mass, &pos2, &vel2); + const double recovered_r = vec3_magnitude(pos2); + const double recovered_v = vec3_magnitude(vel2); + + SECTION("elements round-trip: semi-major axis") { + const double da = fabs(recovered.semi_major_axis - orig.semi_major_axis); + INFO("Original a: " << orig.semi_major_axis); + INFO("Recovered a: " << recovered.semi_major_axis); + INFO("Error: " << da << " m"); + REQUIRE_THAT(da, WithinAbs(0.0, A_TOL)); + } + + SECTION("elements round-trip: eccentricity") { + const double de = fabs(recovered.eccentricity - orig.eccentricity); + INFO("Original e: " << orig.eccentricity); + INFO("Recovered e: " << recovered.eccentricity); + INFO("Error: " << de); + REQUIRE_THAT(de, WithinAbs(0.0, E_TOL)); + } + + SECTION("elements round-trip: true anomaly") { + const double dnu = fabs(recovered.true_anomaly - orig.true_anomaly); + INFO("Original nu: " << orig.true_anomaly); + INFO("Recovered nu: " << recovered.true_anomaly); + INFO("Error: " << dnu); + REQUIRE_THAT(dnu, WithinAbs(0.0, ANG_TOL)); + } + + SECTION("elements round-trip: inclination") { + const double dinc = fabs(recovered.inclination - orig.inclination); + INFO("Original inc: " << orig.inclination); + INFO("Recovered inc: " << recovered.inclination); + INFO("Error: " << dinc); + REQUIRE_THAT(dinc, WithinAbs(0.0, ANG_TOL)); + } + + SECTION("radius preservation") { + INFO("Original r: " << expected_r); + INFO("Recovered r: " << recovered_r); + REQUIRE_THAT(recovered_r, WithinAbs(expected_r, R_TOL)); + } + + SECTION("velocity magnitude preservation") { + INFO("Original v: " << expected_v); + INFO("Recovered v: " << recovered_v); + REQUIRE_THAT(recovered_v, WithinAbs(expected_v, V_TOL)); + } + + destroy_simulation(sim); +} diff --git a/tests/test_cartesian_to_elements_basic.toml b/tests/test_cartesian_to_elements_basic.toml new file mode 100644 index 0000000..4e66204 --- /dev/null +++ b/tests/test_cartesian_to_elements_basic.toml @@ -0,0 +1,15 @@ +# Basic elliptical orbit: a=15000km, e=0.5, zero inclination + +[[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 = "Test_Spacecraft" +mass = 1000.0 +parent_index = 0 +orbit = { semi_major_axis = 1.5e7, eccentricity = 0.5, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }