From cfb2c92bd4ac2e970ae5cd8726c99c17fe7fc832 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Wed, 28 Jan 2026 11:44:15 -0500 Subject: [PATCH] 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 --- docs/planning/molniya-orbit-test-plan.md | 248 +++++++++++++++++++++++ tests/test_inclined_orbits.cpp | 203 +++++++++++++++++++ tests/test_inclined_orbits.toml | 35 ++++ 3 files changed, 486 insertions(+) create mode 100644 docs/planning/molniya-orbit-test-plan.md create mode 100644 tests/test_inclined_orbits.cpp create mode 100644 tests/test_inclined_orbits.toml diff --git a/docs/planning/molniya-orbit-test-plan.md b/docs/planning/molniya-orbit-test-plan.md new file mode 100644 index 0000000..44c526f --- /dev/null +++ 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 diff --git a/tests/test_inclined_orbits.cpp b/tests/test_inclined_orbits.cpp new file mode 100644 index 0000000..72f64a2 --- /dev/null +++ b/tests/test_inclined_orbits.cpp @@ -0,0 +1,203 @@ +#include +#include +#include "../src/physics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include + +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); +} diff --git a/tests/test_inclined_orbits.toml b/tests/test_inclined_orbits.toml new file mode 100644 index 0000000..eba1ec1 --- /dev/null +++ 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 +}