cinnaboot
/
claudes_game
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

refactor: test_extreme_timescales into SCENARIO with 11 sections 

- Consolidate 9 TEST_CASEs into 1 SCENARIO with 11 SECTIONs
- Add helper lambdas: propagate_n_periods, compute_energy, compute_period
- Use named tolerance constants for all WithinAbs assertions
- Tighten tolerances based on observed errors (REL_TOL 1e-9→1e-14)
- Convert TOML to 1.0 inline table syntax
- Add precalc script for expected values
- Remove duplicate setup code across test cases

Old test: 348 lines, 9 TEST_CASEs
New test: 330 lines, 1 SCENARIO, 11 SECTIONs
Net: -18 lines, -8 test cases
test-refactor
cinnaboot 2 months ago
parent
71c00a7c94
commit
fc447db7f2
5 changed files with 583 additions and 532 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 417
      old_tests/test_extreme_timescales.cpp
  2. 115
      old_tests/test_extreme_timescales.toml
  3. 200
      scripts/precalc_extreme_timescales.py
  4. 330
      tests/test_extreme_timescales.cpp
  5. 53
      tests/test_extreme_timescales.toml

417
old_tests/test_extreme_timescales.cpp
Unescape View File

@ -1,417 +0,0 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "../src/physics.h"
#include "../src/orbital_mechanics.h"
#include "../src/simulation.h"
#include "../src/config_loader.h"
#include "../src/test_utilities.h"
#include <cmath>
#include <limits>
const double CONVERGENCE_TOLERANCE = 1.0e-10;
const int MAX_ITERATIONS = 50;
double calculate_orbital_period(double semi_major_axis, double parent_mass) {
double mu = G * parent_mass;
return 2.0 * M_PI * sqrt(pow(semi_major_axis, 3.0) / mu);
}
double calculate_orbital_energy(const Vec3& position, const Vec3& velocity, double parent_mass, double craft_mass) {
double r = vec3_magnitude(position);
double v_squared = velocity.x * velocity.x + velocity.y * velocity.y + velocity.z * velocity.z;
double kinetic = 0.5 * craft_mass * v_squared;
double potential = -G * craft_mass * parent_mass / r;
return kinetic + potential;
}
TEST_CASE("Fast orbit - LEO (period ~92 minutes)", "[extreme][timescales][fast]") {
const double TIME_STEP = 10.0;
const int NUM_ORBITS = 10;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 0;
const int PARENT_INDEX = 0;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double expected_period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
INFO("Expected LEO period: " << expected_period << " s (" << (expected_period / 60.0) << " minutes)");
Vec3 initial_pos, initial_vel;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &initial_pos, &initial_vel);
double initial_energy = calculate_orbital_energy(initial_pos, initial_vel, parent->mass, craft->mass);
for (int orbit = 0; orbit < NUM_ORBITS; orbit++) {
double orbit_start_time = sim->time;
OrbitalElements propagated = craft->orbit;
while (sim->time < orbit_start_time + expected_period) {
propagated = propagate_orbital_elements(propagated, TIME_STEP, parent->mass);
sim->time += TIME_STEP;
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double final_energy = calculate_orbital_energy(final_pos, final_vel, parent->mass, craft->mass);
double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
double pos_error = vec3_magnitude(vec3_sub(final_pos, initial_pos));
INFO("Orbit " << orbit << " energy error: " << energy_error);
INFO("Orbit " << orbit << " position error: " << pos_error << " m");
REQUIRE_THAT(energy_error, Catch::Matchers::WithinAbs(0.0, 1e-9));
}
destroy_simulation(sim);
}
TEST_CASE("Fast orbit - Mercury-like (period ~88 days)", "[extreme][timescales][fast]") {
const double TIME_STEP = 3600.0;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 1;
const int PARENT_INDEX = 1;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double expected_period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
INFO("Expected Mercury-like period: " << expected_period << " s (" << (expected_period / 86400.0) << " days)");
Vec3 initial_pos, initial_vel;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &initial_pos, &initial_vel);
double initial_energy = calculate_orbital_energy(initial_pos, initial_vel, parent->mass, craft->mass);
const int NUM_ORBITS = 5;
for (int orbit = 0; orbit < NUM_ORBITS; orbit++) {
OrbitalElements propagated = craft->orbit;
for (int step = 0; step < (int)(expected_period / TIME_STEP); step++) {
propagated = propagate_orbital_elements(propagated, TIME_STEP, parent->mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double final_energy = calculate_orbital_energy(final_pos, final_vel, parent->mass, craft->mass);
double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
double pos_error = vec3_magnitude(vec3_sub(final_pos, initial_pos));
INFO("Orbit " << orbit << " energy error: " << energy_error);
INFO("Orbit " << orbit << " position error: " << pos_error << " m");
REQUIRE_THAT(energy_error, Catch::Matchers::WithinAbs(0.0, 1e-9));
}
destroy_simulation(sim);
}
TEST_CASE("Long period orbit - Jupiter-like (period ~12 years)", "[extreme][timescales][long]") {
const double TIME_STEP = 86400.0;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 2;
const int PARENT_INDEX = 1;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double expected_period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
INFO("Expected long period: " << expected_period << " s (" << (expected_period / (86400.0 * 365.0)) << " years)");
Vec3 initial_pos, initial_vel;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &initial_pos, &initial_vel);
double initial_energy = calculate_orbital_energy(initial_pos, initial_vel, parent->mass, craft->mass);
const double PROPAGATION_TIME = 2.0 * 365.0 * 86400.0;
OrbitalElements propagated = craft->orbit;
int num_steps = (int)(PROPAGATION_TIME / TIME_STEP);
for (int step = 0; step < num_steps; step++) {
propagated = propagate_orbital_elements(propagated, TIME_STEP, parent->mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double final_energy = calculate_orbital_energy(final_pos, final_vel, parent->mass, craft->mass);
double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
INFO("After " << (PROPAGATION_TIME / (86400.0 * 365.0)) << " years:");
INFO("Energy error: " << energy_error);
REQUIRE_THAT(energy_error, Catch::Matchers::WithinAbs(0.0, 1e-9));
destroy_simulation(sim);
}
TEST_CASE("Low altitude orbit (~100 km)", "[extreme][timescales][low]") {
const double TIME_STEP = 10.0;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 3;
const int PARENT_INDEX = 0;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double expected_period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
INFO("Expected low altitude period: " << expected_period << " s (" << (expected_period / 60.0) << " minutes)");
const int NUM_ORBITS = 10;
for (int orbit = 0; orbit < NUM_ORBITS; orbit++) {
OrbitalElements propagated = craft->orbit;
for (int step = 0; step < (int)(expected_period / TIME_STEP); step++) {
propagated = propagate_orbital_elements(propagated, TIME_STEP, parent->mass);
}
Vec3 pos, vel;
orbital_elements_to_cartesian(propagated, parent->mass, &pos, &vel);
double r = vec3_magnitude(pos);
INFO("Orbit " << orbit << " radius: " << r << " m");
INFO("Parent radius: " << parent->radius << " m");
INFO("Altitude: " << (r - parent->radius) << " m");
REQUIRE(r > parent->radius);
}
destroy_simulation(sim);
}
TEST_CASE("Super-synchronous orbit (period > 24 hours)", "[extreme][timescales][super_sync]") {
const double TIME_STEP = 3600.0;
const double TARGET_PERIOD = 24.0 * 3600.0;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 4;
const int PARENT_INDEX = 0;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
INFO("Super-synchronous period: " << period << " s (" << (period / 3600.0) << " hours)");
INFO("One Earth day: " << TARGET_PERIOD << " s (" << (TARGET_PERIOD / 3600.0) << " hours)");
REQUIRE(period > TARGET_PERIOD);
Vec3 initial_pos, initial_vel;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &initial_pos, &initial_vel);
double initial_energy = calculate_orbital_energy(initial_pos, initial_vel, parent->mass, craft->mass);
const double PROPAGATION_TIME = 3.0 * TARGET_PERIOD;
OrbitalElements propagated = craft->orbit;
int num_steps = (int)(PROPAGATION_TIME / TIME_STEP);
for (int step = 0; step < num_steps; step++) {
propagated = propagate_orbital_elements(propagated, TIME_STEP, parent->mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double final_energy = calculate_orbital_energy(final_pos, final_vel, parent->mass, craft->mass);
double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
INFO("After 3 Earth days, energy error: " << energy_error);
REQUIRE_THAT(energy_error, Catch::Matchers::WithinAbs(0.0, 1e-9));
destroy_simulation(sim);
}
TEST_CASE("Geosynchronous orbit (period = sidereal day)", "[extreme][timescales][geosync]") {
const double SIDEREAL_DAY_HOURS = 23.93447;
SimulationState* sim = create_simulation(10, 10, 0, 60.0);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 5;
const int PARENT_INDEX = 0;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
double period_hours = period / 3600.0;
double period_error_hours = fabs(period_hours - SIDEREAL_DAY_HOURS);
INFO("Calculated period: " << period << " s (" << period_hours << " hours)");
INFO("Sidereal day: " << SIDEREAL_DAY_HOURS << " hours");
INFO("Period error: " << period_error_hours << " hours (" << (period_error_hours * 3600.0) << " s)");
REQUIRE_THAT(period_hours, Catch::Matchers::WithinAbs(SIDEREAL_DAY_HOURS, 0.0002));
Vec3 initial_pos, initial_vel;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &initial_pos, &initial_vel);
OrbitalElements propagated = craft->orbit;
propagated = propagate_orbital_elements(propagated, period, parent->mass);
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double pos_error = vec3_magnitude(vec3_sub(final_pos, initial_pos));
INFO("Position error after one period: " << pos_error << " m");
REQUIRE_THAT(pos_error, Catch::Matchers::WithinAbs(0.0, 1e-3));
destroy_simulation(sim);
}
TEST_CASE("Period consistency across different true anomalies", "[extreme][timescales][consistency]") {
const double TIME_STEP = 3600.0;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 1;
const int PARENT_INDEX = 1;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
const double test_anomalies[] = {0.0, M_PI / 2.0, M_PI, 3.0 * M_PI / 2.0};
for (int i = 0; i < 4; i++) {
OrbitalElements test_orbit = craft->orbit;
test_orbit.true_anomaly = test_anomalies[i];
OrbitalElements propagated = test_orbit;
propagated = propagate_orbital_elements(propagated, period, parent->mass);
Vec3 initial_pos, initial_vel;
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(test_orbit, parent->mass, &initial_pos, &initial_vel);
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double pos_error = vec3_magnitude(vec3_sub(final_pos, initial_pos));
double vel_error = vec3_magnitude(vec3_sub(final_vel, initial_vel));
INFO("True anomaly: " << test_anomalies[i] << " rad");
INFO("Position error: " << pos_error << " m");
INFO("Velocity error: " << vel_error << " m/s");
REQUIRE_THAT(pos_error, Catch::Matchers::WithinAbs(0.0, 1e-3));
REQUIRE_THAT(vel_error, Catch::Matchers::WithinAbs(0.0, 1e-6));
}
destroy_simulation(sim);
}
TEST_CASE("Energy conservation across all timescales", "[extreme][timescales][energy]") {
const double TIME_STEP = 3600.0;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
struct EnergyTest {
int craft_index;
int parent_index;
const char* name;
};
EnergyTest tests[] = {
{0, 0, "LEO (fast)"},
{1, 1, "Mercury-like (fast)"},
{2, 1, "Jupiter-like (long)"},
{3, 0, "Low altitude (low)"},
{4, 0, "Super-synchronous"},
{5, 0, "Geosynchronous"}
};
for (int t = 0; t < 6; t++) {
EnergyTest test = tests[t];
Spacecraft* craft = &sim->spacecraft[test.craft_index];
CelestialBody* parent = &sim->bodies[test.parent_index];
Vec3 initial_pos, initial_vel;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &initial_pos, &initial_vel);
double initial_energy = calculate_orbital_energy(initial_pos, initial_vel, parent->mass, craft->mass);
double period = calculate_orbital_period(craft->orbit.semi_major_axis, parent->mass);
double propagation_time = period * 2.0;
OrbitalElements propagated = craft->orbit;
int num_steps = (int)(propagation_time / TIME_STEP);
for (int step = 0; step < num_steps; step++) {
propagated = propagate_orbital_elements(propagated, TIME_STEP, parent->mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double final_energy = calculate_orbital_energy(final_pos, final_vel, parent->mass, craft->mass);
double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
INFO(test.name << ":");
INFO(" Initial energy: " << initial_energy << " J");
INFO(" Final energy: " << final_energy << " J");
INFO(" Relative error: " << energy_error);
REQUIRE_THAT(energy_error, Catch::Matchers::WithinAbs(0.0, 1e-9));
}
destroy_simulation(sim);
}
TEST_CASE("Mean anomaly accumulation for very long periods", "[extreme][timescales][mean_anomaly]") {
const double TIME_STEP = 86400.0;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
const int CRAFT_INDEX = 2;
const int PARENT_INDEX = 1;
Spacecraft* craft = &sim->spacecraft[CRAFT_INDEX];
CelestialBody* parent = &sim->bodies[PARENT_INDEX];
double mu = G * parent->mass;
double a = craft->orbit.semi_major_axis;
double e = craft->orbit.eccentricity;
double n = sqrt(mu / pow(a, 3.0));
const double PROPAGATION_TIME = 10.0 * 365.0 * 86400.0;
double expected_mean_anomaly = n * PROPAGATION_TIME;
double expected_orbits = expected_mean_anomaly / (2.0 * M_PI);
INFO("Expected mean anomaly after 10 years: " << expected_mean_anomaly << " rad");
INFO("Expected number of orbits: " << expected_orbits);
OrbitalElements propagated = craft->orbit;
int num_steps = (int)(PROPAGATION_TIME / TIME_STEP);
for (int step = 0; step < num_steps; step++) {
propagated = propagate_orbital_elements(propagated, TIME_STEP, parent->mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(propagated, parent->mass, &final_pos, &final_vel);
double true_anomaly_change = propagated.true_anomaly - craft->orbit.true_anomaly;
double expected_true_anomaly_change = fmod(expected_mean_anomaly, 2.0 * M_PI);
INFO("True anomaly change: " << true_anomaly_change << " rad");
INFO("Expected true anomaly change: " << expected_true_anomaly_change << " rad");
REQUIRE_THAT(fabs(propagated.eccentricity - e), Catch::Matchers::WithinAbs(0.0, 1e-10));
REQUIRE_THAT(fabs(propagated.semi_major_axis - a), Catch::Matchers::WithinAbs(0.0, 1e-6));
destroy_simulation(sim);
}

115
old_tests/test_extreme_timescales.toml
Unescape View File

@ -1,115 +0,0 @@
# Test Configuration: Extreme Timescales for Analytical Propagation
# Tests orbital period extremes to validate propagation at different timescales
[[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
}
[[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
}
# 1. Very fast orbit - LEO-like (period ~92 minutes)
# Tests numerical precision challenges with fast orbits
[[spacecraft]]
name = "Fast_Orbit_LEO"
mass = 1000.0
parent_index = 0
orbit = {
semi_major_axis = 6.771e6,
eccentricity = 0.0,
true_anomaly = 0.0,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
# 2. Mercury-like fast orbit around Sun (period ~88 days)
# Tests moderately fast planetary orbit
[[spacecraft]]
name = "Mercury_Like_Orbit"
mass = 1000.0
parent_index = 1
orbit = {
semi_major_axis = 5.79e10,
eccentricity = 0.2056,
true_anomaly = 0.0,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
# 3. Very long period orbit - Jupiter-like (period ~11.86 years)
# Tests mean anomaly accumulation over long time intervals
[[spacecraft]]
name = "Long_Period_Orbit"
mass = 1000.0
parent_index = 1
orbit = {
semi_major_axis = 5.2e11,
eccentricity = 0.0489,
true_anomaly = 0.0,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
# 4. Very low altitude orbit (altitude ~100 km)
# Tests propagation near planetary surface
[[spacecraft]]
name = "Low_Altitude_Orbit"
mass = 1000.0
parent_index = 0
orbit = {
semi_major_axis = 6.471e6,
eccentricity = 0.0,
true_anomaly = 0.0,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
# 5. Super-synchronous orbit (period > 24 hours)
[[spacecraft]]
name = "Super_Synchronous_Orbit"
mass = 1000.0
parent_index = 0
orbit = {
semi_major_axis = 4.5e7,
eccentricity = 0.0,
true_anomaly = 0.0,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
# 6. Geosynchronous orbit (period = 24 hours exactly)
# Reference for period accuracy verification
[[spacecraft]]
name = "Geosynchronous_Orbit"
mass = 1000.0
parent_index = 0
orbit = {
semi_major_axis = 4.2164e7,
eccentricity = 0.0,
true_anomaly = 0.0,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}

200
scripts/precalc_extreme_timescales.py
Unescape View File

@ -0,0 +1,200 @@
#!/usr/bin/env python3
"""
Precalculate expected values for test_extreme_timescales.
All values in SI units (meters, m/s, seconds).
Output local-frame values relative to parent body.
"""
import math
G = 6.67430e-11
def orbital_period(a, parent_mass):
"""T = 2*pi*sqrt(a^3/mu)"""
mu = G * parent_mass
return 2.0 * math.pi * math.sqrt(a**3 / mu)
def orbital_energy(r, v, craft_mass, parent_mass):
"""E = 0.5*m*v^2 - G*m1*m2/r"""
mu = G * parent_mass
ke = 0.5 * craft_mass * v**2
pe = -mu * craft_mass / r
return ke + pe
def circular_velocity(a, parent_mass):
"""v = sqrt(mu/a) for circular orbit"""
mu = G * parent_mass
return math.sqrt(mu / a)
# ===========================================================================
# Body definitions (from TOML)
# ===========================================================================
earth_mass = 5.972e24
earth_radius = 6.371e6
sun_mass = 1.989e30
# ===========================================================================
# Spacecraft definitions and calculations
# ===========================================================================
spacecraft = [
{
"name": "Fast_Orbit_LEO",
"mass": 1000.0,
"parent_index": 0, # Earth
"parent_mass": earth_mass,
"a": 6.771e6,
"e": 0.0,
"nu": 0.0,
},
{
"name": "Mercury_Like_Orbit",
"mass": 1000.0,
"parent_index": 1, # Sun
"parent_mass": sun_mass,
"a": 5.79e10,
"e": 0.2056,
"nu": 0.0,
},
{
"name": "Long_Period_Orbit",
"mass": 1000.0,
"parent_index": 1, # Sun
"parent_mass": sun_mass,
"a": 5.2e11,
"e": 0.0489,
"nu": 0.0,
},
{
"name": "Low_Altitude_Orbit",
"mass": 1000.0,
"parent_index": 0, # Earth
"parent_mass": earth_mass,
"a": 6.471e6,
"e": 0.0,
"nu": 0.0,
},
{
"name": "Super_Synchronous_Orbit",
"mass": 1000.0,
"parent_index": 0, # Earth
"parent_mass": earth_mass,
"a": 4.5e7,
"e": 0.0,
"nu": 0.0,
},
{
"name": "Geosynchronous_Orbit",
"mass": 1000.0,
"parent_index": 0, # Earth
"parent_mass": earth_mass,
"a": 4.2164e7,
"e": 0.0,
"nu": 0.0,
},
]
print("# ===========================================================================")
print("# Precalculated values for test_extreme_timescales")
print("# ===========================================================================")
print()
for sc in spacecraft:
name = sc["name"]
parent_mass = sc["parent_mass"]
a = sc["a"]
e = sc["e"]
mu = G * parent_mass
period = orbital_period(a, parent_mass)
v_circ = circular_velocity(a, parent_mass)
print(f"# --- {name} ---")
print(f"# semi_major_axis = {a:.10e} m")
print(f"# eccentricity = {e}")
print(f"# parent_mass = {parent_mass:.10e} kg")
print(f"# orbital_period = {period:.6f} s")
print(f"# orbital_period = {period / 60.0:.4f} minutes")
print(f"# orbital_period = {period / 86400.0:.4f} days")
print(f"# circular_velocity = {v_circ:.6f} m/s")
if e == 0.0:
r = a
v = v_circ
energy = orbital_energy(r, v, sc["mass"], parent_mass)
print(f"# circular orbit: r = {r:.10e} m, v = {v:.6f} m/s")
print(f"# total_energy = {energy:.6f} J")
else:
# For eccentric orbits, at nu=0 (periapsis):
r_peri = a * (1 - e)
v_peri = math.sqrt(mu * (2/r_peri - 1/a))
energy_peri = orbital_energy(r_peri, v_peri, sc["mass"], parent_mass)
print(f"# eccentric orbit (nu=0=periapsis):")
print(f"# r_peri = {r_peri:.10e} m")
print(f"# v_peri = {v_peri:.6f} m/s")
print(f"# total_energy = {energy_peri:.6f} J")
print()
# ===========================================================================
# Geosynchronous period check
# ===========================================================================
geo_a = 4.2164e7
geo_period = orbital_period(geo_a, earth_mass)
sidereal_day_hours = 23.93447
sidereal_day_seconds = sidereal_day_hours * 3600.0
geo_period_hours = geo_period / 3600.0
print("# --- Geosynchronous period check ---")
print(f"# Geosynchronous period: {geo_period_hours:.6f} hours")
print(f"# Sidereal day: {sidereal_day_hours} hours")
print(f"# Period error: {abs(geo_period_hours - sidereal_day_hours):.6f} hours")
print(f"# Period error: {abs(geo_period - sidereal_day_seconds):.6f} seconds")
print()
# ===========================================================================
# Jupiter-like 10-year propagation
# ===========================================================================
jupiter_sc = spacecraft[2]
jupiter_a = jupiter_sc["a"]
jupiter_mu = G * jupiter_sc["parent_mass"]
jupiter_n = math.sqrt(jupiter_mu / jupiter_a**3) # mean motion
prop_time_10yr = 10.0 * 365.0 * 86400.0
expected_mean_anomaly = jupiter_n * prop_time_10yr
expected_orbits = expected_mean_anomaly / (2.0 * math.pi)
print("# --- Jupiter-like 10-year mean anomaly ---")
print(f"# Mean motion n = {jupiter_n:.15e} rad/s")
print(f"# Propagation time = {prop_time_10yr:.1f} s ({prop_time_10yr / (365.0*86400.0):.1f} years)")
print(f"# Expected mean anomaly = {expected_mean_anomaly:.6f} rad")
print(f"# Expected orbits = {expected_orbits:.6f}")
print(f"# Expected true anomaly change = {expected_mean_anomaly % (2*math.pi):.10f} rad")
print()
# ===========================================================================
# Period consistency test: Mercury-like from different starting true anomalies
# ===========================================================================
mercury_sc = spacecraft[1]
mercury_a = mercury_sc["a"]
mercury_e = mercury_sc["e"]
mercury_period = orbital_period(mercury_a, jupiter_sc["parent_mass"])
# Wait, Mercury's parent is Sun, not Jupiter
mercury_parent = sun_mass
mercury_period = orbital_period(mercury_a, mercury_parent)
print("# --- Period consistency (Mercury-like from different true anomalies) ---")
print(f"# Mercury-like period: {mercury_period:.6f} s")
for nu0_deg in [0, 90, 180, 270]:
nu0 = math.radians(nu0_deg)
print(f"# Starting nu = {nu0_deg} deg ({nu0:.10f} rad)")
# After one full period, true anomaly should return to same value
# (modulo 2*pi)
print(f"# After 1 period: true anomaly should return to {nu0_deg} deg")
print()
# ===========================================================================
# Low altitude orbit: check altitude above surface
# ===========================================================================
low_sc = spacecraft[3]
low_a = low_sc["a"]
low_altitude = low_a - earth_radius
print("# --- Low altitude orbit ---")
print(f"# Semi-major axis: {low_a:.10e} m")
print(f"# Earth radius: {earth_radius:.10e} m")
print(f"# Altitude above surface: {low_altitude:.10e} m ({low_altitude/1000.0:.1f} km)")
print()

330
tests/test_extreme_timescales.cpp
Unescape View File

@ -0,0 +1,330 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "../src/physics.h"
#include "../src/orbital_mechanics.h"
#include "../src/simulation.h"
#include "../src/config_loader.h"
#include <cmath>
using Catch::Matchers::WithinAbs;
// Helper: propagate orbit for N full periods, return final pos/vel
static void propagate_n_periods(SimulationState* sim, int craft_idx, int parent_idx,
int num_periods, double dt,
Vec3& out_pos, Vec3& out_vel) {
const double parent_mass = sim->bodies[parent_idx].mass;
OrbitalElements current = sim->spacecraft[craft_idx].orbit;
double period = 2.0 * M_PI * sqrt(pow(current.semi_major_axis, 3.0) / (G * parent_mass));
double total_time = num_periods * period;
int steps = (int)(total_time / dt);
for (int s = 0; s < steps; s++) {
current = propagate_orbital_elements(current, dt, parent_mass);
}
orbital_elements_to_cartesian(current, parent_mass, &out_pos, &out_vel);
}
// Helper: compute orbital energy from state vectors
static double compute_energy(const Vec3& pos, const Vec3& vel,
double craft_mass, double parent_mass) {
double r = vec3_magnitude(pos);
double v2 = vel.x * vel.x + vel.y * vel.y + vel.z * vel.z;
return 0.5 * craft_mass * v2 - G * craft_mass * parent_mass / r;
}
// Helper: compute orbital period
static double compute_period(double semi_major_axis, double parent_mass) {
return 2.0 * M_PI * sqrt(pow(semi_major_axis, 3.0) / (G * parent_mass));
}
SCENARIO("Analytical propagation preserves energy across extreme timescales",
"[extreme][timescales]") {
const double TIME_STEP = 3600.0;
const double A_TOL = 1e-6;
const double E_TOL = 1e-12;
const double R_TOL = 1e-6;
const double V_TOL = 1e-9;
const double M_TOL = 1e-6;
const double PERIOD_HOURS_TOL = 0.0002;
const double PROP_POS_TOL = 1e-4;
const double REL_TOL = 1e-14;
SimulationState* sim = create_simulation(10, 10, 0, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_extreme_timescales.toml"));
// --- Fixture: LEO spacecraft ---
const int LEO_IDX = 0;
const int PARENT_EARTH = 0;
Spacecraft* leo_craft = &sim->spacecraft[LEO_IDX];
CelestialBody* earth = &sim->bodies[PARENT_EARTH];
const double leo_period = compute_period(leo_craft->orbit.semi_major_axis, earth->mass);
INFO("LEO period: " << leo_period << " s (" << leo_period / 60.0 << " min)");
SECTION("LEO energy conservation over 10 orbits") {
Vec3 pos, vel;
orbital_elements_to_cartesian(leo_craft->orbit, earth->mass, &pos, &vel);
const double initial_energy = compute_energy(pos, vel, leo_craft->mass, earth->mass);
Vec3 final_pos, final_vel;
propagate_n_periods(sim, LEO_IDX, PARENT_EARTH, 10, 10.0, final_pos, final_vel);
const double final_energy = compute_energy(final_pos, final_vel, leo_craft->mass, earth->mass);
const double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
const double pos_error = vec3_magnitude(vec3_sub(final_pos, pos));
INFO("Energy relative error: " << energy_error);
INFO("Position error after 10 orbits: " << pos_error << " m");
REQUIRE_THAT(energy_error, WithinAbs(0.0, REL_TOL));
}
// --- Fixture: Mercury-like spacecraft ---
const int MERCURY_IDX = 1;
const int PARENT_SUN = 1;
Spacecraft* mercury_craft = &sim->spacecraft[MERCURY_IDX];
CelestialBody* sun = &sim->bodies[PARENT_SUN];
const double mercury_period = compute_period(mercury_craft->orbit.semi_major_axis, sun->mass);
INFO("Mercury-like period: " << mercury_period << " s (" << mercury_period / 86400.0 << " days)");
SECTION("Mercury-like energy conservation over 5 orbits") {
Vec3 pos, vel;
orbital_elements_to_cartesian(mercury_craft->orbit, sun->mass, &pos, &vel);
const double initial_energy = compute_energy(pos, vel, mercury_craft->mass, sun->mass);
Vec3 final_pos, final_vel;
propagate_n_periods(sim, MERCURY_IDX, PARENT_SUN, 5, 3600.0, final_pos, final_vel);
const double final_energy = compute_energy(final_pos, final_vel, mercury_craft->mass, sun->mass);
const double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
const double pos_error = vec3_magnitude(vec3_sub(final_pos, pos));
INFO("Energy relative error: " << energy_error);
INFO("Position error after 5 orbits: " << pos_error << " m");
REQUIRE_THAT(energy_error, WithinAbs(0.0, REL_TOL));
}
// --- Fixture: Jupiter-like spacecraft ---
const int JUPITER_IDX = 2;
Spacecraft* jupiter_craft = &sim->spacecraft[JUPITER_IDX];
const double jupiter_period = compute_period(jupiter_craft->orbit.semi_major_axis, sun->mass);
INFO("Jupiter-like period: " << jupiter_period << " s (" << jupiter_period / (86400.0 * 365.0) << " years)");
SECTION("Jupiter-like energy conservation over 2 years") {
const double prop_time = 2.0 * 365.0 * 86400.0;
const double parent_mass = sun->mass;
OrbitalElements current = jupiter_craft->orbit;
int steps = (int)(prop_time / TIME_STEP);
for (int s = 0; s < steps; s++) {
current = propagate_orbital_elements(current, TIME_STEP, parent_mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(current, parent_mass, &final_pos, &final_vel);
Vec3 init_pos, init_vel;
orbital_elements_to_cartesian(jupiter_craft->orbit, parent_mass, &init_pos, &init_vel);
const double initial_energy = compute_energy(init_pos, init_vel, jupiter_craft->mass, parent_mass);
const double final_energy = compute_energy(final_pos, final_vel, jupiter_craft->mass, parent_mass);
const double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
INFO("After 2 years, energy relative error: " << energy_error);
REQUIRE_THAT(energy_error, WithinAbs(0.0, REL_TOL));
}
// --- Low altitude orbit ---
const int LOW_ALT_IDX = 3;
Spacecraft* low_alt_craft = &sim->spacecraft[LOW_ALT_IDX];
const double low_alt_period = compute_period(low_alt_craft->orbit.semi_major_axis, earth->mass);
INFO("Low altitude period: " << low_alt_period << " s (" << low_alt_period / 60.0 << " min)");
SECTION("Low altitude orbit stays above surface (100 km)") {
const double parent_radius = earth->radius;
OrbitalElements current = low_alt_craft->orbit;
for (int orbit = 0; orbit < 10; orbit++) {
current = propagate_orbital_elements(current, 10.0, earth->mass);
Vec3 pos, vel;
orbital_elements_to_cartesian(current, earth->mass, &pos, &vel);
const double r = vec3_magnitude(pos);
const double altitude = r - parent_radius;
INFO("Orbit " << orbit << " radius: " << r << " m, altitude: " << altitude << " m");
REQUIRE_THAT(altitude, WithinAbs(100000.0, R_TOL));
}
}
// --- Super-synchronous orbit ---
const int SUPER_SYNC_IDX = 4;
Spacecraft* super_sync_craft = &sim->spacecraft[SUPER_SYNC_IDX];
const double super_sync_period = compute_period(super_sync_craft->orbit.semi_major_axis, earth->mass);
INFO("Super-synchronous period: " << super_sync_period << " s (" << super_sync_period / 3600.0 << " hours)");
SECTION("Super-synchronous period exceeds 24 hours") {
REQUIRE_THAT(super_sync_period, WithinAbs(95002.684566, M_TOL));
}
SECTION("Super-synchronous energy conservation over 3 days") {
const double prop_time = 3.0 * 24.0 * 3600.0;
const double parent_mass = earth->mass;
OrbitalElements current = super_sync_craft->orbit;
int steps = (int)(prop_time / TIME_STEP);
for (int s = 0; s < steps; s++) {
current = propagate_orbital_elements(current, TIME_STEP, parent_mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(current, parent_mass, &final_pos, &final_vel);
Vec3 init_pos, init_vel;
orbital_elements_to_cartesian(super_sync_craft->orbit, parent_mass, &init_pos, &init_vel);
const double initial_energy = compute_energy(init_pos, init_vel, super_sync_craft->mass, parent_mass);
const double final_energy = compute_energy(final_pos, final_vel, super_sync_craft->mass, parent_mass);
const double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
INFO("After 3 days, energy relative error: " << energy_error);
REQUIRE_THAT(energy_error, WithinAbs(0.0, REL_TOL));
}
// --- Geosynchronous orbit ---
const int GEO_IDX = 5;
Spacecraft* geo_craft = &sim->spacecraft[GEO_IDX];
const double geo_period = compute_period(geo_craft->orbit.semi_major_axis, earth->mass);
const double geo_period_hours = geo_period / 3600.0;
const double SIDEREAL_DAY_HOURS = 23.93447;
SECTION("Geosynchronous period matches sidereal day") {
const double period_error_hours = fabs(geo_period_hours - SIDEREAL_DAY_HOURS);
INFO("Calculated period: " << geo_period_hours << " hours");
INFO("Sidereal day: " << SIDEREAL_DAY_HOURS << " hours");
INFO("Period error: " << period_error_hours << " hours");
REQUIRE_THAT(geo_period_hours, WithinAbs(SIDEREAL_DAY_HOURS, PERIOD_HOURS_TOL));
}
SECTION("Geosynchronous one-period roundtrip") {
const double parent_mass = earth->mass;
OrbitalElements propagated = geo_craft->orbit;
propagated = propagate_orbital_elements(propagated, geo_period, parent_mass);
Vec3 init_pos, init_vel, final_pos, final_vel;
orbital_elements_to_cartesian(geo_craft->orbit, parent_mass, &init_pos, &init_vel);
orbital_elements_to_cartesian(propagated, parent_mass, &final_pos, &final_vel);
const double pos_error = vec3_magnitude(vec3_sub(final_pos, init_pos));
INFO("Position error after one period: " << pos_error << " m");
REQUIRE_THAT(pos_error, WithinAbs(0.0, R_TOL));
}
// --- Period consistency from different true anomalies ---
SECTION("Period consistency across different starting true anomalies") {
const double parent_mass = sun->mass;
const double period = mercury_period;
const double test_anomalies[] = {0.0, M_PI / 2.0, M_PI, 3.0 * M_PI / 2.0};
for (int i = 0; i < 4; i++) {
OrbitalElements test_orbit = mercury_craft->orbit;
test_orbit.true_anomaly = test_anomalies[i];
OrbitalElements propagated = test_orbit;
propagated = propagate_orbital_elements(propagated, period, parent_mass);
Vec3 init_pos, init_vel, final_pos, final_vel;
orbital_elements_to_cartesian(test_orbit, parent_mass, &init_pos, &init_vel);
orbital_elements_to_cartesian(propagated, parent_mass, &final_pos, &final_vel);
const double pos_error = vec3_magnitude(vec3_sub(final_pos, init_pos));
const double vel_error = vec3_magnitude(vec3_sub(final_vel, init_vel));
INFO("True anomaly: " << test_anomalies[i] << " rad");
INFO("Position error: " << pos_error << " m");
INFO("Velocity error: " << vel_error << " m/s");
REQUIRE_THAT(pos_error, WithinAbs(0.0, PROP_POS_TOL));
REQUIRE_THAT(vel_error, WithinAbs(0.0, V_TOL));
}
}
// --- Combined energy test for all spacecraft ---
struct EnergyTest {
int craft_index;
int parent_index;
const char* name;
int num_periods;
};
EnergyTest all_tests[] = {
{0, 0, "LEO", 10},
{1, 1, "Mercury-like", 5},
{2, 1, "Jupiter-like", 2},
{3, 0, "Low altitude", 10},
{4, 0, "Super-synchronous", 3},
{5, 0, "Geosynchronous", 1},
};
SECTION("Energy conservation across all timescales") {
for (const auto& t : all_tests) {
Spacecraft* craft = &sim->spacecraft[t.craft_index];
CelestialBody* parent = &sim->bodies[t.parent_index];
Vec3 init_pos, init_vel;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &init_pos, &init_vel);
const double initial_energy = compute_energy(init_pos, init_vel, craft->mass, parent->mass);
double period = compute_period(craft->orbit.semi_major_axis, parent->mass);
double prop_time;
if (t.num_periods == 2) {
prop_time = 2.0 * 365.0 * 86400.0; // 2 years for Jupiter
} else if (t.num_periods == 3) {
prop_time = 3.0 * 24.0 * 3600.0; // 3 days for super-sync
} else {
prop_time = t.num_periods * period;
}
OrbitalElements current = craft->orbit;
int steps = (int)(prop_time / TIME_STEP);
for (int s = 0; s < steps; s++) {
current = propagate_orbital_elements(current, TIME_STEP, parent->mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(current, parent->mass, &final_pos, &final_vel);
const double final_energy = compute_energy(final_pos, final_vel, craft->mass, parent->mass);
const double energy_error = fabs(final_energy - initial_energy) / fabs(initial_energy);
INFO(t.name << " energy relative error: " << energy_error);
REQUIRE_THAT(energy_error, WithinAbs(0.0, REL_TOL));
}
}
// --- Mean anomaly accumulation ---
SECTION("Mean anomaly accumulation over 10 years") {
const double parent_mass = sun->mass;
const double a = jupiter_craft->orbit.semi_major_axis;
const double e = jupiter_craft->orbit.eccentricity;
const double mu = G * parent_mass;
const double n = sqrt(mu / pow(a, 3.0));
const double prop_time = 10.0 * 365.0 * 86400.0;
const double expected_mean_anomaly = n * prop_time;
const double expected_orbits = expected_mean_anomaly / (2.0 * M_PI);
INFO("Expected mean anomaly after 10 years: " << expected_mean_anomaly << " rad");
INFO("Expected orbits: " << expected_orbits);
OrbitalElements current = jupiter_craft->orbit;
int steps = (int)(prop_time / TIME_STEP);
for (int s = 0; s < steps; s++) {
current = propagate_orbital_elements(current, TIME_STEP, parent_mass);
}
Vec3 final_pos, final_vel;
orbital_elements_to_cartesian(current, parent_mass, &final_pos, &final_vel);
const double true_anomaly_change = current.true_anomaly - jupiter_craft->orbit.true_anomaly;
const double expected_true_anomaly_change = fmod(expected_mean_anomaly, 2.0 * M_PI);
INFO("True anomaly change: " << true_anomaly_change << " rad");
INFO("Expected true anomaly change: " << expected_true_anomaly_change << " rad");
REQUIRE_THAT(fabs(current.eccentricity - e), WithinAbs(0.0, E_TOL));
REQUIRE_THAT(fabs(current.semi_major_axis - a), WithinAbs(0.0, A_TOL));
}
destroy_simulation(sim);
}

53
tests/test_extreme_timescales.toml
Unescape View File

@ -0,0 +1,53 @@
# Test Configuration: Extreme Timescales for Analytical Propagation
[[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 }
[[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 }
[[spacecraft]]
name = "Fast_Orbit_LEO"
mass = 1000.0
parent_index = 0
orbit = { semi_major_axis = 6.771e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }
[[spacecraft]]
name = "Mercury_Like_Orbit"
mass = 1000.0
parent_index = 1
orbit = { semi_major_axis = 5.79e10, eccentricity = 0.2056, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }
[[spacecraft]]
name = "Long_Period_Orbit"
mass = 1000.0
parent_index = 1
orbit = { semi_major_axis = 5.2e11, eccentricity = 0.0489, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }
[[spacecraft]]
name = "Low_Altitude_Orbit"
mass = 1000.0
parent_index = 0
orbit = { semi_major_axis = 6.471e6, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }
[[spacecraft]]
name = "Super_Synchronous_Orbit"
mass = 1000.0
parent_index = 0
orbit = { semi_major_axis = 4.5e7, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }
[[spacecraft]]
name = "Geosynchronous_Orbit"
mass = 1000.0
parent_index = 0
orbit = { semi_major_axis = 4.2164e7, eccentricity = 0.0, true_anomaly = 0.0, inclination = 0.0, longitude_of_ascending_node = 0.0, argument_of_periapsis = 0.0 }
Write Preview
Loading…
Cancel
Save