diff --git a/continue.md b/continue.md index c13e9f3..56b0802 100644 --- a/continue.md +++ b/continue.md @@ -27,24 +27,26 @@ - Energy functions (KE, PE, total) - Hyperbolic propagation -## 1. Structure +## Refactoring Rules + +### 1. Structure - One `SCENARIO("description")` per logical test group, with `[tag1][tag2]` annotations - Shared fixture: all constants, structs, and variables declared between `SCENARIO` opening and first `SECTION` - Precompute expected values analytically at fixture level (use `scripts/*.py` for complex simulations) -## 2. Duplication elimination +### 2. Duplication Elimination - Use lambdas that capture the fixture for repeated setup→call→assert patterns - Reuse shared structs in-place (mutate fields rather than recreating) - Single-line `SECTION`s when body is one statement: `SECTION("name") { helper(arg); }` -## 3. Assertions +### 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.) - Replace qualitative checks (`a > b`) with quantitative (`WithinAbs(expected, tol)`) - `INFO("label: " << value)` for debugging context -## 4. Precalc scripts +### 4. Precalc Scripts - For each test file, create `scripts/precalc_.py` that computes expected values. - **Check the capability matrix first** — only use sim_engine.py for implemented features. - **Always output local-frame values** (distances from parent, not global from origin). @@ -54,32 +56,46 @@ - Run with: `python3 scripts/precalc_.py` - If sim_engine.py lacks a feature, use analytical formulas instead (**but notify the user what feature was missing**) -## 5. Process per file -- **Pre-check**: Verify the test file has a TOML config in `old_tests/`. If it doesn't, skip — it's likely hardcoded. -- **Pre-check**: Check if the test is already in `tests/` (already refactored). Skip if so. -- **Pre-check**: Check the capability matrix — if the test needs SOI, maneuvers, rendezvous, etc., flag this before starting. +## Refactoring Procedure + +### Step 1: Refactor +- Verify the test file has a TOML config in `old_tests/`. If it doesn't, skip — it's likely hardcoded. +- Check if the test is already in `tests/` (already refactored). Skip if so. +- Check the capability matrix — if the test needs SOI, maneuvers, rendezvous, etc., flag this before starting. - Process **one test file at a time**. - Create `scripts/precalc_.py` and run it to get expected values. - Copy from `old_tests/` to `tests/`, rewrite using the pattern from `test_true_anomaly_roundtrip.cpp`. - Rewrite TOML configs to TOML 1.0 inline table syntax (single-line `{}`). +- Follow all rules in sections 1-4 above. + +### Step 2: Tighten Tolerances - Build and verify: `make test-build` then `./build/orbit_test '[tag]' -s`. - Run full suite: `make test`. -- **Verify no other tests broke** — the new test shouldn't affect unrelated tests. -- Note any sim_engine.py features that were missing for this test (add to capability matrix if needed). -- notify the user about any broken tests (e.g., OrbitTracker bugs) +- 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. + +### Step 3: Code Review +- Remove unused includes (only include what's actually used). +- Remove unused variables (e.g., `const double mu = G * M_sun;` if never referenced). +- Look for repeated initialization patterns — extract into helper lambdas (`make_elements`, `convert_and_recover`). +- Use `const` for all fixture data and recovered results. +- Replace C-style arrays with `std::array` where appropriate. +- Ensure consistent tolerance usage (no hardcoded `1e-4` when `ANG_TOL_COARSE` exists). +- Check for `compare_vec3` availability in `test_utilities` instead of 6 individual `REQUIRE_THAT` calls. +- Run full suite again: `make test`. - **Always ask for review** before moving to the next file. - **Only commit when asked.** ---- - ## Test Refactoring Status ### 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) ### Can Refactor Now (sim_engine.py supports all features needed) - `test_cartesian_to_elements_basic` — element conversion round-trip -- `test_cartesian_to_elements_advanced` — advanced conversion cases - `test_parabolic_orbit` — parabolic propagation via Barker's - `test_extreme_eccentricity` — high-eccentricity orbits - `test_extreme_orientation_mixed` — extreme inclinations/eccentricities @@ -100,4 +116,3 @@ ### Skip (Hardcoded / No TOML) - `test_integration` — hardcoded vector tests, no TOML config - diff --git a/src/test_utilities.cpp b/src/test_utilities.cpp index cac2ce6..bbc7a7c 100644 --- a/src/test_utilities.cpp +++ b/src/test_utilities.cpp @@ -45,6 +45,12 @@ OrbitalMetrics calculate_orbital_metrics(CelestialBody* body, CelestialBody* par return metrics; } +bool compare_vec3(Vec3 a, Vec3 b, double tolerance) { + return fabs(a.x - b.x) <= tolerance && + fabs(a.y - b.y) <= tolerance && + fabs(a.z - b.z) <= tolerance; +} + OrbitTracker* create_orbit_tracker_with_min_time(int body_index, double min_time_seconds) { OrbitTracker* tracker = (OrbitTracker*)malloc(sizeof(OrbitTracker)); tracker->body_index = body_index; diff --git a/src/test_utilities.h b/src/test_utilities.h index 4267ae9..5e20f0d 100644 --- a/src/test_utilities.h +++ b/src/test_utilities.h @@ -49,4 +49,6 @@ void destroy_orbit_tracker(OrbitTracker* tracker); int dump_simulation_state(SimulationState* sim, const char* label, char* buffer, int buffer_size); +bool compare_vec3(Vec3 a, Vec3 b, double tolerance); + #endif diff --git a/tests/test_cartesian_to_elements_advanced.cpp b/tests/test_cartesian_to_elements_advanced.cpp new file mode 100644 index 0000000..9397f0f --- /dev/null +++ b/tests/test_cartesian_to_elements_advanced.cpp @@ -0,0 +1,222 @@ +#include +#include +#include +#include +#include +#include "../src/orbital_mechanics.h" +#include "../src/test_utilities.h" + +using Catch::Matchers::WithinAbs; + +SCENARIO("Cartesian to Elements - Advanced conversion tests", + "[orbital_mechanics][cartesian][elements]") { + const double M_sun = 1.989e30; + const double A_TOL = 1e-2; + const double E_TOL = 1e-4; + const double ANG_TOL = 1e-6; + const double ANG_TOL_COARSE = 1e-4; + + auto convert_and_recover = [&](const OrbitalElements& elements) { + Vec3 pos, vel; + orbital_elements_to_cartesian(elements, M_sun, &pos, &vel); + return cartesian_to_orbital_elements(pos, vel, M_sun); + }; + + auto make_elements = [&](double a, double e, double nu, double inc, + double lon_anode, double arg_peri) { + OrbitalElements el = {}; + el.semi_major_axis = a; + el.eccentricity = e; + el.true_anomaly = nu; + el.inclination = inc; + el.longitude_of_ascending_node = lon_anode; + el.argument_of_periapsis = arg_peri; + return el; + }; + + SECTION("eccentricity spectrum: circular to highly hyperbolic") { + const double r = 1.496e11; + const double v_circular = sqrt(G * M_sun / r); + const Vec3 pos_circ = {r, 0.0, 0.0}; + const Vec3 vel_circ = {0.0, v_circular, 0.0}; + + const OrbitalElements circular = make_elements(r, 0.0, 0.0, 0.0, 0.0, 0.0); + + Vec3 converted_pos, converted_vel; + orbital_elements_to_cartesian(circular, M_sun, &converted_pos, &converted_vel); + + const OrbitalElements recovered_circ = + cartesian_to_orbital_elements(converted_pos, converted_vel, M_sun); + REQUIRE_THAT(recovered_circ.eccentricity, WithinAbs(0.0, 1e-10)); + REQUIRE_THAT(recovered_circ.semi_major_axis, WithinAbs(r, A_TOL)); + REQUIRE(compare_vec3(pos_circ, converted_pos, A_TOL)); + REQUIRE(compare_vec3(vel_circ, converted_vel, 1e-6)); + + // Near-circular (e=0.001) + const OrbitalElements near_circ = make_elements(1.496e11, 0.001, 0.5, 0.0, 0.0, 0.0); + const OrbitalElements rec_near_circ = convert_and_recover(near_circ); + REQUIRE_THAT(rec_near_circ.eccentricity, WithinAbs(0.001, 1e-6)); + REQUIRE_THAT(rec_near_circ.semi_major_axis, WithinAbs(1.496e11, A_TOL)); + + // Elliptical (e=0.5) + const OrbitalElements elliptical = make_elements(1.0e11, 0.5, 0.8, 0.0, 0.0, 0.0); + const OrbitalElements rec_elliptical = convert_and_recover(elliptical); + REQUIRE_THAT(rec_elliptical.eccentricity, WithinAbs(0.5, E_TOL)); + REQUIRE_THAT(rec_elliptical.semi_major_axis, WithinAbs(1.0e11, A_TOL)); + + // Highly elliptical (e=0.95) + const OrbitalElements high_ell = make_elements(1.0e11, 0.95, 0.1, 0.0, 0.0, 0.0); + const OrbitalElements rec_high_ell = convert_and_recover(high_ell); + REQUIRE_THAT(rec_high_ell.eccentricity, WithinAbs(0.95, 1e-3)); + REQUIRE_THAT(rec_high_ell.semi_major_axis, WithinAbs(1.0e11, A_TOL)); + + // Near-parabolic (e=0.999) + const OrbitalElements near_par = make_elements(1.0e11, 0.999, 0.05, 0.0, 0.0, 0.0); + const OrbitalElements rec_near_par = convert_and_recover(near_par); + REQUIRE_THAT(rec_near_par.eccentricity, WithinAbs(0.999, 1e-3)); + + // Parabolic (e=1.0) + OrbitalElements parabolic = {}; + parabolic.semi_latus_rectum = 1.0e11; + parabolic.eccentricity = 1.0; + parabolic.true_anomaly = 0.5; + parabolic.inclination = 0.0; + parabolic.longitude_of_ascending_node = 0.0; + parabolic.argument_of_periapsis = 0.0; + const OrbitalElements rec_parabolic = convert_and_recover(parabolic); + REQUIRE_THAT(rec_parabolic.eccentricity, WithinAbs(1.0, 1e-2)); + REQUIRE_THAT(rec_parabolic.semi_latus_rectum, WithinAbs(1.0e11, A_TOL)); + + // Hyperbolic (e=2.0) + const OrbitalElements hyper = make_elements(-1.0e11, 2.0, 0.5, 0.0, 0.0, 0.0); + const OrbitalElements rec_hyper = convert_and_recover(hyper); + REQUIRE_THAT(rec_hyper.eccentricity, WithinAbs(2.0, 1e-3)); + REQUIRE_THAT(rec_hyper.semi_major_axis, WithinAbs(-1.0e11, A_TOL)); + + // Highly hyperbolic (e=10.0) + const OrbitalElements high_hyper = make_elements(-1.0e10, 10.0, 0.8, 0.0, 0.0, 0.0); + const OrbitalElements rec_high_hyper = convert_and_recover(high_hyper); + REQUIRE_THAT(rec_high_hyper.eccentricity, WithinAbs(10.0, 1e-3)); + REQUIRE_THAT(rec_high_hyper.semi_major_axis, WithinAbs(-1.0e10, A_TOL)); + } + + SECTION("inclination: zero, polar, and retrograde") { + // Zero inclination (equatorial) + const OrbitalElements eq = make_elements(1.0e11, 0.3, 0.5, 0.0, 0.0, 0.0); + const OrbitalElements rec_eq = convert_and_recover(eq); + REQUIRE_THAT(rec_eq.inclination, WithinAbs(0.0, ANG_TOL)); + REQUIRE_THAT(rec_eq.eccentricity, WithinAbs(0.3, E_TOL)); + + // 90-degree inclination (polar) + const OrbitalElements polar = make_elements(1.0e11, 0.2, 0.6, M_PI / 2.0, 0.5, 0.3); + const OrbitalElements rec_polar = convert_and_recover(polar); + REQUIRE_THAT(rec_polar.inclination, WithinAbs(M_PI / 2.0, ANG_TOL_COARSE)); + REQUIRE_THAT(rec_polar.longitude_of_ascending_node, WithinAbs(0.5, ANG_TOL_COARSE)); + REQUIRE_THAT(rec_polar.argument_of_periapsis, WithinAbs(0.3, ANG_TOL_COARSE)); + + // 180-degree inclination (retrograde) + const OrbitalElements retro = make_elements(1.0e11, 0.2, 0.6, M_PI, 0.5, 0.3); + const OrbitalElements rec_retro = convert_and_recover(retro); + REQUIRE_THAT(rec_retro.inclination, WithinAbs(M_PI, ANG_TOL_COARSE)); + } + + SECTION("true anomaly at key orbital positions") { + struct nu_test { + double nu; + double expected_nu; + const char* label; + }; + std::vector tests = { + {0.0, 0.0, "periapsis"}, + {M_PI, M_PI, "apoapsis"}, + {M_PI / 2.0, M_PI / 2.0, "quadrature +90"}, + {-M_PI / 2.0, 3.0 * M_PI / 2.0, "quadrature -90"}, + {3.0 * M_PI / 2.0, 3.0 * M_PI / 2.0, "quadrature +270"}, + {-3.0 * M_PI / 2.0, M_PI / 2.0, "quadrature -270"}, + }; + + for (const auto& t : tests) { + const OrbitalElements elements = make_elements(1.0e11, 0.5, t.nu, 0.0, 0.0, 0.0); + const OrbitalElements recovered = convert_and_recover(elements); + INFO("Test: " << t.label << " (input nu=" << t.nu << ")"); + REQUIRE_THAT(recovered.true_anomaly, WithinAbs(t.expected_nu, ANG_TOL)); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, E_TOL)); + } + } + + SECTION("quadrature at various eccentricities") { + struct e_test { + double e; + double e_tol; + double nu_tol; + }; + std::vector tests = { + {0.9, 1e-3, 1e-5}, + {0.1, 1e-5, 1e-6}, + }; + + for (const auto& t : tests) { + const OrbitalElements elements = make_elements(1.0e11, t.e, M_PI / 2.0, 0.0, 0.0, 0.0); + const OrbitalElements recovered = convert_and_recover(elements); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(t.e, t.e_tol)); + REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, A_TOL)); + REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, t.nu_tol)); + } + } + + SECTION("large true anomaly values") { + struct large_nu_test { + double nu; + double expected_nu; + double tol; + const char* label; + }; + std::vector tests = { + {5.0, 5.0, 1e-6, "nu=5.0"}, + {-5.0, 1.28318530717958623, 1e-6, "nu=-5.0"}, + {10.0, 10.0 - 2.0 * M_PI, 1e-5, "nu=10.0"}, + }; + + for (const auto& t : tests) { + const OrbitalElements elements = make_elements(1.0e11, 0.5, t.nu, 0.0, 0.0, 0.0); + const OrbitalElements recovered = convert_and_recover(elements); + INFO("Test: " << t.label); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, E_TOL)); + REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, A_TOL)); + REQUIRE_THAT(recovered.true_anomaly, WithinAbs(t.expected_nu, t.tol)); + } + } + + SECTION("3D orientation with quadrature point") { + const OrbitalElements elements = make_elements(1.0e11, 0.5, M_PI / 2.0, + M_PI / 3.0, M_PI / 4.0, M_PI / 6.0); + const OrbitalElements recovered = convert_and_recover(elements); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, E_TOL)); + REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, A_TOL)); + REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-5)); + REQUIRE_THAT(recovered.inclination, WithinAbs(M_PI / 3.0, ANG_TOL_COARSE)); + REQUIRE_THAT(recovered.longitude_of_ascending_node, WithinAbs(M_PI / 4.0, ANG_TOL_COARSE)); + REQUIRE_THAT(recovered.argument_of_periapsis, WithinAbs(M_PI / 6.0, ANG_TOL_COARSE)); + } + + SECTION("multiple true anomaly points in sequence") { + std::array true_anomalies = {0.0, M_PI / 4.0, M_PI / 2.0, + 3.0 * M_PI / 4.0, M_PI}; + + for (double nu : true_anomalies) { + const OrbitalElements elements = make_elements(1.0e11, 0.5, nu, 0.0, 0.0, 0.0); + const OrbitalElements recovered = convert_and_recover(elements); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(0.5, E_TOL)); + REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(1.0e11, A_TOL)); + REQUIRE_THAT(recovered.true_anomaly, WithinAbs(nu, ANG_TOL)); + } + } + + SECTION("hyperbolic orbit at quadrature point") { + const OrbitalElements elements = make_elements(-1.0e11, 2.0, M_PI / 2.0, 0.0, 0.0, 0.0); + const OrbitalElements recovered = convert_and_recover(elements); + REQUIRE_THAT(recovered.eccentricity, WithinAbs(2.0, 1e-3)); + REQUIRE_THAT(recovered.semi_major_axis, WithinAbs(-1.0e11, A_TOL)); + REQUIRE_THAT(recovered.true_anomaly, WithinAbs(M_PI / 2.0, 1e-5)); + } +}