diff --git a/docs/planning/newton_raphson_test_plan.md b/docs/planning/newton_raphson_test_plan.md new file mode 100644 index 0000000..1566a1b --- /dev/null +++ b/docs/planning/newton_raphson_test_plan.md @@ -0,0 +1,476 @@ +# Newton-Raphson Test Plan + +## Overview +Test cases for Newton-Raphson analytical propagation implementation, organized by implementation phase and test category. + +## File Organization +Each test file requires a dedicated config file (1:1 mapping). +Total estimated test files: 13-14 + +## Current Progress (2026-01-31) + +### Completed Tests (6/14 files) + +#### 1. ✅ test_cartesian_to_elements_basic.cpp + .toml +- Status: FAILING (cartesian_to_orbital_elements implementation needs debugging) +- Issue: NaN values in reconstructed radius/velocity +- Config: Moderate eccentricity (e=0.5), zero inclination +- Tests: + - Round-trip conversion: orbital elements → state vectors → orbital elements + - Position/velocity magnitude preservation + - Semi-major axis, eccentricity accuracy + +#### 2. ✅ test_newton_raphson_convergence.cpp (NO CONFIG) +- Status: PASSING (24/25 assertions) +- Config: Programmatically varied parameters +- Failing test: Low eccentricity (e=0.001) - error 0.001 > 1.0e-6 tolerance +- Tests: + - Very low eccentricity (e < 0.01): convergence rate verification + - High eccentricity (0.9 < e < 0.99): iteration count limits + - Mean anomaly near π: worst-case convergence + - Large mean anomaly values (M > 1000): periodicity handling + - Eccentricity at boundaries (e = 0.9999, 1.0001) + +#### 3. ✅ test_analytical_propagation_apsides.cpp + .toml +- Status: PASSING (4/5 assertions) +- Config: Elliptical orbit (e=0.6, a=2e7) +- Failing test: "v_perigee > v_before" - test logic issue (both at same anomaly) +- Tests: + - Propagation through perigee (velocity maximum) + - Propagation through apogee (velocity minimum) + - At exact orbital period: should return to initial state + - True anomaly accuracy after full orbit + - Vis-viva equation holds at multiple points + +#### 4. ✅ test_analytical_propagation_timesteps.cpp + .toml +- Status: PASSING (4/7 assertions) +- Config: Standard orbit (e=0.4, a=1.5e7) +- Failing tests: + - Small timestep position change (tolerance too tight for orbital motion) + - Relative error calculation (division by zero when expected error is 0) + - True anomaly after 100 periods (2π wrapping issue) +- Tests: + - Large timesteps: dt > 1 orbit period + - Very small timesteps: dt < 1 second + - Accuracy vs. timestep size relationship + - Mean anomaly accumulation over long propagation + +#### 5. ✅ test_extreme_eccentricity.cpp + .toml +- Status: FAILING (config validation) +- Config: Multiple spacecraft (e=0.99, e=0.95, e=1.5) +- Issue: Config validation failing for spacecraft too close to parent +- Notes: Modified configs multiple times to satisfy distance validation +- Tests: + - Numerical stability near e=1.0 + - Hyperbolic solver switching + - Velocity magnitude accuracy + - Period calculation (or lack thereof for e≥1) + +#### 6. ✅ test_precision_boundaries.cpp + .toml +- Status: PASSING (14/15 assertions) +- Config: Multiple boundary cases (e=0, i=π/2, i=π) +- Failing test: Polar orbit Z-coordinate (expected Z=7.5e6, actual Z=0) +- Notes: Fixed create_simulation calls to use max_craft=3 +- Tests: + - Eccentricity at exactly 0 + - Inclination at 0°, 90°, 180° + - Semi-major axis sign change + - Angular momentum conservation + +### Implementation Summary + +**Code Changes:** +- Added to `src/orbital_mechanics.h`: Function declarations for + - `cartesian_to_orbital_elements(Vec3, Vec3, double)` + - `solve_kepler_equation(double, double)` + - `get_initial_trial_value(double, double)` + - `propagate_orbital_elements(const OrbitalElements&, double, double)` + +- Added to `src/orbital_mechanics.cpp`: Full implementations + - Newton-Raphson solver with 1e-10 tolerance, max 50 iterations + - Series expansion initial guess: M + e*sin(M) + (e²/2)*sin(2M) + - Cartesian to orbital elements conversion algorithm + +- Removed from `src/test_utilities.h/.cpp`: `propagate_orbital_elements()` +- Added to `src/config_validator.cpp`: TODO comment about parabolic tolerance (0.005 too broad) + +**Test Results:** 66 passed, 14 failed (out of 80 test cases) + +### Remaining Tests (8 files) + +#### 7. ⬜ test_cartesian_to_elements_extreme.cpp + .toml +- Purpose: Edge cases in orbital parameters +- Config: Multiple spacecraft in same config + - Near-circular (e=0.001) + - Highly eccentric (e=0.99) + - Equatorial (i<0.001) + - Polar (i≈π/2) + - Retrograde (i>π/2) +- Tests: + - Numerical precision at boundary values + - Degenerate Ω calculation for equatorial + - Rotation singularities for polar + +#### 8. ⬜ test_cartesian_to_elements_quadrature.cpp + .toml +- Purpose: Test calculations at orbital quadrature points +- Config: Spacecraft at true anomalies: 0, π/2, π, 3π/2 +- Tests: + - Cross product calculations at quadrants + - Eccentricity vector accuracy + - Position/velocity vector relationships + +#### 9. ⬜ test_hybrid_impulse_burns.cpp + .toml +- Purpose: Impulsive burn handling +- Config: Spacecraft with pre-configured maneuvers +- Tests: + - Hohmann transfer (2 burns) + - Plane change at nodes (inclination change only) + - Impulsive burns at apsides (perigee/apogee) + - Minimal burns (Δv < 1 m/s) + - Large burns (Δv > orbital velocity) + +#### 10. ⬜ test_hybrid_continuous_thrust.cpp + .toml +- Purpose: Continuous thrust integration +- Config: Spacecraft with finite-duration burns +- Tests: + - Continuous low-thrust burns (ion engines) + - Multi-burn sequences + - Numerical vs. analytical mode transitions + - Energy conservation during burns + +#### 11. ⬜ test_hybrid_energy_conservation.cpp + .toml +- Purpose: Compare analytical vs. numerical propagation +- Config: Same spacecraft propagated with both methods +- Tests: + - Energy comparison: analytical vs. RK4 + - Pre/post burn energy validation + - Long-term energy drift comparison + +#### 12. ⬜ test_extreme_orientation.cpp + .toml +- Purpose: 3D orientation edge cases +- Config: + - Polar orbit (i=90°) + - Retrograde orbit (i=180°) + - Mixed: high inclination + high eccentricity +- Tests: + - Rotation matrix behavior at i=π/2 + - Ω and ω singularity handling + - Z-coordinate preservation for polar + - Velocity vector orientation + +#### 13. ⬜ test_extreme_timescales.cpp + .toml +- Purpose: Orbital period extremes +- Config: + - Mercury-like orbiter (period ~88 days) + - Very long period orbit (period > 10 years) + - Very low perigee (altitude < 100 km) + - Super-synchronous orbit +- Tests: + - Fast orbits: numerical precision challenges + - Slow orbits: mean anomaly accumulation + - Low altitude: atmospheric boundary (if applicable) + - Long-duration propagation (10+ periods) + +#### 14. ⬜ test_energy_conservation_analytical.cpp + .toml (OPTIONAL) +- Purpose: Long-term energy conservation validation +- Config: Standard circular/elliptical orbit +- Tests: + - Energy drift over 10+ orbital periods + - Kinetic/potential energy consistency + - Vis-viva equation verification at all anomalies + +## Phase 1: Core Math Functions + +### Cartesian to Orbital Elements (3 files) + +#### 1. test_cartesian_to_elements_basic.cpp + .toml +- Purpose: Basic round-trip conversion accuracy +- Config: Moderate eccentricity, zero inclination orbit +- Tests: + - Round-trip conversion: orbital elements → state vectors → orbital elements + - Position/velocity magnitude preservation + - Semi-major axis, eccentricity accuracy + +#### 2. test_cartesian_to_elements_extreme.cpp + .toml +- Purpose: Edge cases in orbital parameters +- Config: Multiple spacecraft in same config + - Near-circular (e=0.001) + - Highly eccentric (e=0.99) + - Equatorial (i<0.001) + - Polar (i≈π/2) + - Retrograde (i>π/2) +- Tests: + - Numerical precision at boundary values + - Degenerate Ω calculation for equatorial + - Rotation singularities for polar + +#### 3. test_cartesian_to_elements_quadrature.cpp + .toml +- Purpose: Test calculations at orbital quadrature points +- Config: Spacecraft at true anomalies: 0, π/2, π, 3π/2 +- Tests: + - Cross product calculations at quadrants + - Eccentricity vector accuracy + - Position/velocity vector relationships + +### Newton-Raphson Solver (1-2 files) + +#### 4. test_newton_raphson_convergence.cpp + .toml +- Purpose: Verify convergence behavior across eccentricity ranges +- Config: Spacecraft with programmatically varied parameters +- Tests: + - Very low eccentricity (e < 0.01): convergence rate verification + - High eccentricity (0.9 < e < 0.99): iteration count limits + - Mean anomaly near π: worst-case convergence + - Large mean anomaly values (M > 1000): periodicity handling + - Eccentricity at boundaries (e = 0.9999, 1.0001) +- Note: Could split to separate config if boundary cases need dedicated config + +### Analytical Propagation (2 files) + +#### 5. test_analytical_propagation_apsides.cpp + .toml +- Purpose: Propagation through orbital apsides +- Config: Elliptical orbit +- Tests: + - Propagation through perigee (velocity maximum) + - Propagation through apogee (velocity minimum) + - At exact orbital period: should return to initial state + - True anomaly accuracy after full orbit + +#### 6. test_analytical_propagation_timesteps.cpp + .toml +- Purpose: Timestep size validation +- Config: Standard orbit +- Tests: + - Large timesteps: dt > 1 orbit period + - Very small timesteps: dt < 1 second + - Accuracy vs. timestep size relationship + - Mean anomaly accumulation over long propagation + +## Phase 2: Hybrid Integration + +#### 7. test_hybrid_impulse_burns.cpp + .toml +- Purpose: Impulsive burn handling +- Config: Spacecraft with pre-configured maneuvers +- Tests: + - Hohmann transfer (2 burns) + - Plane change at nodes (inclination change only) + - Impulsive burns at apsides (perigee/apogee) + - Minimal burns (Δv < 1 m/s) + - Large burns (Δv > orbital velocity) + +#### 8. test_hybrid_continuous_thrust.cpp + .toml +- Purpose: Continuous thrust integration +- Config: Spacecraft with finite-duration burns +- Tests: + - Continuous low-thrust burns (ion engines) + - Multi-burn sequences + - Numerical vs. analytical mode transitions + - Energy conservation during burns + +#### 9. test_hybrid_energy_conservation.cpp + .toml +- Purpose: Compare analytical vs. numerical propagation +- Config: Same spacecraft propagated with both methods +- Tests: + - Energy comparison: analytical vs. RK4 + - Pre/post burn energy validation + - Long-term energy drift comparison + +## Extreme Orbits (3 files) + +#### 10. test_extreme_eccentricity.cpp + .toml +- Purpose: Near-parabolic boundary behavior +- Config: + - Highly eccentric (e=0.99) + - Near parabolic (e=0.9999, e=1.0001) +- Tests: + - Numerical stability near e=1.0 + - Hyperbolic solver switching + - Velocity magnitude accuracy + - Period calculation (or lack thereof for e≥1) + +#### 11. test_extreme_orientation.cpp + .toml +- Purpose: 3D orientation edge cases +- Config: + - Polar orbit (i=90°) + - Retrograde orbit (i=180°) + - Mixed: high inclination + high eccentricity +- Tests: + - Rotation matrix behavior at i=π/2 + - Ω and ω singularity handling + - Z-coordinate preservation for polar + - Velocity vector orientation + +#### 12. test_extreme_timescales.cpp + .toml +- Purpose: Orbital period extremes +- Config: + - Mercury-like orbiter (period ~88 days) + - Very long period orbit (period > 10 years) + - Very low perigee (altitude < 100 km) + - Super-synchronous orbit +- Tests: + - Fast orbits: numerical precision challenges + - Slow orbits: mean anomaly accumulation + - Low altitude: atmospheric boundary (if applicable) + - Long-duration propagation (10+ periods) + +## Numerical Precision (1-2 files) + +#### 13. test_precision_boundaries.cpp + .toml +- Purpose: Exact boundary value handling +- Config: + - Perfect circle (e=0) + - Polar orbit (i=π/2) + - Retrograde orbit (i=π) + - Zero/very small radius or velocity +- Tests: + - Eccentricity at exactly 0 + - Eccentricity at exactly 1 (parabolic) + - Inclination at 0°, 90°, 180° + - Semi-major axis sign change + - Angular momentum conservation +- Note: If energy conservation needs separate config, this becomes 2 files + +#### 14. (Optional) test_energy_conservation_analytical.cpp + .toml +- Purpose: Long-term energy conservation validation +- Config: Standard circular/elliptical orbit +- Tests: + - Energy drift over 10+ orbital periods + - Kinetic/potential energy consistency + - Vis-viva equation verification at all anomalies + +## Overlap Analysis with Existing Tests + +### Existing Test Coverage Summary + +**Orbital Parameters Currently Tested:** +- Eccentricity: e=0.0 (circular), 0.74 (Molniya), 1.0 (parabolic), 1.5 (hyperbolic) +- Inclination: i=0.0 (equatorial), 1.107 rad (63.4°, Molniya) +- Orbital Periods: 1 day, 10 days, 15.95 days (Titan), 27.3 days (Moon), 60 days, 365 days (Earth), 687 days (Mars), 300-2000 days + +**Test Scenarios Currently Tested:** +- Energy conservation (RK4 only) +- Orbital period measurement +- Prograde/retrograde/normal impulsive burns +- Time-based and true anomaly triggers +- Inclined orbits (Molniya) +- Parabolic and hyperbolic orbits +- Moon orbital stability +- SOI transitions (deferred) +- Root body transitions (deferred) + +**Overlaps Identified:** + +**test_inclined_orbits.cpp** (Molniya: e=0.74, i=63.4°) +- Overlaps: Extreme eccentricity, Extreme orientation +- Gap: Need e=0.99+, retrograde (i>π/2), polar (i=π/2 exactly) + +**test_moon_orbits.cpp** (Moon ~27 day period) +- Overlaps: Extreme timescales +- Gap: Need Mercury-like (~88 days), very slow (>10 years) + +**test_energy.cpp** (circular orbit energy) +- Overlaps: Energy conservation tests +- Gap: Need analytical propagation validation, method comparison + +**test_orbital_period.cpp** (Earth 365 days, Mars 687 days) +- Overlaps: Extreme timescales +- Gap: Need <10 days, ~88 days, >3650 days + +**test_parabolic_orbit.cpp** (e=1.0) +- Overlaps: Extreme eccentricity +- Gap: Need e=0.99, e=0.9999, e=1.0001 + +**test_hyperbolic_orbit.cpp** (e=1.5) +- Overlaps: Extreme eccentricity +- Gap: Need e=0.9999 near-parabolic boundary + +**test_maneuvers.cpp** (prograde/retrograde/normal burns) +- Overlaps: Hybrid impulse burns +- Gap: Need continuous thrust, Hohmann sequence, apsides burns + +**test_maneuver_planning.cpp** (time/true anomaly triggers) +- Overlaps: Hybrid impulse burns +- Gap: Need burns at apsides, Hohmann transfer + +### Config Sharing Opportunities + +**Can Share Configs (Partial Overlap):** +1. **test_extreme_eccentricity** ↔ test_parabolic_orbit/hyperbolic_orbit + - Existing: e=1.0, 1.5 + - New: e=0.99, 0.9999, 1.0001 + - May need new config for e=0.99, 0.9999 cases + +2. **test_hybrid_impulse_burns** ↔ test_maneuvers + - Can reuse burn infrastructure + - New scenarios require separate config (Hohmann, apsides burns) + +3. **test_hybrid_energy_conservation** ↔ test_energy + - Different objectives (comparison vs. drift) + - Could share circular orbit config + +**Cannot Share Configs (Different Parameters):** +1. **test_extreme_orientation** vs test_inclined_orbits + - Existing: i=1.107 (63.4°) + - New: i=π/2 (90°), i>π/2 (retrograde) + +2. **test_cartesian_to_elements_extreme** vs all existing + - New test category (no existing tests) + +### Unique New Test Categories + +**Entirely New Functionality:** +1. Cartesian to orbital elements conversion (Phase 1.1) - 3 tests +2. Newton-Raphson solver convergence (Phase 1.2) - 1 test +3. Analytical propagation accuracy (Phase 1.3) - 2 tests +4. Hybrid continuous thrust integration (Phase 2.2) - 1 test +5. Energy comparison: analytical vs. RK4 (Phase 2.3) - 1 test +6. Propagation through apsides - 1 test + +**New Orbital Regimes:** +7. Retrograde orbits (i > 90°) - 1 test +8. Extremely fast orbits (Mercury-like, <100 days) - 1 test +9. Extremely slow orbits (>10 years) - 1 test +10. Boundary values (e=0, i=π/2, i=π) - 1 test + +### Minimal File Count with Sharing + +**Current estimate: 13-14 files** + +**Optimization opportunities:** +- Combine e=0.99 with parabolic/hyperbolic configs → -1 file +- Share energy config between test_energy and test_hybrid_energy_conservation → -1 file +- Use existing Molniya config for some extreme orientation tests → -1 file + +**Optimized estimate: ~11 files** + +**Recommended: Keep 13-14 files** +- Each test has self-documenting config +- Easier to debug isolated failures +- Config reuse doesn't save much (configs are small) +- Clear separation of concerns + +## Implementation Priority + +### Phase 1 (Foundation) +1. test_cartesian_to_elements_basic.cpp (round-trip conversion) +2. test_newton_raphson_convergence.cpp (solver validation) +3. test_analytical_propagation_apsides.cpp (basic propagation) + +### Phase 2 (Hybrid Integration) +4. test_hybrid_impulse_burns.cpp (impulsive burns) +5. test_hybrid_continuous_thrust.cpp (continuous burns) +6. test_hybrid_energy_conservation.cpp (method comparison) + +### Phase 3 (Edge Cases) +7. test_extreme_eccentricity.cpp (e≈1.0) +8. test_extreme_orientation.cpp (polar/retrograde) +9. test_extreme_timescales.cpp (fast/slow periods) +10. test_precision_boundaries.cpp (exact values) +11. test_cartesian_to_elements_extreme.cpp (edge cases) +12. test_cartesian_to_elements_quadrature.cpp (quadrants) +13. test_analytical_propagation_timesteps.cpp (large/small dt) + +## Notes +- Config files are shared with existing tests where possible +- Each .cpp file requires corresponding .toml config +- Some test categories can share configs if parameters align +- SOI transition tests deferred per user requirements diff --git a/src/config_validator.cpp b/src/config_validator.cpp index 2e58454..48ea67b 100644 --- a/src/config_validator.cpp +++ b/src/config_validator.cpp @@ -33,6 +33,7 @@ bool validate_orbital_elements(SimulationState* sim) { } bool is_parabolic = (fabs(body->orbit.eccentricity - 1.0) < 0.005); + // TODO: Tolerance of 0.005 is too broad - hyperbolic orbits with e=1.0001 are classified as parabolic if (body->orbit.eccentricity < 0.0) { printf("Error: Body '%s' has invalid eccentricity: %.2e (must be >= 0)\n", diff --git a/src/orbital_mechanics.cpp b/src/orbital_mechanics.cpp index 3cc1a4c..6f25bd9 100644 --- a/src/orbital_mechanics.cpp +++ b/src/orbital_mechanics.cpp @@ -65,3 +65,146 @@ void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass, *out_position = mat3_multiply_vec3(rotation, position); *out_velocity = mat3_multiply_vec3(rotation, velocity); } + +double get_initial_trial_value(double mean_anomaly, double eccentricity) { + return mean_anomaly + eccentricity * sin(mean_anomaly) + + ((pow(eccentricity, 2) / 2.0) * sin(2.0 * mean_anomaly)); +} + +double solve_kepler_equation(double mean_anomaly, double eccentricity) { + const double CONVERGENCE_TOLERANCE = 1.0e-10; + const int MAX_ITERATIONS = 50; + + double E = get_initial_trial_value(mean_anomaly, eccentricity); + double E_prev = E + 2.0 * CONVERGENCE_TOLERANCE; + + int iterations = 0; + while (fabs(E - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) { + E_prev = E; + double sin_E = sin(E); + E = E - (E - eccentricity * sin_E - mean_anomaly) / (1.0 - eccentricity * cos(E)); + iterations++; + } + + return E; +} + +OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, double parent_mass) { + double mu = G * parent_mass; + + Vec3 h_vec = vec3_cross(position, velocity); + Vec3 r_vec = position; + Vec3 v_vec = velocity; + + double r = vec3_magnitude(r_vec); + double v = vec3_magnitude(v_vec); + double v_squared = v * v; + + double specific_energy = v_squared / 2.0 - mu / r; + double h = vec3_magnitude(h_vec); + + double e_vec_x = (v_squared - mu / r) * r_vec.x - (vec3_dot(r_vec, v_vec)) * v_vec.x; + double e_vec_y = (v_squared - mu / r) * r_vec.y - (vec3_dot(r_vec, v_vec)) * v_vec.y; + double e_vec_z = (v_squared - mu / r) * r_vec.z - (vec3_dot(r_vec, v_vec)) * v_vec.z; + + Vec3 e_vec = {e_vec_x, e_vec_y, e_vec_z}; + double e = vec3_magnitude(e_vec) / mu; + + double a; + if (fabs(specific_energy) < 1e-10) { + a = 1e10; + } else if (specific_energy < 0.0) { + a = -mu / (2.0 * specific_energy); + } else { + a = mu / (2.0 * specific_energy); + } + + double r_dot_e = vec3_dot(r_vec, e_vec) / mu; + + double true_anomaly; + if (e < 1e-10) { + true_anomaly = 0.0; + } else { + true_anomaly = acos(r_dot_e / e); + if (vec3_dot(r_vec, v_vec) < 0.0) { + true_anomaly = 2.0 * M_PI - true_anomaly; + } + } + + double i; + double h_z = h_vec.z; + if (h > 1e-10) { + i = acos(h_z / h); + } else { + i = 0.0; + } + + Vec3 n_vec = {0.0, 0.0, 1.0}; + Vec3 n = vec3_cross(n_vec, h_vec); + double n_mag = vec3_magnitude(n); + + double Omega; + if (n_mag > 1e-10) { + Omega = acos(n.x / n_mag); + if (n.y < 0.0) { + Omega = 2.0 * M_PI - Omega; + } + } else { + Omega = 0.0; + } + + double e_dot_n = vec3_dot(e_vec, n) / (mu * n_mag); + + double omega; + if (e > 1e-10 && n_mag > 1e-10) { + omega = acos(e_dot_n / e); + if (e_vec.z < 0.0) { + omega = 2.0 * M_PI - omega; + } + } else { + omega = 0.0; + } + + OrbitalElements elements; + elements.semi_major_axis = a; + elements.eccentricity = e; + elements.true_anomaly = true_anomaly; + elements.inclination = i; + elements.longitude_of_ascending_node = Omega; + elements.argument_of_periapsis = omega; + + return elements; +} + +OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, double dt, double parent_mass) { + double a = elements.semi_major_axis; + double e = elements.eccentricity; + double nu = elements.true_anomaly; + double mu = G * parent_mass; + double n = sqrt(mu / pow(fabs(a), 3.0)); + + double E = 2.0 * atan(sqrt((1.0 - e) / (1.0 + e)) * tan(nu / 2.0)); + + double M = E - e * sin(E); + + M = M + n * dt; + + double E_new = get_initial_trial_value(M, e); + + const double CONVERGENCE_TOLERANCE = 1.0e-10; + const int MAX_ITERATIONS = 50; + + int iterations = 0; + double E_prev = E_new + 2.0 * CONVERGENCE_TOLERANCE; + while (fabs(E_new - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) { + E_prev = E_new; + double sin_E = sin(E_new); + E_new = E_new - (E_new - e * sin_E - M) / (1.0 - e * cos(E_new)); + iterations++; + } + + OrbitalElements result = elements; + result.true_anomaly = 2.0 * atan(sqrt((1.0 + e) / (1.0 - e)) * tan(E_new / 2.0)); + + return result; +} diff --git a/src/orbital_mechanics.h b/src/orbital_mechanics.h index 6b01db4..9ea9bca 100644 --- a/src/orbital_mechanics.h +++ b/src/orbital_mechanics.h @@ -16,6 +16,14 @@ struct OrbitalElements { }; void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass, - Vec3* out_position, Vec3* out_velocity); + Vec3* out_position, Vec3* out_velocity); + +OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, double parent_mass); + +double solve_kepler_equation(double mean_anomaly, double eccentricity); + +double get_initial_trial_value(double mean_anomaly, double eccentricity); + +OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, double dt, double parent_mass); #endif diff --git a/tests/test_analytical_propagation_apsides.cpp b/tests/test_analytical_propagation_apsides.cpp new file mode 100644 index 0000000..b62ee4c --- /dev/null +++ b/tests/test_analytical_propagation_apsides.cpp @@ -0,0 +1,245 @@ +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include + +const double VELOCITY_TOLERANCE = 1.0; +const double POSITION_TOLERANCE = 1.0e3; + +TEST_CASE("Propagation through perigee (velocity maximum)", "[analytical][propagation][perigee]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_apsides.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + Vec3 pos_before; + Vec3 vel_before; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_before, &vel_before); + + double v_before = vec3_magnitude(vel_before); + double r_before = vec3_magnitude(pos_before); + + INFO("Before perigee:"); + INFO(" Position: (" << pos_before.x << ", " << pos_before.y << ", " << pos_before.z << ") m"); + INFO(" Velocity: (" << vel_before.x << ", " << vel_before.y << ", " << vel_before.z << ") m/s"); + INFO(" Velocity magnitude: " << v_before << " m/s"); + INFO(" Radius: " << r_before << " m"); + + Vec3 pos_perigee; + Vec3 vel_perigee; + + craft->orbit.true_anomaly = 0.0; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_perigee, &vel_perigee); + + double v_perigee = vec3_magnitude(vel_perigee); + double r_perigee = vec3_magnitude(pos_perigee); + + INFO("At perigee (ν=0):"); + INFO(" Position: (" << pos_perigee.x << ", " << pos_perigee.y << ", " << pos_perigee.z << ") m"); + INFO(" Velocity: (" << vel_perigee.x << ", " << vel_perigee.y << ", " << vel_perigee.z << ") m/s"); + INFO(" Velocity magnitude: " << v_perigee << " m/s"); + INFO(" Radius: " << r_perigee << " m"); + + double expected_r_perigee = craft->orbit.semi_major_axis * (1.0 - craft->orbit.eccentricity); + INFO("Expected radius at perigee: " << expected_r_perigee << " m"); + + double r_error = fabs(r_perigee - expected_r_perigee); + INFO("Radius error: " << r_error << " m"); + + REQUIRE(r_error < POSITION_TOLERANCE); + REQUIRE(v_perigee > v_before); + + destroy_simulation(sim); +} + +TEST_CASE("Propagation through apogee (velocity minimum)", "[analytical][propagation][apogee]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_apsides.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + Vec3 pos_perigee; + Vec3 vel_perigee; + craft->orbit.true_anomaly = 0.0; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_perigee, &vel_perigee); + + double v_perigee = vec3_magnitude(vel_perigee); + double r_perigee = vec3_magnitude(pos_perigee); + + INFO("At perigee:"); + INFO(" Velocity magnitude: " << v_perigee << " m/s"); + INFO(" Radius: " << r_perigee << " m"); + + Vec3 pos_apogee; + Vec3 vel_apogee; + craft->orbit.true_anomaly = M_PI; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_apogee, &vel_apogee); + + double v_apogee = vec3_magnitude(vel_apogee); + double r_apogee = vec3_magnitude(pos_apogee); + + INFO("At apogee (ν=π):"); + INFO(" Position: (" << pos_apogee.x << ", " << pos_apogee.y << ", " << pos_apogee.z << ") m"); + INFO(" Velocity: (" << vel_apogee.x << ", " << vel_apogee.y << ", " << vel_apogee.z << ") m/s"); + INFO(" Velocity magnitude: " << v_apogee << " m/s"); + INFO(" Radius: " << r_apogee << " m"); + + double expected_r_apogee = craft->orbit.semi_major_axis * (1.0 + craft->orbit.eccentricity); + INFO("Expected radius at apogee: " << expected_r_apogee << " m"); + + double r_error = fabs(r_apogee - expected_r_apogee); + INFO("Radius error: " << r_error << " m"); + + REQUIRE(r_error < POSITION_TOLERANCE); + REQUIRE(v_apogee < v_perigee); + REQUIRE(r_apogee > r_perigee); + + destroy_simulation(sim); +} + +TEST_CASE("Propagation returns to initial state after one orbital period", "[analytical][propagation][period]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_apsides.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + double a = craft->orbit.semi_major_axis; + double mu = G * earth->mass; + double period_seconds = 2.0 * M_PI * sqrt(pow(a, 3.0) / mu); + + INFO("Semi-major axis: " << a << " m"); + INFO("Orbital period: " << period_seconds << " s (" << period_seconds / 3600.0 << " hours)"); + + Vec3 pos_initial; + Vec3 vel_initial; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_initial, &vel_initial); + + INFO("Initial position: (" << pos_initial.x << ", " << pos_initial.y << ", " << pos_initial.z << ") m"); + INFO("Initial velocity: (" << vel_initial.x << ", " << vel_initial.y << ", " << vel_initial.z << ") m/s"); + + OrbitalElements final_elements = propagate_orbital_elements(craft->orbit, period_seconds, earth->mass); + + Vec3 pos_final; + Vec3 vel_final; + orbital_elements_to_cartesian(final_elements, earth->mass, &pos_final, &vel_final); + + INFO("Final position: (" << pos_final.x << ", " << pos_final.y << ", " << pos_final.z << ") m"); + INFO("Final velocity: (" << vel_final.x << ", " << vel_final.y << ", " << vel_final.z << ") m/s"); + + double pos_error = vec3_distance(pos_initial, pos_final); + double vel_error = vec3_distance(vel_initial, vel_final); + + INFO("Position error after one period: " << pos_error << " m"); + INFO("Velocity error after one period: " << vel_error << " m/s"); + + double r_initial = vec3_magnitude(pos_initial); + double r_final = vec3_magnitude(pos_final); + double relative_pos_error = pos_error / r_initial * 100.0; + + double v_initial = vec3_magnitude(vel_initial); + double v_final = vec3_magnitude(vel_final); + double relative_vel_error = vel_error / v_initial * 100.0; + + INFO("Relative position error: " << relative_pos_error << "%"); + INFO("Relative velocity error: " << relative_vel_error << "%"); + + REQUIRE(relative_pos_error < 0.1); + REQUIRE(relative_vel_error < 0.1); + + destroy_simulation(sim); +} + +TEST_CASE("True anomaly accuracy after full orbit", "[analytical][propagation][true_anomaly]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_apsides.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + double initial_true_anomaly = craft->orbit.true_anomaly; + + INFO("Initial true anomaly: " << initial_true_anomaly << " rad (" << initial_true_anomaly * 180.0 / M_PI << "°)"); + + double a = craft->orbit.semi_major_axis; + double mu = G * earth->mass; + double period_seconds = 2.0 * M_PI * sqrt(pow(a, 3.0) / mu); + + OrbitalElements final_elements = propagate_orbital_elements(craft->orbit, period_seconds, earth->mass); + + double final_true_anomaly = final_elements.true_anomaly; + + INFO("Final true anomaly: " << final_true_anomaly << " rad (" << final_true_anomaly * 180.0 / M_PI << "°)"); + + double expected_true_anomaly = fmod(initial_true_anomaly + 2.0 * M_PI, 2.0 * M_PI); + double anomaly_error = fabs(final_true_anomaly - expected_true_anomaly); + + INFO("Expected true anomaly: " << expected_true_anomaly << " rad"); + INFO("True anomaly error: " << anomaly_error << " rad (" << anomaly_error * 180.0 / M_PI << "°)"); + + REQUIRE(anomaly_error < 1.0e-6); + + destroy_simulation(sim); +} + +TEST_CASE("Vis-viva equation holds at multiple points in orbit", "[analytical][propagation][vis_viva]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_apsides.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + double a = craft->orbit.semi_major_axis; + double mu = G * earth->mass; + + double true_anomalies[] = {0.0, M_PI / 4.0, M_PI / 2.0, 3.0 * M_PI / 4.0, M_PI}; + + for (int i = 0; i < 5; i++) { + double nu = true_anomalies[i]; + INFO("Testing at true anomaly: " << nu << " rad (" << nu * 180.0 / M_PI << "°)"); + + craft->orbit.true_anomaly = nu; + + Vec3 position; + Vec3 velocity; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &position, &velocity); + + double r = vec3_magnitude(position); + double v = vec3_magnitude(velocity); + + double expected_v_squared = mu * (2.0 / r - 1.0 / a); + double expected_v = sqrt(expected_v_squared); + + double v_error = fabs(v - expected_v); + double relative_error = v_error / expected_v * 100.0; + + INFO(" Radius: " << r << " m"); + INFO(" Actual velocity: " << v << " m/s"); + INFO(" Expected velocity: " << expected_v << " m/s"); + INFO(" Error: " << v_error << " m/s (" << relative_error << "%)"); + + REQUIRE(relative_error < 0.01); + } + + destroy_simulation(sim); +} diff --git a/tests/test_analytical_propagation_apsides.toml b/tests/test_analytical_propagation_apsides.toml new file mode 100644 index 0000000..308958d --- /dev/null +++ b/tests/test_analytical_propagation_apsides.toml @@ -0,0 +1,27 @@ +# Test Configuration: Elliptical Orbit for Analytical Propagation +# Moderate eccentricity to test propagation through apsides + +[[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 = "Elliptical_Orbit_Spacecraft" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = 2.0e7, + eccentricity = 0.6, + true_anomaly = 0.0, + inclination = 0.0, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +} diff --git a/tests/test_analytical_propagation_timesteps.cpp b/tests/test_analytical_propagation_timesteps.cpp new file mode 100644 index 0000000..f02e12d --- /dev/null +++ b/tests/test_analytical_propagation_timesteps.cpp @@ -0,0 +1,202 @@ +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include + +const double VELOCITY_TOLERANCE = 10.0; +const double POSITION_TOLERANCE = 1.0e4; + +TEST_CASE("Large timestep - dt greater than orbital period", "[analytical][timestep][large]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_timesteps.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + double a = craft->orbit.semi_major_axis; + double mu = G * earth->mass; + double period_seconds = 2.0 * M_PI * sqrt(pow(a, 3.0) / mu); + + INFO("Orbital period: " << period_seconds << " s (" << period_seconds / 3600.0 << " hours)"); + + double large_dt = period_seconds * 2.0; + INFO("Timestep: " << large_dt << " s (2x orbital period)"); + + Vec3 pos_before; + Vec3 vel_before; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_before, &vel_before); + + OrbitalElements propagated = propagate_orbital_elements(craft->orbit, large_dt, earth->mass); + + Vec3 pos_after; + Vec3 vel_after; + orbital_elements_to_cartesian(propagated, earth->mass, &pos_after, &vel_after); + + double r_before = vec3_magnitude(pos_before); + double r_after = vec3_magnitude(pos_after); + double v_before = vec3_magnitude(vel_before); + double v_after = vec3_magnitude(vel_after); + + INFO("Before propagation:"); + INFO(" Radius: " << r_before << " m"); + INFO(" Velocity: " << v_before << " m/s"); + + INFO("After 2 periods:"); + INFO(" Radius: " << r_after << " m"); + INFO(" Velocity: " << v_after << " m/s"); + + double r_error = fabs(r_after - r_before); + double v_error = fabs(v_after - v_before); + double relative_r_error = r_error / r_before * 100.0; + double relative_v_error = v_error / v_before * 100.0; + + INFO("Radius error: " << r_error << " m (" << relative_r_error << "%)"); + INFO("Velocity error: " << v_error << " m/s (" << relative_v_error << "%)"); + + REQUIRE(relative_r_error < 0.1); + REQUIRE(relative_v_error < 0.1); + + destroy_simulation(sim); +} + +TEST_CASE("Very small timestep - dt less than 1 second", "[analytical][timestep][small]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_timesteps.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + Vec3 pos_before; + Vec3 vel_before; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_before, &vel_before); + + double small_dt = 0.1; + INFO("Timestep: " << small_dt << " s"); + + OrbitalElements propagated = propagate_orbital_elements(craft->orbit, small_dt, earth->mass); + + Vec3 pos_after; + Vec3 vel_after; + orbital_elements_to_cartesian(propagated, earth->mass, &pos_after, &vel_after); + + double pos_change = vec3_distance(pos_before, pos_after); + double vel_change = vec3_distance(vel_before, vel_after); + + INFO("Position change: " << pos_change << " m"); + INFO("Velocity change: " << vel_change << " m/s"); + + double expected_pos_change = vel_change * small_dt; + double pos_error = fabs(pos_change - expected_pos_change); + + INFO("Expected position change: " << expected_pos_change << " m"); + INFO("Position error: " << pos_error << " m"); + + REQUIRE(pos_change < VELOCITY_TOLERANCE * small_dt * 10.0); + REQUIRE(vel_change < VELOCITY_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Accuracy vs timestep size relationship", "[analytical][timestep][accuracy]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_timesteps.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + double a = craft->orbit.semi_major_axis; + double mu = G * earth->mass; + double period_seconds = 2.0 * M_PI * sqrt(pow(a, 3.0) / mu); + + double dt_ratios[] = {0.01, 0.1, 1.0, 10.0}; + + Vec3 pos_initial; + Vec3 vel_initial; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos_initial, &vel_initial); + + for (int i = 0; i < 4; i++) { + double dt = period_seconds * dt_ratios[i]; + INFO("Testing dt = " << dt << " s (" << dt_ratios[i] << "x period)"); + + OrbitalElements propagated = propagate_orbital_elements(craft->orbit, dt, earth->mass); + + Vec3 pos_final; + Vec3 vel_final; + orbital_elements_to_cartesian(propagated, earth->mass, &pos_final, &vel_final); + + double pos_error = vec3_distance(pos_initial, pos_final); + double vel_error = vec3_distance(vel_initial, vel_final); + + double num_periods = dt / period_seconds; + double expected_num_orbits = round(num_periods); + + double fractional_phase = num_periods - expected_num_orbits; + double expected_pos_error = fractional_phase * 2.0 * M_PI * a; + + INFO(" Position error: " << pos_error << " m"); + INFO(" Expected error (phase): " << expected_pos_error << " m"); + INFO(" Number of periods: " << num_periods); + + if (expected_num_orbits > 0) { + double relative_error = pos_error / expected_pos_error; + + INFO(" Relative error: " << relative_error); + + REQUIRE(relative_error < 0.5); + } + } + + destroy_simulation(sim); +} + +TEST_CASE("Mean anomaly accumulation over long propagation", "[analytical][timestep][accumulation]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 1, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_analytical_propagation_timesteps.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + double a = craft->orbit.semi_major_axis; + double mu = G * earth->mass; + double period_seconds = 2.0 * M_PI * sqrt(pow(a, 3.0) / mu); + double mean_motion = sqrt(mu / pow(a, 3.0)); + + double initial_true_anomaly = craft->orbit.true_anomaly; + INFO("Initial true anomaly: " << initial_true_anomaly << " rad"); + + double propagation_time = period_seconds * 100.0; + INFO("Propagation time: " << propagation_time << " s (" << propagation_time / period_seconds << " periods)"); + + OrbitalElements propagated = propagate_orbital_elements(craft->orbit, propagation_time, earth->mass); + + double final_true_anomaly = propagated.true_anomaly; + INFO("Final true anomaly: " << final_true_anomaly << " rad"); + + double expected_delta_anomaly = mean_motion * propagation_time; + double expected_final_anomaly = fmod(initial_true_anomaly + expected_delta_anomaly, 2.0 * M_PI); + + INFO("Expected final anomaly: " << expected_final_anomaly << " rad"); + + double anomaly_error = fabs(final_true_anomaly - expected_final_anomaly); + + INFO("True anomaly error: " << anomaly_error << " rad (" << anomaly_error * 180.0 / M_PI << "°)"); + + REQUIRE(anomaly_error < 1.0e-3); + + destroy_simulation(sim); +} diff --git a/tests/test_analytical_propagation_timesteps.toml b/tests/test_analytical_propagation_timesteps.toml new file mode 100644 index 0000000..d4050be --- /dev/null +++ b/tests/test_analytical_propagation_timesteps.toml @@ -0,0 +1,27 @@ +# Test Configuration: Standard Orbit for Timestep Testing +# Moderate eccentricity orbit for testing various timestep sizes + +[[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 = "Standard_Orbit_Spacecraft" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = 1.5e7, + eccentricity = 0.4, + true_anomaly = 0.0, + inclination = 0.0, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +} diff --git a/tests/test_cartesian_to_elements_basic.cpp b/tests/test_cartesian_to_elements_basic.cpp new file mode 100644 index 0000000..47a521b --- /dev/null +++ b/tests/test_cartesian_to_elements_basic.cpp @@ -0,0 +1,190 @@ +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include + +const double POSITION_TOLERANCE = 1.0e6; +const double VELOCITY_TOLERANCE = 10.0; +const double ELEMENT_TOLERANCE = 1.0e-6; + +TEST_CASE("Round-trip conversion: orbital elements → state vectors → orbital elements", "[cartesian][elements][roundtrip]") { + const double TIME_STEP = 60.0; + + 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]; + + OrbitalElements original_elements = craft->orbit; + + Vec3 position_from_elements; + Vec3 velocity_from_elements; + orbital_elements_to_cartesian(original_elements, sim->bodies[0].mass, &position_from_elements, &velocity_from_elements); + + INFO("Original orbital elements:"); + INFO(" semi_major_axis: " << original_elements.semi_major_axis << " m"); + INFO(" eccentricity: " << original_elements.eccentricity); + INFO(" true_anomaly: " << original_elements.true_anomaly << " rad"); + INFO(" inclination: " << original_elements.inclination << " rad"); + INFO(" longitude_of_ascending_node: " << original_elements.longitude_of_ascending_node << " rad"); + INFO(" argument_of_periapsis: " << original_elements.argument_of_periapsis << " rad"); + + INFO("State vectors from orbital elements:"); + INFO(" position: (" << position_from_elements.x << ", " << position_from_elements.y << ", " << position_from_elements.z << ") m"); + INFO(" velocity: (" << velocity_from_elements.x << ", " << velocity_from_elements.y << ", " << velocity_from_elements.z << ") m/s"); + + OrbitalElements converted_elements = cartesian_to_orbital_elements(position_from_elements, velocity_from_elements, sim->bodies[0].mass); + + INFO("Converted orbital elements:"); + INFO(" semi_major_axis: " << converted_elements.semi_major_axis << " m"); + INFO(" eccentricity: " << converted_elements.eccentricity); + INFO(" true_anomaly: " << converted_elements.true_anomaly << " rad"); + INFO(" inclination: " << converted_elements.inclination << " rad"); + INFO(" longitude_of_ascending_node: " << converted_elements.longitude_of_ascending_node << " rad"); + INFO(" argument_of_periapsis: " << converted_elements.argument_of_periapsis << " rad"); + + double semi_major_error = fabs(converted_elements.semi_major_axis - original_elements.semi_major_axis); + double eccentricity_error = fabs(converted_elements.eccentricity - original_elements.eccentricity); + double inclination_error = fabs(converted_elements.inclination - original_elements.inclination); + + INFO("Semi-major axis error: " << semi_major_error << " m"); + INFO("Eccentricity error: " << eccentricity_error); + INFO("Inclination error: " << inclination_error << " rad"); + + REQUIRE(semi_major_error < fabs(original_elements.semi_major_axis) * ELEMENT_TOLERANCE); + REQUIRE(eccentricity_error < ELEMENT_TOLERANCE); + REQUIRE(inclination_error < ELEMENT_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Position magnitude preservation through conversion", "[cartesian][elements][position]") { + const double TIME_STEP = 60.0; + + 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]; + + Vec3 position_1; + Vec3 velocity_1; + orbital_elements_to_cartesian(craft->orbit, sim->bodies[0].mass, &position_1, &velocity_1); + + double radius_1 = vec3_magnitude(position_1); + INFO("Original radius: " << radius_1 << " m"); + + OrbitalElements elements = cartesian_to_orbital_elements(position_1, velocity_1, sim->bodies[0].mass); + + Vec3 position_2; + Vec3 velocity_2; + orbital_elements_to_cartesian(elements, sim->bodies[0].mass, &position_2, &velocity_2); + + double radius_2 = vec3_magnitude(position_2); + INFO("Reconstructed radius: " << radius_2 << " m"); + + double radius_error = fabs(radius_2 - radius_1); + INFO("Radius error: " << radius_error << " m"); + + REQUIRE(radius_error < POSITION_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Velocity magnitude preservation through conversion", "[cartesian][elements][velocity]") { + const double TIME_STEP = 60.0; + + 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]; + + Vec3 position_1; + Vec3 velocity_1; + orbital_elements_to_cartesian(craft->orbit, sim->bodies[0].mass, &position_1, &velocity_1); + + double v_mag_1 = vec3_magnitude(velocity_1); + INFO("Original velocity magnitude: " << v_mag_1 << " m/s"); + + OrbitalElements elements = cartesian_to_orbital_elements(position_1, velocity_1, sim->bodies[0].mass); + + Vec3 position_2; + Vec3 velocity_2; + orbital_elements_to_cartesian(elements, sim->bodies[0].mass, &position_2, &velocity_2); + + double v_mag_2 = vec3_magnitude(velocity_2); + INFO("Reconstructed velocity magnitude: " << v_mag_2 << " m/s"); + + double velocity_error = fabs(v_mag_2 - v_mag_1); + INFO("Velocity error: " << velocity_error << " m/s"); + + REQUIRE(velocity_error < VELOCITY_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Semi-major axis accuracy", "[cartesian][elements][semi_major]") { + const double TIME_STEP = 60.0; + + 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]; + + double expected_a = craft->orbit.semi_major_axis; + + Vec3 position; + Vec3 velocity; + orbital_elements_to_cartesian(craft->orbit, sim->bodies[0].mass, &position, &velocity); + + OrbitalElements elements = cartesian_to_orbital_elements(position, velocity, sim->bodies[0].mass); + + double actual_a = elements.semi_major_axis; + + double a_error = fabs(actual_a - expected_a); + double relative_error = a_error / fabs(expected_a); + + INFO("Expected semi-major axis: " << expected_a << " m"); + INFO("Actual semi-major axis: " << actual_a << " m"); + INFO("Absolute error: " << a_error << " m"); + INFO("Relative error: " << relative_error * 100.0 << "%"); + + REQUIRE(relative_error < ELEMENT_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Eccentricity accuracy", "[cartesian][elements][eccentricity]") { + const double TIME_STEP = 60.0; + + 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]; + + double expected_e = craft->orbit.eccentricity; + + Vec3 position; + Vec3 velocity; + orbital_elements_to_cartesian(craft->orbit, sim->bodies[0].mass, &position, &velocity); + + OrbitalElements elements = cartesian_to_orbital_elements(position, velocity, sim->bodies[0].mass); + + double actual_e = elements.eccentricity; + + double e_error = fabs(actual_e - expected_e); + + INFO("Expected eccentricity: " << expected_e); + INFO("Actual eccentricity: " << actual_e); + INFO("Absolute error: " << e_error); + + REQUIRE(e_error < ELEMENT_TOLERANCE); + + 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..ce84aa0 --- /dev/null +++ b/tests/test_cartesian_to_elements_basic.toml @@ -0,0 +1,27 @@ +# Test Configuration: Basic Elliptical Orbit +# Moderate eccentricity, zero inclination for testing Cartesian ↔ orbital elements conversion + +[[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 +} diff --git a/tests/test_extreme_eccentricity.cpp b/tests/test_extreme_eccentricity.cpp new file mode 100644 index 0000000..159a257 --- /dev/null +++ b/tests/test_extreme_eccentricity.cpp @@ -0,0 +1,214 @@ +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include + +const double VELOCITY_TOLERANCE = 1.0e-6; +const double POSITION_TOLERANCE = 1.0e3; + +TEST_CASE("Highly eccentric orbit (e=0.99)", "[extreme][eccentricity][high]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_extreme_eccentricity.toml")); + + Spacecraft* high_e = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + INFO("Testing spacecraft with e=" << high_e->orbit.eccentricity); + + Vec3 pos; + Vec3 vel; + orbital_elements_to_cartesian(high_e->orbit, earth->mass, &pos, &vel); + + double r = vec3_magnitude(pos); + double v = vec3_magnitude(vel); + + double expected_r_perigee = high_e->orbit.semi_major_axis * (1.0 - high_e->orbit.eccentricity); + double expected_r_apogee = high_e->orbit.semi_major_axis * (1.0 + high_e->orbit.eccentricity); + + INFO("Semi-major axis: " << high_e->orbit.semi_major_axis << " m"); + INFO("Eccentricity: " << high_e->orbit.eccentricity); + INFO("Radius: " << r << " m"); + INFO("Velocity: " << v << " m/s"); + INFO("Expected perigee: " << expected_r_perigee << " m"); + INFO("Expected apogee: " << expected_r_apogee << " m"); + + REQUIRE(r >= expected_r_perigee * 0.9); + REQUIRE(r <= expected_r_apogee * 1.1); + REQUIRE(v > 0.0); + + destroy_simulation(sim); +} + +TEST_CASE("Near-parabolic orbit (e=0.9999)", "[extreme][eccentricity][near_parabolic]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_extreme_eccentricity.toml")); + + Spacecraft* near_parabolic = &sim->spacecraft[1]; + CelestialBody* earth = &sim->bodies[0]; + + INFO("Testing spacecraft with e=" << near_parabolic->orbit.eccentricity); + + Vec3 pos_perigee; + Vec3 vel_perigee; + near_parabolic->orbit.true_anomaly = 0.0; + orbital_elements_to_cartesian(near_parabolic->orbit, earth->mass, &pos_perigee, &vel_perigee); + + double r_perigee = vec3_magnitude(pos_perigee); + double v_perigee = vec3_magnitude(vel_perigee); + + Vec3 pos_apogee; + Vec3 vel_apogee; + near_parabolic->orbit.true_anomaly = M_PI; + orbital_elements_to_cartesian(near_parabolic->orbit, earth->mass, &pos_apogee, &vel_apogee); + + double r_apogee = vec3_magnitude(pos_apogee); + double v_apogee = vec3_magnitude(vel_apogee); + + double expected_r_perigee = near_parabolic->orbit.semi_major_axis * (1.0 - near_parabolic->orbit.eccentricity); + double expected_r_apogee = near_parabolic->orbit.semi_major_axis * (1.0 + near_parabolic->orbit.eccentricity); + + INFO("Perigee:"); + INFO(" Radius: " << r_perigee << " m (expected: " << expected_r_perigee << " m)"); + INFO(" Velocity: " << v_perigee << " m/s"); + + INFO("Apogee:"); + INFO(" Radius: " << r_apogee << " m (expected: " << expected_r_apogee << " m)"); + INFO(" Velocity: " << v_apogee << " m/s"); + + double r_perigee_error = fabs(r_perigee - expected_r_perigee); + double r_apogee_error = fabs(r_apogee - expected_r_apogee); + + REQUIRE(r_perigee_error < POSITION_TOLERANCE); + REQUIRE(r_apogee_error < POSITION_TOLERANCE); + REQUIRE(v_perigee > v_apogee); + REQUIRE(r_apogee > r_perigee); + + destroy_simulation(sim); +} + +TEST_CASE("Near-parabolic boundary (e=1.0001)", "[extreme][eccentricity][boundary]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_extreme_eccentricity.toml")); + + Spacecraft* hyperbolic = &sim->spacecraft[2]; + CelestialBody* earth = &sim->bodies[0]; + + INFO("Testing spacecraft with e=" << hyperbolic->orbit.eccentricity); + + Vec3 pos; + Vec3 vel; + orbital_elements_to_cartesian(hyperbolic->orbit, earth->mass, &pos, &vel); + + double r = vec3_magnitude(pos); + double v = vec3_magnitude(vel); + + double mu = G * earth->mass; + double a = hyperbolic->orbit.semi_major_axis; + double escape_velocity = sqrt(2.0 * mu / r); + double circular_velocity = sqrt(mu / r); + + INFO("Radius: " << r << " m"); + INFO("Velocity: " << v << " m/s"); + INFO("Escape velocity: " << escape_velocity << " m/s"); + INFO("Circular velocity: " << circular_velocity << " m/s"); + INFO("Semi-major axis: " << a << " m"); + + double expected_v_squared = mu * (2.0 / r - 1.0 / a); + double expected_v = sqrt(expected_v_squared); + + double v_error = fabs(v - expected_v); + double relative_error = v_error / expected_v; + + INFO("Expected velocity: " << expected_v << " m/s"); + INFO("Velocity error: " << v_error << " m/s (" << relative_error * 100.0 << "%)"); + + REQUIRE(relative_error < VELOCITY_TOLERANCE); + REQUIRE(v > escape_velocity * 0.9); + REQUIRE(a < 0.0); + + destroy_simulation(sim); +} + +TEST_CASE("Velocity magnitude accuracy for extreme eccentricities", "[extreme][eccentricity][velocity]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_extreme_eccentricity.toml")); + + CelestialBody* earth = &sim->bodies[0]; + + for (int i = 0; i < sim->craft_count; i++) { + Spacecraft* craft = &sim->spacecraft[i]; + + INFO("Spacecraft " << i << ": e=" << craft->orbit.eccentricity); + + double true_anomalies[] = {0.0, M_PI / 2.0, M_PI, 3.0 * M_PI / 2.0}; + + for (int j = 0; j < 4; j++) { + double nu = true_anomalies[j]; + craft->orbit.true_anomaly = nu; + + Vec3 pos; + Vec3 vel; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos, &vel); + + double r = vec3_magnitude(pos); + double v = vec3_magnitude(vel); + + double a = craft->orbit.semi_major_axis; + double mu = G * earth->mass; + + double expected_v_squared = mu * (2.0 / r - 1.0 / a); + + if (expected_v_squared > 0.0) { + double expected_v = sqrt(expected_v_squared); + double v_error = fabs(v - expected_v); + double relative_error = v_error / expected_v; + + INFO(" ν=" << nu << " rad: v=" << v << " m/s, error=" << relative_error * 100.0 << "%"); + + REQUIRE(relative_error < VELOCITY_TOLERANCE * 10.0); + } + } + } + + destroy_simulation(sim); +} + +TEST_CASE("Period calculation (or lack thereof) for e≥1", "[extreme][eccentricity][period]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_extreme_eccentricity.toml")); + + Spacecraft* high_e = &sim->spacecraft[0]; + Spacecraft* near_parabolic = &sim->spacecraft[1]; + Spacecraft* hyperbolic = &sim->spacecraft[2]; + + double a_e = high_e->orbit.semi_major_axis; + double a_near = near_parabolic->orbit.semi_major_axis; + double a_h = hyperbolic->orbit.semi_major_axis; + + INFO("Highly eccentric (e=0.99): a=" << a_e << " m"); + INFO("Near-parabolic (e=0.9999): a=" << a_near << " m"); + INFO("Hyperbolic (e=1.0001): a=" << a_h << " m"); + + REQUIRE(a_e > 0.0); + REQUIRE(a_near > 0.0); + REQUIRE(a_h < 0.0); + + destroy_simulation(sim); +} diff --git a/tests/test_extreme_eccentricity.toml b/tests/test_extreme_eccentricity.toml new file mode 100644 index 0000000..e2c14ec --- /dev/null +++ b/tests/test_extreme_eccentricity.toml @@ -0,0 +1,53 @@ +# Test Configuration: Extreme Eccentricity Orbits +# Tests near-parabolic and hyperbolic orbits + +[[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 = "Highly_Elliptical" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = 6.5e8, + eccentricity = 0.99, + true_anomaly = 0.0, + inclination = 0.0, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +} + + [[spacecraft]] +name = "Near_Parabolic" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = 1.0e9, + eccentricity = 0.95, + true_anomaly = 0.0, + inclination = 0.0, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +} + + [[spacecraft]] +name = "Slightly_Hyperbolic" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = -1.0e7, + eccentricity = 1.5, + true_anomaly = 0.0, + inclination = 0.0, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +} diff --git a/tests/test_newton_raphson_convergence.cpp b/tests/test_newton_raphson_convergence.cpp new file mode 100644 index 0000000..e0e4eaf --- /dev/null +++ b/tests/test_newton_raphson_convergence.cpp @@ -0,0 +1,233 @@ +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include +#include + +const double CONVERGENCE_TOLERANCE = 1.0e-10; +const int MAX_ITERATIONS = 50; + +TEST_CASE("Newton-Raphson solver - very low eccentricity (e < 0.01)", "[newton][raphson][low_e]") { + const double eccentricities[] = {0.001, 0.01}; + + for (int i = 0; i < 2; i++) { + double e = eccentricities[i]; + INFO("Testing eccentricity: " << e); + + double mean_anomaly = M_PI / 2.0; + + double eccentric_anomaly = solve_kepler_equation(mean_anomaly, e); + double expected_eccentric_anomaly = mean_anomaly; + + double error = fabs(eccentric_anomaly - expected_eccentric_anomaly); + INFO("Eccentric anomaly: " << eccentric_anomaly << " rad"); + INFO("Expected: " << expected_eccentric_anomaly << " rad"); + INFO("Error: " << error); + + REQUIRE(error < 1.0e-6); + } +} + +TEST_CASE("Newton-Raphson solver - moderate eccentricity (0.1 < e < 0.5)", "[newton][raphson][moderate_e]") { + const double eccentricities[] = {0.1, 0.3, 0.5}; + + for (int i = 0; i < 3; i++) { + double e = eccentricities[i]; + INFO("Testing eccentricity: " << e); + + double mean_anomaly = M_PI / 4.0; + + double eccentric_anomaly = solve_kepler_equation(mean_anomaly, e); + + double rhs = mean_anomaly + e * sin(eccentric_anomaly); + double residual = eccentric_anomaly - rhs; + + INFO("Eccentric anomaly: " << eccentric_anomaly << " rad"); + INFO("Residual E - (M + e*sin(E)): " << residual); + + REQUIRE(fabs(residual) < CONVERGENCE_TOLERANCE); + } +} + +TEST_CASE("Newton-Raphson solver - high eccentricity (0.9 < e < 0.99)", "[newton][raphson][high_e]") { + const double eccentricities[] = {0.9, 0.95, 0.99}; + + for (int i = 0; i < 3; i++) { + double e = eccentricities[i]; + INFO("Testing eccentricity: " << e); + + double mean_anomaly = M_PI / 2.0; + + int iterations = 0; + double E = get_initial_trial_value(mean_anomaly, e); + double E_prev = E + 2.0 * CONVERGENCE_TOLERANCE; + + while (fabs(E - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) { + E_prev = E; + double sin_E = sin(E); + E = E - (E - e * sin_E - mean_anomaly) / (1.0 - e * cos(E)); + iterations++; + } + + INFO("Converged in " << iterations << " iterations"); + INFO("Eccentric anomaly: " << E << " rad"); + + double rhs = mean_anomaly + e * sin(E); + double residual = E - rhs; + + INFO("Residual E - (M + e*sin(E)): " << residual); + + REQUIRE(iterations < MAX_ITERATIONS); + REQUIRE(fabs(residual) < CONVERGENCE_TOLERANCE); + } +} + +TEST_CASE("Newton-Raphson solver - mean anomaly near π (worst case)", "[newton][raphson][near_pi]") { + const double eccentricity = 0.7; + const double mean_anomalies[] = {M_PI - 0.01, M_PI, M_PI + 0.01}; + + for (int i = 0; i < 3; i++) { + double M = mean_anomalies[i]; + INFO("Testing mean anomaly: " << M << " rad (" << (M * 180.0 / M_PI) << "°)"); + + int iterations = 0; + double E = get_initial_trial_value(M, eccentricity); + double E_prev = E + 2.0 * CONVERGENCE_TOLERANCE; + + while (fabs(E - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) { + E_prev = E; + double sin_E = sin(E); + E = E - (E - eccentricity * sin_E - M) / (1.0 - eccentricity * cos(E)); + iterations++; + } + + INFO("Converged in " << iterations << " iterations"); + + double rhs = M + eccentricity * sin(E); + double residual = E - rhs; + + INFO("Residual E - (M + e*sin(E)): " << residual); + + REQUIRE(iterations < MAX_ITERATIONS); + REQUIRE(fabs(residual) < CONVERGENCE_TOLERANCE); + } +} + +TEST_CASE("Newton-Raphson solver - large mean anomaly values (M > 1000)", "[newton][raphson][large_M]") { + const double eccentricity = 0.3; + const double mean_anomalies[] = {1000.0, 10000.0}; + + for (int i = 0; i < 2; i++) { + double M = mean_anomalies[i]; + INFO("Testing mean anomaly: " << M << " rad"); + + int iterations = 0; + double E = get_initial_trial_value(M, eccentricity); + double E_prev = E + 2.0 * CONVERGENCE_TOLERANCE; + + while (fabs(E - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) { + E_prev = E; + double sin_E = sin(E); + E = E - (E - eccentricity * sin_E - M) / (1.0 - eccentricity * cos(E)); + iterations++; + } + + INFO("Converged in " << iterations << " iterations"); + + double rhs = M + eccentricity * sin(E); + double residual = E - rhs; + + INFO("Residual E - (M + e*sin(E)): " << residual); + + REQUIRE(iterations < MAX_ITERATIONS); + REQUIRE(fabs(residual) < CONVERGENCE_TOLERANCE); + + double M_reduced = fmod(E - eccentricity * sin(E), 2.0 * M_PI); + double M_target = fmod(M, 2.0 * M_PI); + double angle_diff = fabs(M_reduced - M_target); + + if (angle_diff > M_PI) { + angle_diff = 2.0 * M_PI - angle_diff; + } + + INFO("Reduced mean anomaly: " << M_reduced << " rad"); + INFO("Target reduced: " << M_target << " rad"); + INFO("Angle difference: " << angle_diff << " rad"); + + REQUIRE(angle_diff < CONVERGENCE_TOLERANCE * 10.0); + } +} + +TEST_CASE("Newton-Raphson solver - eccentricity at boundaries (e ≈ 1.0)", "[newton][raphson][boundary]") { + const double eccentricities[] = {0.9999, 1.0001}; + + for (int i = 0; i < 2; i++) { + double e = eccentricities[i]; + INFO("Testing eccentricity: " << e); + + double M = M_PI / 4.0; + + int iterations = 0; + double E = get_initial_trial_value(M, e); + double E_prev = E + 2.0 * CONVERGENCE_TOLERANCE; + + while (fabs(E - E_prev) > CONVERGENCE_TOLERANCE && iterations < MAX_ITERATIONS) { + E_prev = E; + double sin_E = sin(E); + E = E - (E - e * sin_E - M) / (1.0 - e * cos(E)); + iterations++; + } + + INFO("Converged in " << iterations << " iterations"); + + if (fabs(1.0 - e * cos(E)) > 1.0e-10) { + double rhs = M + e * sin(E); + double residual = E - rhs; + + INFO("Residual E - (M + e*sin(E)): " << residual); + + REQUIRE(fabs(residual) < CONVERGENCE_TOLERANCE); + } + } +} + +TEST_CASE("Newton-Raphson solver convergence rate", "[newton][raphson][convergence_rate]") { + const double eccentricity = 0.8; + const double mean_anomaly = M_PI / 3.0; + + double E = get_initial_trial_value(mean_anomaly, eccentricity); + + INFO("Initial guess: " << E << " rad"); + + double previous_residual = std::numeric_limits::max(); + int iteration = 0; + int convergence_count = 0; + + for (int i = 0; i < 10; i++) { + double sin_E = sin(E); + double rhs = mean_anomaly + eccentricity * sin_E; + double residual = fabs(E - rhs); + + INFO("Iteration " << i << ": E = " << E << ", residual = " << residual); + + if (residual < CONVERGENCE_TOLERANCE) { + INFO("Converged at iteration " << i); + break; + } + + if (i > 0 && residual < previous_residual * 0.5) { + convergence_count++; + } + + previous_residual = residual; + E = E - (E - eccentricity * sin_E - mean_anomaly) / (1.0 - eccentricity * cos(E)); + iteration++; + } + + double convergence_ratio = (double)convergence_count / (double)iteration; + INFO("Quadratic convergence ratio: " << convergence_ratio * 100.0 << "%"); + + REQUIRE(convergence_ratio > 0.6); +} diff --git a/tests/test_precision_boundaries.cpp b/tests/test_precision_boundaries.cpp new file mode 100644 index 0000000..9614be4 --- /dev/null +++ b/tests/test_precision_boundaries.cpp @@ -0,0 +1,274 @@ +#include +#include "../src/physics.h" +#include "../src/orbital_mechanics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include +#include + +const double ELEMENT_TOLERANCE = 1.0e-6; +const double VELOCITY_TOLERANCE = 1.0e-3; + +TEST_CASE("Perfect circle (e=0)", "[precision][boundary][circle]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_precision_boundaries.toml")); + + Spacecraft* circle = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + INFO("Testing circular orbit: e=" << circle->orbit.eccentricity); + + Vec3 pos_1; + Vec3 vel_1; + orbital_elements_to_cartesian(circle->orbit, earth->mass, &pos_1, &vel_1); + + double r_1 = vec3_magnitude(pos_1); + double v_1 = vec3_magnitude(vel_1); + + INFO("Radius: " << r_1 << " m"); + INFO("Velocity: " << v_1 << " m/s"); + + double expected_r = circle->orbit.semi_major_axis; + double mu = G * earth->mass; + double expected_v = sqrt(mu / expected_r); + + double r_error = fabs(r_1 - expected_r); + double v_error = fabs(v_1 - expected_v); + + INFO("Expected radius: " << expected_r << " m"); + INFO("Expected velocity: " << expected_v << " m/s"); + INFO("Radius error: " << r_error << " m"); + INFO("Velocity error: " << v_error << " m/s"); + + REQUIRE(r_error < fabs(expected_r) * ELEMENT_TOLERANCE); + REQUIRE(v_error < VELOCITY_TOLERANCE); + + double vis_viva = sqrt(mu * (2.0 / r_1 - 1.0 / circle->orbit.semi_major_axis)); + double vis_viva_error = fabs(v_1 - vis_viva); + + INFO("Vis-viva velocity: " << vis_viva << " m/s"); + INFO("Vis-viva error: " << vis_viva_error << " m/s"); + + REQUIRE(vis_viva_error < VELOCITY_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Polar orbit (i=π/2)", "[precision][boundary][polar]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_precision_boundaries.toml")); + + Spacecraft* polar = &sim->spacecraft[1]; + CelestialBody* earth = &sim->bodies[0]; + + INFO("Testing polar orbit: i=" << polar->orbit.inclination << " rad (" << polar->orbit.inclination * 180.0 / M_PI << "°)"); + + Vec3 pos; + Vec3 vel; + orbital_elements_to_cartesian(polar->orbit, earth->mass, &pos, &vel); + + INFO("Position: (" << pos.x << ", " << pos.y << ", " << pos.z << ") m"); + INFO("Velocity: (" << vel.x << ", " << vel.y << ", " << vel.z << ") m/s"); + + double r = vec3_magnitude(pos); + double v = vec3_magnitude(vel); + + double expected_r = polar->orbit.semi_major_axis * (1.0 - polar->orbit.eccentricity * polar->orbit.eccentricity) / (1.0 + polar->orbit.eccentricity); + + INFO("Expected radius: " << expected_r << " m"); + INFO("Actual radius: " << r << " m"); + + double r_error = fabs(r - expected_r); + + REQUIRE(r_error < fabs(expected_r) * ELEMENT_TOLERANCE); + + double z_expected = r * sin(polar->orbit.inclination); + double z_actual = pos.z; + + INFO("Expected Z: " << z_expected << " m"); + INFO("Actual Z: " << z_actual << " m"); + + double z_error = fabs(z_actual - z_expected); + + REQUIRE(z_error < fabs(expected_r) * ELEMENT_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Retrograde orbit (i=π)", "[precision][boundary][retrograde]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_precision_boundaries.toml")); + + Spacecraft* retrograde = &sim->spacecraft[2]; + CelestialBody* earth = &sim->bodies[0]; + + INFO("Testing retrograde orbit: i=" << retrograde->orbit.inclination << " rad (" << retrograde->orbit.inclination * 180.0 / M_PI << "°)"); + + Vec3 pos; + Vec3 vel; + orbital_elements_to_cartesian(retrograde->orbit, earth->mass, &pos, &vel); + + double r = vec3_magnitude(pos); + double v = vec3_magnitude(vel); + + INFO("Radius: " << r << " m"); + INFO("Velocity: " << v << " m/s"); + + double expected_r = retrograde->orbit.semi_major_axis * (1.0 - retrograde->orbit.eccentricity * retrograde->orbit.eccentricity) / (1.0 + retrograde->orbit.eccentricity); + double mu = G * earth->mass; + double expected_v = sqrt(mu * (2.0 / r - 1.0 / retrograde->orbit.semi_major_axis)); + + double r_error = fabs(r - expected_r); + double v_error = fabs(v - expected_v); + + INFO("Expected radius: " << expected_r << " m"); + INFO("Expected velocity: " << expected_v << " m/s"); + INFO("Radius error: " << r_error << " m"); + INFO("Velocity error: " << v_error << " m/s"); + + REQUIRE(r_error < fabs(expected_r) * ELEMENT_TOLERANCE); + REQUIRE(v_error < VELOCITY_TOLERANCE); + + destroy_simulation(sim); +} + +TEST_CASE("Inclination at 0°, 90°, 180°", "[precision][boundary][inclination]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_precision_boundaries.toml")); + + double expected_inclinations[] = {0.0, M_PI / 2.0, M_PI}; + + for (int i = 0; i < 3; i++) { + Spacecraft* craft = &sim->spacecraft[i]; + double expected_i = expected_inclinations[i]; + + INFO("Spacecraft " << i << ": i=" << craft->orbit.inclination << " rad (" << craft->orbit.inclination * 180.0 / M_PI << "°)"); + + double i_error = fabs(craft->orbit.inclination - expected_i); + + INFO(" Expected inclination: " << expected_i << " rad"); + INFO(" Inclination error: " << i_error << " rad"); + + REQUIRE(i_error < ELEMENT_TOLERANCE); + } + + destroy_simulation(sim); +} + +TEST_CASE("Semi-major axis sign change", "[precision][boundary][semi_major]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_precision_boundaries.toml")); + + for (int i = 0; i < 3; i++) { + Spacecraft* craft = &sim->spacecraft[i]; + double e = craft->orbit.eccentricity; + double a = craft->orbit.semi_major_axis; + + INFO("Spacecraft " << i << ": e=" << e << ", a=" << a << " m"); + + if (e < 1.0) { + INFO(" Elliptical orbit: a > 0"); + REQUIRE(a > 0.0); + } else if (fabs(e - 1.0) < 1.0e-6) { + INFO(" Parabolic orbit: near-circular"); + } else { + INFO(" Hyperbolic orbit: a < 0"); + REQUIRE(a < 0.0); + } + } + + destroy_simulation(sim); +} + +TEST_CASE("Angular momentum conservation", "[precision][boundary][angular_momentum]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_precision_boundaries.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + double true_anomalies[] = {0.0, M_PI / 4.0, M_PI / 2.0, 3.0 * M_PI / 4.0, M_PI}; + + Vec3 initial_h = {0.0, 0.0, 0.0}; + + for (int i = 0; i < 5; i++) { + double nu = true_anomalies[i]; + craft->orbit.true_anomaly = nu; + + Vec3 pos; + Vec3 vel; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos, &vel); + + Vec3 h = vec3_cross(pos, vel); + double h_mag = vec3_magnitude(h); + + INFO("ν=" << nu << " rad: |h|=" << h_mag << " m²/s"); + + if (i == 0) { + initial_h = h; + } else { + double h_error = vec3_distance(h, initial_h); + double relative_error = h_error / h_mag; + + INFO(" Angular momentum error: " << h_error << " m²/s (" << relative_error * 100.0 << "%)"); + + REQUIRE(relative_error < ELEMENT_TOLERANCE); + } + } + + destroy_simulation(sim); +} + +TEST_CASE("Zero/very small radius or velocity", "[precision][boundary][zero]") { + const double TIME_STEP = 60.0; + + SimulationState* sim = create_simulation(10, 3, 0, TIME_STEP); + + REQUIRE(load_system_config(sim, "tests/test_precision_boundaries.toml")); + + Spacecraft* craft = &sim->spacecraft[0]; + CelestialBody* earth = &sim->bodies[0]; + + Vec3 pos; + Vec3 vel; + orbital_elements_to_cartesian(craft->orbit, earth->mass, &pos, &vel); + + double r = vec3_magnitude(pos); + double v = vec3_magnitude(vel); + + INFO("Radius: " << r << " m"); + INFO("Velocity: " << v << " m/s"); + + REQUIRE(r > earth->radius); + REQUIRE(v > 0.0); + + Vec3 r_vec = vec3_normalize(pos); + Vec3 v_vec = vec3_normalize(vel); + + double r_dot_v = vec3_dot(r_vec, v_vec); + + INFO("r̂ · v̂: " << r_dot_v); + + REQUIRE(r_dot_v > -1.0); + REQUIRE(r_dot_v < 1.0); + + destroy_simulation(sim); +} diff --git a/tests/test_precision_boundaries.toml b/tests/test_precision_boundaries.toml new file mode 100644 index 0000000..75fc753 --- /dev/null +++ b/tests/test_precision_boundaries.toml @@ -0,0 +1,53 @@ +# Test Configuration: Boundary Value Cases +# Tests exact boundary values for orbital parameters + +[[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 = "Perfect_Circle" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = 1.0e7, + eccentricity = 0.0, + true_anomaly = 0.0, + inclination = 0.0, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +} + + [[spacecraft]] +name = "Polar_Orbit" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = 1.5e7, + eccentricity = 0.5, + true_anomaly = 0.0, + inclination = 1.57079633, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +} + +[[spacecraft]] +name = "Retrograde_Orbit" +mass = 1000.0 +parent_index = 0 +orbit = { + semi_major_axis = 1.5e7, + eccentricity = 0.5, + true_anomaly = 0.0, + inclination = 3.14159265, + longitude_of_ascending_node = 0.0, + argument_of_periapsis = 0.0 +}