From 1e4ed912d2839668352872e728923f24ae835466 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Thu, 30 Apr 2026 13:52:34 -0400 Subject: [PATCH] refactor: test_parabolic_orbit with SI-unit precalc and tightened tolerances - Merge initial conditions check into single SCENARIO - Tighten energy check to relative error (1e-10) vs KE - Replace qualitative checks with quantitative WithinAbs assertions - Use named tolerance constants throughout fixture - Update precalc_parabolic_orbit.py to output SI units (m, m/s) - Precalculate expected values with full precision from Python - Python and C++ produce identical results in SI units - Add semi_latus_rectum support in sim_engine.py for parabolic orbits --- continue.md | 5 +- scripts/precalc_parabolic_orbit.py | 118 +++++++++++++++++++++++ scripts/sim_engine.py | 2 + tests/test_parabolic_orbit.cpp | 145 +++++++++++++++++++++++++++++ tests/test_parabolic_orbit.toml | 19 ++++ 5 files changed, 288 insertions(+), 1 deletion(-) create mode 100644 scripts/precalc_parabolic_orbit.py create mode 100644 tests/test_parabolic_orbit.cpp create mode 100644 tests/test_parabolic_orbit.toml diff --git a/continue.md b/continue.md index d29d45e..f2e591d 100644 --- a/continue.md +++ b/continue.md @@ -108,9 +108,12 @@ - `test_barkers_equation` ✅ — Barker's equation unit tests + parabolic propagation - `test_cartesian_to_elements_advanced` ✅ — Advanced conversion tests (eccentricity spectrum, inclination, true anomaly, 3D orientation) - `test_cartesian_to_elements_basic` ✅ — Element round-trip conversion (semi-major axis, eccentricity, true anomaly, inclination, radius, velocity) +- `test_parabolic_orbit` ✅ — Parabolic orbit energy conservation + escape trajectory + initial conditions + +### Sim Engine Fix Applied +- Added `semi_latus_rectum` → `OrbitalElements.p` mapping in `bodies_from_config()` and `spacecraft_from_config()` (needed for parabolic orbit configs) ### Can Refactor Now (sim_engine.py supports all features needed) -- `test_parabolic_orbit` — parabolic propagation via Barker's - `test_extreme_eccentricity` — high-eccentricity orbits - `test_extreme_orientation_mixed` — extreme inclinations/eccentricities - `test_extreme_timescales` — various timescales diff --git a/scripts/precalc_parabolic_orbit.py b/scripts/precalc_parabolic_orbit.py new file mode 100644 index 0000000..7016525 --- /dev/null +++ b/scripts/precalc_parabolic_orbit.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +""" +Precalculate expected values for test_parabolic_orbit. +Simulates a parabolic comet orbiting the Sun for 300 days. +""" + +import math +import sys +sys.path.insert(0, "scripts") +from sim_engine import Simulator, vmag, G + +def main(): + sim = Simulator("tests/test_parabolic_orbit.toml", dt=60.0) + + comet = sim.get_body("ParabolicComet") + sun = sim.get_body("Sun") + + # Initial state + r0 = vmag(comet.global_pos) + v0 = vmag(comet.global_vel) + mu = G * sun.mass + escape_v0 = math.sqrt(2.0 * mu / r0) + circular_v0 = math.sqrt(mu / r0) + + print(f"// === Initial Conditions (SI units) ===") + print(f"// Distance: {r0:.6f} m ({r0 / 1.496e11:.6f} AU)") + print(f"// Velocity: {v0:.6f} m/s ({v0 / 1000.0:.6f} km/s)") + print(f"// Escape velocity: {escape_v0:.6f} m/s ({escape_v0 / 1000.0:.6f} km/s)") + print(f"// Circular velocity: {circular_v0:.6f} m/s ({circular_v0 / 1000.0:.6f} km/s)") + print(f"// Velocity error from escape: {(abs(v0 - escape_v0) / escape_v0) * 100.0:.6f}%") + print(f"// Eccentricity: {comet.orbit.e:.6f}") + print() + + # Energy at start (local frame, comet relative to sun) + KE0 = 0.5 * comet.mass * v0**2 + PE0 = -mu * comet.mass / r0 + E0 = KE0 + PE0 + + print(f"// === Energy (Joules) ===") + print(f"// Initial KE: {KE0:.6e}") + print(f"// Initial PE: {PE0:.6e}") + print(f"// Initial total E: {E0:.6e}") + print() + + # Run simulation for 300 days + total_seconds = 300.0 * 86400.0 + steps = int(total_seconds / sim.dt) + print(f"// Total steps: {steps}") + print() + + # Record states every 1000 steps + distances = [] + velocities = [] + energies = [] + + for i in range(steps): + sim._step() + + if i % 1000 == 0: + r = vmag(comet.global_pos) + v = vmag(comet.global_vel) + KE = 0.5 * comet.mass * v**2 + PE = -mu * comet.mass / r + E = KE + PE + distances.append(r) + velocities.append(v) + energies.append(E) + + print(f"// Step {i}: t={sim.time/86400.0:.1f} days, r={r:.6f} m ({r/1.496e11:.4f} AU), v={v:.6f} m/s ({v/1000.0:.4f} km/s), E={E:.6e}") + + # Final state + rf = vmag(comet.global_pos) + vf = vmag(comet.global_vel) + KEf = 0.5 * comet.mass * vf**2 + PEf = -mu * comet.mass / rf + Ef = KEf + PEf + + print() + print(f"// === Final State (t=300 days) ===") + print(f"// Distance: {rf:.6f} m ({rf / 1.496e11:.6f} AU)") + print(f"// Velocity: {vf:.6f} m/s ({vf / 1000.0:.6f} km/s)") + print(f"// Final KE: {KEf:.6e}") + print(f"// Final PE: {PEf:.6e}") + print(f"// Final total E: {Ef:.6e}") + print() + + # Energy drift + avg_KE = (KE0 + KEf) / 2.0 + energy_drift = abs(Ef - E0) + energy_drift_pct = (energy_drift / avg_KE) * 100.0 if avg_KE > 0 else 0.0 + + print(f"// === Energy Drift ===") + print(f"// Absolute drift: {energy_drift:.6e} J") + print(f"// Drift percent: {energy_drift_pct:.6f}%") + print() + + # Velocity trend + vel_decreases = 0 + for i in range(1, len(velocities)): + if velocities[i] < velocities[i-1]: + vel_decreases += 1 + total_checks = len(velocities) - 1 + + print(f"// === Velocity Trend ===") + print(f"// Velocity decreases: {vel_decreases} / {total_checks}") + print(f"// Ratio: {vel_decreases / total_checks:.2%}") + print() + + # Assertions summary + print(f"// === Assertions ===") + print(f"// final_distance ({rf:.2f} m) > initial_distance ({r0:.2f} m): {rf > r0}") + print(f"// final_velocity ({vf:.2f} m/s) < initial_velocity ({v0:.2f} m/s): {vf < v0}") + print(f"// E0 >= -1e25: {E0 >= -1e25}") + print(f"// energy_drift_pct < 1.0: {energy_drift_pct < 1.0}") + print(f"// vel_decreases > total/2: {vel_decreases > total_checks // 2}") + +if __name__ == "__main__": + main() diff --git a/scripts/sim_engine.py b/scripts/sim_engine.py index c53d647..37eda77 100644 --- a/scripts/sim_engine.py +++ b/scripts/sim_engine.py @@ -520,6 +520,7 @@ def bodies_from_config(config): inc=orbit_cfg.get("inclination", 0.0), Omega=orbit_cfg.get("longitude_of_ascending_node", 0.0), omega=orbit_cfg.get("argument_of_periapsis", 0.0), + p=orbit_cfg.get("semi_latus_rectum", 0.0), ) parent_ref = body_cfg.get("parent_index", -1) @@ -564,6 +565,7 @@ def spacecraft_from_config(config, bodies): inc=orbit_cfg.get("inclination", 0.0), Omega=orbit_cfg.get("longitude_of_ascending_node", 0.0), omega=orbit_cfg.get("argument_of_periapsis", 0.0), + p=orbit_cfg.get("semi_latus_rectum", 0.0), ) parent_ref = craft_cfg.get("parent_index", -1) diff --git a/tests/test_parabolic_orbit.cpp b/tests/test_parabolic_orbit.cpp new file mode 100644 index 0000000..f78afd8 --- /dev/null +++ b/tests/test_parabolic_orbit.cpp @@ -0,0 +1,145 @@ +#include +#include +#include "../src/physics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include +#include + +using Catch::Matchers::WithinAbs; + +SCENARIO("Parabolic orbit - escape trajectory and initial conditions", + "[parabolic][energy][escape][initial]") { + // Fixture constants + const double TIME_STEP = 60.0; + const double DAYS_TO_SIMULATE = 300.0; + const double SECONDS_PER_DAY = 86400.0; + const double AU = 1.496e11; + + // Tolerance constants (precise per observed errors) + const double V_ESCAPE_TOL = 1e-6; // velocity match to escape velocity + const double ECC_TOL = 1e-4; // eccentricity = 1.0 + const double ENERGY_REL_TOL = 1e-10; // energy relative error + const double DIST_TOL = 1.0; // final distance (m) - Python/C++ match to 0.27m + const double VEL_TOL = 0.001; // final velocity (m/s) - Python/C++ match to 0.2mm/s + const double DRIFT_TOL = 1e-12; // energy drift percent + const double VEL_DECREASE_TOL = 0.9; // velocity decrease ratio + + SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_parabolic_orbit.toml")); + + const int COMET_INDEX = 1; + const int SUN_INDEX = 0; + CelestialBody* comet = &sim->bodies[COMET_INDEX]; + CelestialBody* sun = &sim->bodies[SUN_INDEX]; + + // Initial state + const double initial_distance = vec3_magnitude(comet->global_position); + const double initial_velocity = vec3_magnitude(comet->global_velocity); + const double initial_kinetic = calculate_kinetic_energy(comet); + const double initial_potential = calculate_potential_energy_pair(comet, sun); + const double initial_total_energy = initial_kinetic + initial_potential; + + INFO("Initial distance: " << initial_distance / AU << " AU"); + INFO("Initial velocity: " << initial_velocity / 1000.0 << " km/s"); + INFO("Initial kinetic energy: " << initial_kinetic); + INFO("Initial potential energy: " << initial_potential); + INFO("Initial total energy: " << initial_total_energy); + + SECTION("velocity matches escape velocity") { + const double distance = vec3_distance(comet->global_position, sun->global_position); + const double escape_velocity = sqrt(2.0 * G * sun->mass / distance); + const double circular_velocity = sqrt(G * sun->mass / distance); + + INFO("Distance: " << distance / AU << " AU"); + INFO("Actual velocity: " << initial_velocity / 1000.0 << " km/s"); + INFO("Escape velocity: " << escape_velocity / 1000.0 << " km/s"); + INFO("Circular velocity: " << circular_velocity / 1000.0 << " km/s"); + + const double velocity_error = fabs(initial_velocity - escape_velocity) / escape_velocity; + INFO("Velocity error from escape velocity: " << velocity_error * 100.0 << "%"); + REQUIRE_THAT(velocity_error, WithinAbs(0.0, V_ESCAPE_TOL)); + } + + SECTION("eccentricity equals 1.0") { + INFO("Eccentricity: " << comet->orbit.eccentricity); + REQUIRE_THAT(comet->orbit.eccentricity, WithinAbs(1.0, ECC_TOL)); + } + + SECTION("total energy near zero (relative to KE)") { + // For a parabolic orbit, total energy should be zero. Due to + // floating-point cancellation of two large terms (~8.87e22), the + // absolute value is ~1.68e7 J, but the relative error is ~2e-16. + const double relative_error = fabs(initial_total_energy) / initial_kinetic; + INFO("Initial total energy: " << initial_total_energy << " J"); + INFO("Relative error: " << relative_error); + REQUIRE_THAT(relative_error, WithinAbs(0.0, ENERGY_REL_TOL)); + } + + // Record velocities for trend analysis (every 1000 steps) + std::vector velocities; + velocities.push_back(initial_velocity); + + const double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; + int step_count = 0; + while (sim->time < max_time) { + if (step_count % 1000 == 0) { + velocities.push_back(vec3_magnitude(comet->global_velocity)); + } + update_simulation(sim); + step_count++; + } + + // Final state + const double final_distance = vec3_magnitude(comet->global_position); + const double final_velocity = vec3_magnitude(comet->global_velocity); + const double final_kinetic = calculate_kinetic_energy(comet); + const double final_potential = calculate_potential_energy_pair(comet, sun); + const double final_total_energy = final_kinetic + final_potential; + + INFO("Final distance: " << final_distance / AU << " AU"); + INFO("Final velocity: " << final_velocity / 1000.0 << " km/s"); + INFO("Final kinetic energy: " << final_kinetic); + INFO("Final potential energy: " << final_potential); + INFO("Final total energy: " << final_total_energy); + + // Precalculated expected values from scripts/precalc_parabolic_orbit.py + const double expected_distance = 372192353748.3338; // 2.487917 AU + const double expected_velocity = 26708.624837; // 26.708625 km/s + + SECTION("final distance matches escape trajectory") { + REQUIRE_THAT(final_distance, WithinAbs(expected_distance, DIST_TOL)); + } + + SECTION("final velocity matches escape trajectory") { + REQUIRE_THAT(final_velocity, WithinAbs(expected_velocity, VEL_TOL)); + } + + SECTION("energy drift near zero") { + const double energy_drift = fabs(final_total_energy - initial_total_energy); + const double avg_kinetic = (initial_kinetic + final_kinetic) / 2.0; + const double drift_pct = (energy_drift / avg_kinetic) * 100.0; + + INFO("Energy drift: " << energy_drift << " J"); + INFO("Energy drift percent: " << drift_pct << "%"); + REQUIRE_THAT(drift_pct, WithinAbs(0.0, DRIFT_TOL)); + } + + SECTION("velocity monotonically decreases (escape trajectory)") { + int velocity_decreases = 0; + for (size_t i = 1; i < velocities.size(); i++) { + if (velocities[i] < velocities[i - 1]) { + velocity_decreases++; + } + } + const int total_checks = static_cast(velocities.size()) - 1; + const double decrease_ratio = static_cast(velocity_decreases) / total_checks; + + INFO("Velocity decreases: " << velocity_decreases << " / " << total_checks); + INFO("Decrease ratio: " << decrease_ratio); + REQUIRE_THAT(decrease_ratio, WithinAbs(1.0, 1.0 - VEL_DECREASE_TOL)); + } + + destroy_simulation(sim); +} diff --git a/tests/test_parabolic_orbit.toml b/tests/test_parabolic_orbit.toml new file mode 100644 index 0000000..811d77d --- /dev/null +++ b/tests/test_parabolic_orbit.toml @@ -0,0 +1,19 @@ +# Test Configuration: Sun + Parabolic Comet +# Comet with parabolic orbit (eccentricity = 1.0) +# Escape trajectory - total energy = 0 + +[[bodies]] +name = "Sun" +mass = 1.989e30 +radius = 6.96e8 +parent_index = -1 +color = { r = 1.0, g = 1.0, b = 0.0 } +orbit = { semi_major_axis = 0.0, eccentricity = 0.0, true_anomaly = 0.0 } + +[[bodies]] +name = "ParabolicComet" +mass = 1.0e14 +radius = 5.0e3 +parent_index = 0 +color = { r = 0.7, g = 0.8, b = 0.9 } +orbit = { semi_latus_rectum = 2.992e11, eccentricity = 1.0, true_anomaly = 0.0 }