Add Molniya orbit test cases and document config loader bug

- Create test configuration for Molniya orbits with Earth as root body - Implement test suite for highly inclined orbits: * Position verification at multiple true anomalies * Orbital period verification * Generic inclined orbit test * Inclination parameter preservation - Document critical bug in spacecraft initialization: * Config loader incorrectly adds parent radius to semi_major_axis * Affects all spacecraft using semi_major_axis directly * Causes significant position errors (1-11M meters) * Molniya tests fail due to this bug, not test code
5 months ago · cfb2c92bd4
3 changed files with 486 additions and 0 deletions
--- a/docs/planning/molniya-orbit-test-plan.md
+++ b/docs/planning/molniya-orbit-test-plan.md
@ -0,0 +1,248 @@
 # Plan: Add Molniya Orbit Test Case
 ## Overview
 Add test cases for highly inclined orbits, specifically focusing on Molniya orbits and generic inclined orbit behavior.
 ## Background - Molniya Orbit Properties
 From Wikipedia research, Molniya orbits are designed for high-latitude coverage:
 - **Orbital Period**: ~718 minutes (~12 hours, half a sidereal day)
 - **Eccentricity**: ~0.74
 - **Inclination**: 63.4° (critical value that prevents perigee precession)
 - **Argument of Perigee**: 270° (apogee at northernmost point)
 - **Perigee altitude**: ~600 km
 - **Apogee altitude**: ~39,700 km
 - **Semi-major axis**: ~26,600 km
 ## Current Codebase Status
 ### Earth as Root Body
 ✅ **SUPPORTED**
 - Validator checks: `parent_index < body_index or -1`
 - Earth can be root if it's first body (index 0)
 - No special validation prevents this
 ### 3D Orientation Support
 ⚠️ **DEFINED BUT NOT APPLIED**
 - `inclination`, `longitude_of_ascending_node`, `argument_of_periapsis` exist in `OrbitalElements` struct
 - `orbital_elements_to_cartesian()` in `src/orbital_mechanics.cpp` only produces 2D orbits (x, y, z=0)
 - Documented as "deferred implementation" in technical_reference.md line 114
 - Config parser supports loading these parameters (src/config_loader.cpp)
 - Test approach: Tests will expect 3D behavior and fail with `[!mayfail]` tag
 ## Implementation Plan
 ### Step 1: Create Test Configuration File
 **File**: `tests/test_molniya.toml`
 ```toml
 # Test Configuration: Molniya Orbit
 # Earth as root body with highly elliptical, highly inclined satellite orbit
 [[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 = "Molniya_Satellite"
 mass = 1000.0
 parent_index = 0
 orbit = {
    semi_major_axis = 26540000.0,
    eccentricity = 0.74,
    true_anomaly = 0.0,
    inclination = 1.107,           # 63.4° in radians
    longitude_of_ascending_node = 0.0,
    argument_of_periapsis = 4.71    # 270° in radians
 }
 ```
 **Key points**:
 - Earth is root body (index 0, parent_index = -1)
 - Spacecraft uses Molniya orbital parameters
 - 3D orientation parameters included even though not applied yet
 ### Step 2: Create Test File
 **File**: `tests/test_inclined_orbits.cpp`
 **Test Cases**:
 1. **Molniya Position Verification** (tagged `[!mayfail]`)
   - Test satellite position at true_anomaly = 0 (perigee)
   - Test satellite position at true_anomaly = π/2
   - Test satellite position at true_anomaly = π (apogee)
   - Test satellite position at true_anomaly = 3π/2
   - Verify radii match polar equation: r = a(1-e²)/(1+e·cos(ν))
   - Check z-coordinate ≠ 0 for inclined orbits (will fail until 3D implemented)
   - Tolerance: ±10 km (10,000 meters)
 2. **Molniya Orbital Period** (tagged `[!mayfail]`)
   - Verify period matches theoretical calculation from Kepler's 3rd law
   - Theoretical period: T = 2π√(a³/μ)
   - Use `OrbitTracker` with minimum time to prevent false completion
   - Expected: ~11.96 hours
   - Tolerance: ±600 seconds (10 minutes)
 3. **Generic Inclined Orbit** (tagged `[!mayfail]`)
   - Test moderate inclination (45°) with e=0.5, a=10,000 km
   - Verify z-coordinate is non-zero (will fail until 3D implemented)
   - Verify position magnitude matches orbital radius
   - Simpler case to test 3D orientation without Molniya complexity
 4. **Inclination Parameter Preservation**
   - Verify config loader correctly reads and preserves inclination value
   - This should PASS (just checking config loading, not physics)
   - Ensures 3D parameters are being stored correctly
 **Constants**:
 - `POSITION_TOLERANCE_METERS = 10000.0` (±10 km per user request)
 - `PERIOD_TOLERANCE_SECONDS = 600.0` (±10 minutes)
 ### Step 3: Update Build System
 Add test file to existing test build (Catch2 automatically includes all `test_*.cpp` files)
 ### Step 4: Expected Behavior
 #### Before 3D Implementation
 - All tests with `[!mayfail]` tag will run but failure is acceptable
 - Position tests: Pass radius check, fail on z-coordinate check
 - Period test: May pass (period independent of inclination for same orbital energy)
 - Generic inclined test: Fail on z-coordinate check
 - Config preservation test: **PASS** (just checks config loading)
 #### After 3D Implementation
 - All position tests should **PASS**
 - Period test should **PASS**
 - Generic inclined test should **PASS**
 - Config preservation test should **PASS**
 - Remove `[!mayfail]` tags from passing tests
 ## Implementation Findings
 ### Bug Discovered in Spacecraft Initialization
 **Location**: `src/config_loader.cpp` lines 257-259
 **Code**:
 ```cpp
 if (!is_parabolic && craft->parent_index >= 0 && craft->parent_index < sim->body_count) {
    CelestialBody* parent = &sim->bodies[craft->parent_index];
    craft->orbit.semi_major_axis += parent->radius;  // BUG: Always adds parent radius
 }
 ```
 **Problem**:
 - For **spacecraft**, code unconditionally adds parent radius to `semi_major_axis`
 - For **bodies**, this logic does NOT exist (bodies use semi_major_axis as-is)
 - This creates an inconsistency between body and spacecraft orbital definitions
 - The logic was intended for `altitude` parameter convenience but is applied even when `semi_major_axis` is specified directly
 **Impact on Molniya Tests**:
 - Config specifies: `semi_major_axis = 26,540,000 m`
 - Actual used: `26,540,000 + 6,371,000 = 32,911,000 m`
 - Position errors: 1.7M m at perigee, 11.1M m at apogee
 - Orbital period: ~16.5 hours instead of expected ~12 hours
 - All radius-based tests fail with errors far exceeding ±10 km tolerance
 **Root Cause**:
 The config loader has these paths (simplified):
 ```cpp
 // For bodies:
 if (semi_major_axis.type == TOML_FP64) {
    body->orbit.semi_major_axis = semi_major.u.fp64;  // Used directly
 } else if (altitude.type == TOML_FP64) {
    body->orbit.semi_major_axis = altitude.u.fp64;  // Stored as altitude
 }
 // For spacecraft:
 if (semi_major_axis.type == TOML_FP64) {
    craft->orbit.semi_major_axis = semi_major.u.fp64;  // Stored as SMA
 } else if (altitude.type == TOML_FP64) {
    craft->orbit.semi_major_axis = altitude.u.fp64;  // Stored as altitude
 }
 // AFTER loading all spacecraft (line 253-260):
 for (int i = 0; i < sim->craft_count; i++) {
    Spacecraft* craft = &sim->spacecraft[i];
    bool is_parabolic = (fabs(craft->orbit.eccentricity - 1.0) < 0.005);
    if (!is_parabolic && craft->parent_index >= 0) {
        CelestialBody* parent = &sim->bodies[craft->parent_index];
        craft->orbit.semi_major_axis += parent->radius;  // ALWAYS adds!
    }
 }
 ```
 **Required Fix**:
 The logic should only add parent radius when `altitude` was specified:
 - Option A: Track which parameter was used during parsing
 - Option B: Check if semi_major_axis < parent radius (altitude would be small)
 - Option C: Only apply this transformation for altitude-based initialization
 ### Additional Issues Found
 #### OrbitTracker Not Completing
 - **Symptom**: Molniya period test never completes orbit (tracker->orbit_completed = false)
 - **Possible causes**:
  1. Minimum time threshold too restrictive (currently 0.01 days = 14.4 minutes)
  2. Quadrant transition logic failing for highly elliptical orbits
  3. Angle comparison tolerance too tight for large orbital variations
 - **Status**: Not investigated yet, secondary priority to config bug
 #### Test File Name Mismatch
 - Original plan: `tests/test_molniya.toml`
 - Actual created: `tests/test_inclined_orbits.toml` (matches test file name)
 - Reason: User preference to have config filename match test filename
 - **Status**: ✅ RESOLVED (updated config in final implementation)
 ## Test Verification Formula
 ### Orbital Radius at True Anomaly
 ```
 r = a(1-e²) / (1 + e·cos(ν))
 ```
 Where:
 - a = semi-major_axis
 - e = eccentricity
 - ν = true_anomaly
 ### Orbital Period (Kepler's 3rd Law)
 ```
 T = 2π√(a³/μ)
 ```
 Where:
 - μ = G·M (standard gravitational parameter)
 - M = parent body mass
 ## Questions Resolved
 1. **Tolerance**: ±10 km (10,000 meters) - chosen as moderate tolerance for RK4 integration
 2. **Test Scope**: Both Molniya + generic inclined orbit tests for completeness
 3. **3D Status**: Tests expect 3D and fail with `[!mayfail]` tag, documenting expected behavior
 ## Next Steps After This Plan
 1. Implement 3D orientation in `src/orbital_mechanics.cpp`
   - Apply rotation matrices for inclination, RAAN, and argument of periapsis
   - Update `orbital_elements_to_cartesian()` to include 3D transformations
 2. Remove `[!mayfail]` tags from tests that now pass
 3. Add additional 3D-specific tests:
   - Verify correct orientation of orbital plane
   - Test ascending/descending node crossings
   - Verify argument of periapsis positioning
 ## Files to Create
 - `tests/test_molniya.toml` (new test config)
 - `tests/test_inclined_orbits.cpp` (new test file)
 ## Files to Potentially Modify (future 3D implementation)
 - `src/orbital_mechanics.cpp` - add 3D rotation logic
 - `src/orbital_mechanics.h` - potentially expose rotation helper functions
--- a/tests/test_inclined_orbits.cpp
+++ b/tests/test_inclined_orbits.cpp
@ -0,0 +1,203 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 #include "../src/physics.h"
 #include "../src/simulation.h"
 #include "../src/config_loader.h"
 #include "../src/test_utilities.h"
 #include <cmath>
 const double POSITION_TOLERANCE_METERS = 10000.0;
 const double PERIOD_TOLERANCE_SECONDS = 600.0;
 TEST_CASE("Molniya orbit - position verification at multiple true anomalies", "[inclined][molniya][!mayfail]") {
    const double TIME_STEP = 60.0;
    const double SECONDS_PER_DAY = 86400.0;
    const double SEMI_MAJOR_AXIS = 26540000.0;
    const double ECCENTRICITY = 0.74;
    const double EARTH_MASS = 5.972e24;
    const double MU = G * EARTH_MASS;
    SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP);
    REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml"));
    Spacecraft* molniya = &sim->spacecraft[0];
    CelestialBody* earth = &sim->bodies[0];
    SECTION("Position at perigee (true_anomaly = 0)") {
        double expected_radius = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY);
        double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position));
        double radius_error = fabs(actual_radius - expected_radius);
        INFO("Expected radius at perigee: " << expected_radius << " m");
        INFO("Actual radius: " << actual_radius << " m");
        INFO("Error: " << radius_error << " m");
        REQUIRE(radius_error < POSITION_TOLERANCE_METERS);
        CHECK(molniya->global_position.z != 0.0);
        INFO("Z-coordinate should be non-zero for inclined orbit (currently deferred)");
    }
    SECTION("Position at true_anomaly = π/2 (90°)") {
        molniya->orbit.true_anomaly = M_PI / 2.0;
        initialize_orbital_objects(sim);
        double expected_radius = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY * ECCENTRICITY) / (1.0 + ECCENTRICITY * cos(M_PI / 2.0));
        double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position));
        double radius_error = fabs(actual_radius - expected_radius);
        INFO("Expected radius at ν=π/2: " << expected_radius << " m");
        INFO("Actual radius: " << actual_radius << " m");
        INFO("Error: " << radius_error << " m");
        REQUIRE(radius_error < POSITION_TOLERANCE_METERS);
        CHECK(molniya->global_position.z != 0.0);
    }
    SECTION("Position at apogee (true_anomaly = π)") {
        molniya->orbit.true_anomaly = M_PI;
        initialize_orbital_objects(sim);
        double expected_radius = SEMI_MAJOR_AXIS * (1.0 + ECCENTRICITY);
        double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position));
        double radius_error = fabs(actual_radius - expected_radius);
        INFO("Expected radius at apogee: " << expected_radius << " m");
        INFO("Actual radius: " << actual_radius << " m");
        INFO("Error: " << radius_error << " m");
        REQUIRE(radius_error < POSITION_TOLERANCE_METERS);
        CHECK(molniya->global_position.z != 0.0);
        INFO("At apogee, satellite should be at northernmost point (max z)");
    }
    SECTION("Position at true_anomaly = 3π/2 (270°)") {
        molniya->orbit.true_anomaly = 3.0 * M_PI / 2.0;
        initialize_orbital_objects(sim);
        double expected_radius = SEMI_MAJOR_AXIS * (1.0 - ECCENTRICITY * ECCENTRICITY) / (1.0 + ECCENTRICITY * cos(3.0 * M_PI / 2.0));
        double actual_radius = vec3_magnitude(vec3_sub(molniya->global_position, earth->global_position));
        double radius_error = fabs(actual_radius - expected_radius);
        INFO("Expected radius at ν=3π/2: " << expected_radius << " m");
        INFO("Actual radius: " << actual_radius << " m");
        INFO("Error: " << radius_error << " m");
        REQUIRE(radius_error < POSITION_TOLERANCE_METERS);
        CHECK(molniya->global_position.z != 0.0);
        INFO("At ν=270°, satellite should be at southernmost point (min z)");
    }
    destroy_simulation(sim);
 }
 TEST_CASE("Molniya orbit - orbital period verification", "[inclined][molniya][period][!mayfail]") {
    const double TIME_STEP = 60.0;
    const double SECONDS_PER_HOUR = 3600.0;
    const double MAX_SIMULATION_HOURS = 15.0;
    SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP);
    REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml"));
    Spacecraft* molniya = &sim->spacecraft[0];
    CelestialBody* earth = &sim->bodies[0];
    double semi_major_axis = molniya->orbit.semi_major_axis;
    double mu = G * earth->mass;
    double theoretical_period_seconds = 2.0 * M_PI * sqrt(pow(semi_major_axis, 3) / mu);
    double theoretical_period_hours = theoretical_period_seconds / SECONDS_PER_HOUR;
    INFO("Semi-major axis: " << semi_major_axis << " m");
    INFO("Theoretical period from Kepler's 3rd law: " << theoretical_period_hours << " hours");
    OrbitTracker* tracker = create_orbit_tracker_with_min_time(0, 0.01);
    double max_time = MAX_SIMULATION_HOURS * SECONDS_PER_HOUR;
    while (sim->time < max_time && !tracker->orbit_completed) {
        update_simulation(sim);
        update_orbit_tracker(tracker, (CelestialBody*)molniya, earth, sim->time);
    }
    REQUIRE(tracker->orbit_completed);
    double measured_period_hours = tracker->time_at_completion / SECONDS_PER_HOUR;
    double period_error_hours = fabs(measured_period_hours - theoretical_period_hours);
    INFO("Measured period: " << measured_period_hours << " hours");
    INFO("Period error: " << period_error_hours << " hours");
    INFO("Period error: " << (period_error_hours / theoretical_period_hours * 100.0) << "%");
    REQUIRE(period_error_hours * SECONDS_PER_HOUR < PERIOD_TOLERANCE_SECONDS);
    destroy_orbit_tracker(tracker);
    destroy_simulation(sim);
 }
 TEST_CASE("Generic inclined orbit - moderate inclination", "[inclined][generic][!mayfail]") {
    const double TIME_STEP = 60.0;
    const double SEMI_MAJOR_AXIS = 10000000.0;
    const double ECCENTRICITY = 0.5;
    const double INCLINATION_DEG = 45.0;
    const double INCLINATION_RAD = INCLINATION_DEG * M_PI / 180.0;
    SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP);
    REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml"));
    Spacecraft* craft = &sim->spacecraft[0];
    CelestialBody* earth = &sim->bodies[0];
    craft->orbit.semi_major_axis = SEMI_MAJOR_AXIS;
    craft->orbit.eccentricity = ECCENTRICITY;
    craft->orbit.true_anomaly = 0.0;
    craft->orbit.inclination = INCLINATION_RAD;
    craft->orbit.longitude_of_ascending_node = 0.0;
    craft->orbit.argument_of_periapsis = 0.0;
    initialize_orbital_objects(sim);
    SECTION("Z-coordinate is non-zero for inclined orbit") {
        double z_position = craft->global_position.z;
        INFO("Z-coordinate: " << z_position << " m");
        REQUIRE(z_position != 0.0);
    }
    SECTION("Position magnitude matches orbital radius") {
        double position_vector_mag = vec3_magnitude(craft->global_position);
        double orbital_radius = vec3_magnitude(vec3_sub(craft->global_position, earth->global_position));
        double magnitude_error = fabs(position_vector_mag - orbital_radius);
        INFO("Position vector magnitude: " << position_vector_mag << " m");
        INFO("Orbital radius: " << orbital_radius << " m");
        INFO("Error: " << magnitude_error << " m");
        REQUIRE(magnitude_error < POSITION_TOLERANCE_METERS);
    }
    destroy_simulation(sim);
 }
 TEST_CASE("Inclined orbit - inclination parameter is preserved", "[inclined][config]") {
    const double TIME_STEP = 60.0;
    const double EXPECTED_INCLINATION_RAD = 1.107;
    const double EXPECTED_INCLINATION_DEG = EXPECTED_INCLINATION_RAD * 180.0 / M_PI;
    SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP);
    REQUIRE(load_system_config(sim, "tests/test_inclined_orbits.toml"));
    Spacecraft* molniya = &sim->spacecraft[0];
    INFO("Loaded inclination: " << (molniya->orbit.inclination * 180.0 / M_PI) << " degrees");
    INFO("Expected inclination: " << EXPECTED_INCLINATION_DEG << " degrees");
    REQUIRE_THAT(molniya->orbit.inclination, Catch::Matchers::WithinAbs(EXPECTED_INCLINATION_RAD, 0.01));
    destroy_simulation(sim);
 }
--- a/tests/test_inclined_orbits.toml
+++ b/tests/test_inclined_orbits.toml
@ -0,0 +1,35 @@
 # Test Configuration: Molniya Orbit
 # Earth as root body with highly elliptical, highly inclined satellite orbit
 # Molniya orbit parameters:
 #   - Period: ~718 minutes (~12 hours)
 #   - Eccentricity: 0.74
 #   - Inclination: 63.4°
 #   - Argument of perigee: 270° (apogee at northernmost point)
 #   - Perigee altitude: ~600 km
 #   - Apogee altitude: ~39,700 km
 #   - Semi-major axis: ~26,600 km
 [[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 = "Molniya_Satellite"
 mass = 1000.0
 parent_index = 0
 orbit = {
    semi_major_axis = 26540000.0,
    eccentricity = 0.74,
    true_anomaly = 0.0,
    inclination = 1.107,
    longitude_of_ascending_node = 0.0,
    argument_of_periapsis = 4.71
 }