refactor: test_cartesian_to_elements_advanced + restore compare_vec3 + document procedure

- Rewrite test with SCENARIO/SECTION, helper lambdas, std::array
- Tighten tolerances: A_TOL=1e-2, E_TOL=1e-4, ANG_TOL=1e-6
- Restore compare_vec3() to test_utilities for cleaner assertions
- Document 3-step refactoring procedure in continue.md
- Fix duplicate entry in status tracking

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_<test_name>.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_<test_name>.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_<test_name>.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
-

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;

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

tests/test_cartesian_to_elements_advanced.cpp

@@ -0,0 +1,222 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <cmath>
+#include <array>
+#include <vector>
+#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<nu_test> 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<e_test> 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<large_nu_test> 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<double, 5> 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));
+    }
+}