54 changed files with 301 additions and 521 deletions
@ -0,0 +1,27 @@ |
|||||||
|
# Test Configuration: Sun + Earth (circular orbit) |
||||||
|
# Earth at 1 AU with circular orbit |
||||||
|
# Expected orbital period: ~365 days |
||||||
|
|
||||||
|
[[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 = "Earth" |
||||||
|
mass = 5.972e24 |
||||||
|
radius = 6.371e6 |
||||||
|
parent_index = 0 |
||||||
|
color = { r = 0.0, g = 0.5, b = 1.0 } |
||||||
|
orbit = { |
||||||
|
semi_major_axis = 1.496e11, |
||||||
|
eccentricity = 0.0, |
||||||
|
true_anomaly = 0.0 |
||||||
|
} |
||||||
@ -0,0 +1,102 @@ |
|||||||
|
#include <catch2/catch_test_macros.hpp> |
||||||
|
#include "../src/physics.h" |
||||||
|
#include "../src/simulation.h" |
||||||
|
#include "../src/config_loader.h" |
||||||
|
#include "../src/test_utilities.h" |
||||||
|
#include <cmath> |
||||||
|
|
||||||
|
TEST_CASE("Orbital period - Earth (RK4)", "[period][rk4]") { |
||||||
|
const double TIME_STEP = 60.0; |
||||||
|
const double EXPECTED_PERIOD_DAYS = 365.0; |
||||||
|
const double SECONDS_PER_DAY = 86400.0; |
||||||
|
const double MAX_SIMULATION_DAYS = 400.0; |
||||||
|
|
||||||
|
SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP); |
||||||
|
|
||||||
|
REQUIRE(load_system_config(sim, "tests/test_orbital_period.toml")); |
||||||
|
|
||||||
|
OrbitTracker* tracker = create_orbit_tracker(1); |
||||||
|
|
||||||
|
double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY; |
||||||
|
while (sim->time < max_time && !tracker->orbit_completed) { |
||||||
|
update_simulation(sim); |
||||||
|
update_orbit_tracker(tracker, &sim->bodies[1], &sim->bodies[0], sim->time); |
||||||
|
} |
||||||
|
|
||||||
|
REQUIRE(tracker->orbit_completed); |
||||||
|
|
||||||
|
double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY; |
||||||
|
double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS); |
||||||
|
|
||||||
|
INFO("Expected period: " << EXPECTED_PERIOD_DAYS << " days"); |
||||||
|
INFO("Measured period: " << measured_period_days << " days"); |
||||||
|
INFO("Error: " << period_error_days << " days"); |
||||||
|
|
||||||
|
REQUIRE(period_error_days < 5.0); |
||||||
|
|
||||||
|
destroy_orbit_tracker(tracker); |
||||||
|
destroy_simulation(sim); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("Orbital period - Mars (RK4)", "[period][rk4]") { |
||||||
|
const double TIME_STEP = 60.0; |
||||||
|
const double EXPECTED_PERIOD_DAYS = 687.0; |
||||||
|
const double SECONDS_PER_DAY = 86400.0; |
||||||
|
const double MAX_SIMULATION_DAYS = 750.0; |
||||||
|
|
||||||
|
SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP); |
||||||
|
|
||||||
|
REQUIRE(load_system_config(sim, "tests/test_orbital_period.toml")); |
||||||
|
|
||||||
|
OrbitTracker* tracker = create_orbit_tracker(2); |
||||||
|
|
||||||
|
double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY; |
||||||
|
while (sim->time < max_time && !tracker->orbit_completed) { |
||||||
|
update_simulation(sim); |
||||||
|
update_orbit_tracker(tracker, &sim->bodies[2], &sim->bodies[0], sim->time); |
||||||
|
} |
||||||
|
|
||||||
|
REQUIRE(tracker->orbit_completed); |
||||||
|
|
||||||
|
double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY; |
||||||
|
double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS); |
||||||
|
|
||||||
|
INFO("Expected period: " << EXPECTED_PERIOD_DAYS << " days"); |
||||||
|
INFO("Measured period: " << measured_period_days << " days"); |
||||||
|
INFO("Error: " << period_error_days << " days"); |
||||||
|
|
||||||
|
REQUIRE(period_error_days < 25.0); |
||||||
|
|
||||||
|
destroy_orbit_tracker(tracker); |
||||||
|
destroy_simulation(sim); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("Orbit direction - prograde for zero inclination", "[direction]") { |
||||||
|
const double TIME_STEP = 60.0; |
||||||
|
const double TEST_DURATION_DAYS = 1.0; |
||||||
|
const double SECONDS_PER_DAY = 86400.0; |
||||||
|
const int STEPS = (int)(TEST_DURATION_DAYS * SECONDS_PER_DAY / TIME_STEP); |
||||||
|
|
||||||
|
SimulationState* sim = create_simulation(2, 0, 0, TIME_STEP); |
||||||
|
REQUIRE(load_system_config(sim, "tests/test_energy.toml")); |
||||||
|
|
||||||
|
CelestialBody* sun = &sim->bodies[0]; |
||||||
|
CelestialBody* earth = &sim->bodies[1]; |
||||||
|
|
||||||
|
Vec3 initial_rel_pos = vec3_sub(earth->global_position, sun->global_position); |
||||||
|
double theta_start = atan2(initial_rel_pos.y, initial_rel_pos.x); |
||||||
|
|
||||||
|
for (int i = 0; i < STEPS; i++) { |
||||||
|
update_simulation(sim); |
||||||
|
} |
||||||
|
|
||||||
|
Vec3 final_rel_pos = vec3_sub(earth->global_position, sun->global_position); |
||||||
|
double theta_final = atan2(final_rel_pos.y, final_rel_pos.x); |
||||||
|
|
||||||
|
INFO("Initial angle: " << theta_start << " rad"); |
||||||
|
INFO("Final angle: " << theta_final << " rad"); |
||||||
|
|
||||||
|
REQUIRE(theta_final > theta_start); |
||||||
|
|
||||||
|
destroy_simulation(sim); |
||||||
|
} |
||||||
@ -0,0 +1,39 @@ |
|||||||
|
# Test Configuration: Sun + Earth + Mars (circular orbits) |
||||||
|
# Earth at 1 AU, Mars at 1.5 AU with circular orbits |
||||||
|
# Expected orbital periods: Earth ~365 days, Mars ~687 days |
||||||
|
|
||||||
|
[[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 = "Earth" |
||||||
|
mass = 5.972e24 |
||||||
|
radius = 6.371e6 |
||||||
|
parent_index = 0 |
||||||
|
color = { r = 0.0, g = 0.5, b = 1.0 } |
||||||
|
orbit = { |
||||||
|
semi_major_axis = 1.496e11, |
||||||
|
eccentricity = 0.0, |
||||||
|
true_anomaly = 0.0 |
||||||
|
} |
||||||
|
|
||||||
|
[[bodies]] |
||||||
|
name = "Mars" |
||||||
|
mass = 6.39e23 |
||||||
|
radius = 3.3895e6 |
||||||
|
parent_index = 0 |
||||||
|
color = { r = 0.8, g = 0.3, b = 0.1 } |
||||||
|
orbit = { |
||||||
|
semi_major_axis = 2.244e11, |
||||||
|
eccentricity = 0.0, |
||||||
|
true_anomaly = 0.0 |
||||||
|
} |
||||||
@ -0,0 +1,133 @@ |
|||||||
|
#include <catch2/catch_test_macros.hpp> |
||||||
|
#include <catch2/matchers/catch_matchers_floating_point.hpp> |
||||||
|
#include "../src/physics.h" |
||||||
|
#include "../src/orbital_mechanics.h" |
||||||
|
#include <cmath> |
||||||
|
|
||||||
|
TEST_CASE("True anomaly round-trip conversion at periapsis", "[orbital_elements][true_anomaly]") { |
||||||
|
double parent_mass = 5.972e24; |
||||||
|
|
||||||
|
OrbitalElements elements = {0}; |
||||||
|
elements.semi_major_axis = 7000e3; |
||||||
|
elements.eccentricity = 0.3; |
||||||
|
elements.true_anomaly = 0.0; |
||||||
|
elements.inclination = 0.0; |
||||||
|
elements.longitude_of_ascending_node = 0.0; |
||||||
|
elements.argument_of_periapsis = 0.0; |
||||||
|
|
||||||
|
Vec3 pos, vel; |
||||||
|
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
||||||
|
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
||||||
|
|
||||||
|
INFO("Original true_anomaly: " << elements.true_anomaly); |
||||||
|
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
||||||
|
|
||||||
|
REQUIRE_THAT(reconstructed.true_anomaly, |
||||||
|
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("True anomaly round-trip conversion at apoapsis", "[orbital_elements][true_anomaly]") { |
||||||
|
double parent_mass = 5.972e24; |
||||||
|
|
||||||
|
OrbitalElements elements = {0}; |
||||||
|
elements.semi_major_axis = 7000e3; |
||||||
|
elements.eccentricity = 0.3; |
||||||
|
elements.true_anomaly = M_PI; |
||||||
|
elements.inclination = 0.0; |
||||||
|
elements.longitude_of_ascending_node = 0.0; |
||||||
|
elements.argument_of_periapsis = 0.0; |
||||||
|
|
||||||
|
Vec3 pos, vel; |
||||||
|
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
||||||
|
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
||||||
|
|
||||||
|
INFO("Original true_anomaly: " << elements.true_anomaly); |
||||||
|
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
||||||
|
|
||||||
|
REQUIRE_THAT(reconstructed.true_anomaly, |
||||||
|
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("True anomaly round-trip conversion at 90 degrees", "[orbital_elements][true_anomaly]") { |
||||||
|
double parent_mass = 5.972e24; |
||||||
|
|
||||||
|
OrbitalElements elements = {0}; |
||||||
|
elements.semi_major_axis = 7000e3; |
||||||
|
elements.eccentricity = 0.3; |
||||||
|
elements.true_anomaly = M_PI / 2.0; |
||||||
|
elements.inclination = 0.0; |
||||||
|
elements.longitude_of_ascending_node = 0.0; |
||||||
|
elements.argument_of_periapsis = 0.0; |
||||||
|
|
||||||
|
Vec3 pos, vel; |
||||||
|
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
||||||
|
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
||||||
|
|
||||||
|
INFO("Original true_anomaly: " << elements.true_anomaly); |
||||||
|
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
||||||
|
|
||||||
|
REQUIRE_THAT(reconstructed.true_anomaly, |
||||||
|
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("True anomaly round-trip conversion at 270 degrees", "[orbital_elements][true_anomaly]") { |
||||||
|
double parent_mass = 5.972e24; |
||||||
|
|
||||||
|
OrbitalElements elements = {0}; |
||||||
|
elements.semi_major_axis = 7000e3; |
||||||
|
elements.eccentricity = 0.3; |
||||||
|
elements.true_anomaly = 3.0 * M_PI / 2.0; |
||||||
|
elements.inclination = 0.0; |
||||||
|
elements.longitude_of_ascending_node = 0.0; |
||||||
|
elements.argument_of_periapsis = 0.0; |
||||||
|
|
||||||
|
Vec3 pos, vel; |
||||||
|
orbital_elements_to_cartesian(elements, parent_mass, &pos, &vel); |
||||||
|
OrbitalElements reconstructed = cartesian_to_orbital_elements(pos, vel, parent_mass); |
||||||
|
|
||||||
|
INFO("Original true_anomaly: " << elements.true_anomaly); |
||||||
|
INFO("Reconstructed true_anomaly: " << reconstructed.true_anomaly); |
||||||
|
|
||||||
|
REQUIRE_THAT(reconstructed.true_anomaly, |
||||||
|
Catch::Matchers::WithinAbs(elements.true_anomaly, 0.01)); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("Radius at periapsis matches expected value", "[orbital_elements][sanity]") { |
||||||
|
double parent_mass = 5.972e24; |
||||||
|
|
||||||
|
OrbitalElements peri = {0}; |
||||||
|
peri.semi_major_axis = 7000e3; |
||||||
|
peri.eccentricity = 0.3; |
||||||
|
peri.true_anomaly = 0.0; |
||||||
|
|
||||||
|
Vec3 pos, vel; |
||||||
|
orbital_elements_to_cartesian(peri, parent_mass, &pos, &vel); |
||||||
|
double r_peri = vec3_magnitude(pos); |
||||||
|
double expected_peri = peri.semi_major_axis * (1.0 - peri.eccentricity); |
||||||
|
|
||||||
|
INFO("At true_anomaly=0:"); |
||||||
|
INFO(" Calculated radius: " << r_peri); |
||||||
|
INFO(" Expected: " << expected_peri); |
||||||
|
|
||||||
|
REQUIRE_THAT(r_peri, Catch::Matchers::WithinAbs(expected_peri, 1.0)); |
||||||
|
} |
||||||
|
|
||||||
|
TEST_CASE("Radius at apoapsis matches expected value", "[orbital_elements][sanity]") { |
||||||
|
double parent_mass = 5.972e24; |
||||||
|
|
||||||
|
OrbitalElements apo = {0}; |
||||||
|
apo.semi_major_axis = 7000e3; |
||||||
|
apo.eccentricity = 0.3; |
||||||
|
apo.true_anomaly = M_PI; |
||||||
|
|
||||||
|
Vec3 pos, vel; |
||||||
|
orbital_elements_to_cartesian(apo, parent_mass, &pos, &vel); |
||||||
|
double r_apo = vec3_magnitude(pos); |
||||||
|
double expected_apo = apo.semi_major_axis * (1.0 + apo.eccentricity); |
||||||
|
|
||||||
|
INFO("At true_anomaly=pi:"); |
||||||
|
INFO(" Calculated radius: " << r_apo); |
||||||
|
INFO(" Expected: " << expected_apo); |
||||||
|
|
||||||
|
REQUIRE_THAT(r_apo, Catch::Matchers::WithinAbs(expected_apo, 1.0)); |
||||||
|
} |
||||||
@ -1,73 +0,0 @@ |
|||||||
# Makefile for Informational Tests
|
|
||||||
# These tests are not part of the main test suite
|
|
||||||
# They are for diagnostic and informational purposes only
|
|
||||||
|
|
||||||
BUILD_DIR = ../../build
|
|
||||||
SRC_DIR = ../../src
|
|
||||||
|
|
||||||
# Compiler flags
|
|
||||||
CFLAGS = -Wall -Wextra -g -ggdb3 -std=c++14
|
|
||||||
INCLUDES = -I$(SRC_DIR) \
|
|
||||||
-isystem../../ext/tomlc17/src \
|
|
||||||
-isystem../../ext/raylib/src \
|
|
||||||
-isystem../../ext/raygui/src \
|
|
||||||
-isystem../../ext/raygui/styles
|
|
||||||
|
|
||||||
# Source files (from main project that informational tests depend on)
|
|
||||||
CORE_OBJECTS = $(BUILD_DIR)/test_utilities.o \
|
|
||||||
$(BUILD_DIR)/physics.o \
|
|
||||||
$(BUILD_DIR)/orbital_mechanics.o \
|
|
||||||
$(BUILD_DIR)/simulation.o \
|
|
||||||
$(BUILD_DIR)/config_loader.o \
|
|
||||||
$(BUILD_DIR)/config_validator.o \
|
|
||||||
$(BUILD_DIR)/maneuver.o \
|
|
||||||
$(BUILD_DIR)/spacecraft.o \
|
|
||||||
$(BUILD_DIR)/tomlc17.o
|
|
||||||
|
|
||||||
# Informational test sources
|
|
||||||
INFO_TESTS = test_time_step_stability.cpp
|
|
||||||
INFO_BINS = $(patsubst %.cpp,%,$(INFO_TESTS))
|
|
||||||
|
|
||||||
# Default target - build all informational tests
|
|
||||||
all: $(INFO_BINS) |
|
||||||
@echo "All informational tests built successfully"
|
|
||||||
@echo ""
|
|
||||||
@echo "To run tests:"
|
|
||||||
@echo " ./test_time_step_stability"
|
|
||||||
@echo ""
|
|
||||||
@echo "To run specific test case:"
|
|
||||||
@echo " ./test_time_step_stability '[timestep][stability]'"
|
|
||||||
|
|
||||||
# Build individual test binaries
|
|
||||||
$(INFO_BINS): $(patsubst %.cpp,%.o,$(INFO_TESTS)) $(CORE_OBJECTS) |
|
||||||
g++ $(patsubst %.cpp,%.o,$<) $(CORE_OBJECTS) -o $@ -lCatch2Main -lCatch2 -lm
|
|
||||||
@echo "Built $@"
|
|
||||||
|
|
||||||
# Build test object files
|
|
||||||
test_time_step_stability.o: test_time_step_stability.cpp |
|
||||||
g++ $(CFLAGS) $(INCLUDES) -c $< -o $@
|
|
||||||
|
|
||||||
# Build core objects if they don't exist
|
|
||||||
$(BUILD_DIR)/%.o: |
|
||||||
$(MAKE) -C ../../ build/$*.o
|
|
||||||
|
|
||||||
# Clean built files in informational directory
|
|
||||||
clean: |
|
||||||
rm -f $(INFO_BINS) *.o
|
|
||||||
|
|
||||||
# Clean all (including core build - use with caution!)
|
|
||||||
clean-all: |
|
||||||
$(MAKE) -C ../../ clean
|
|
||||||
rm -f $(INFO_BINS) *.o
|
|
||||||
|
|
||||||
# Run the time step stability test (must run from project root)
|
|
||||||
run-timestep: test_time_step_stability |
|
||||||
@echo "Note: Tests must be run from project root directory:"
|
|
||||||
@echo " From project root: ./tests/informational/$@"
|
|
||||||
|
|
||||||
# Run with verbose output to see binary search progress
|
|
||||||
run-timestep-verbose: test_time_step_stability |
|
||||||
@echo "Note: Tests must be run from project root directory:"
|
|
||||||
@echo " From project root: ./tests/informational/$@ -s \"[timestep][stability]\""
|
|
||||||
|
|
||||||
.PHONY: all clean clean-all run-timestep run-timestep-verbose |
|
||||||
@ -1,72 +0,0 @@ |
|||||||
# Informational Tests |
|
||||||
|
|
||||||
This directory contains tests that are not part of the standard test suite. These tests are provided for informational purposes only and should be run separately. |
|
||||||
|
|
||||||
## Purpose |
|
||||||
|
|
||||||
Informational tests are designed to: |
|
||||||
- Explore simulation boundaries and limits |
|
||||||
- Provide diagnostic information about the simulation |
|
||||||
- Test extreme or edge cases that don't fit into normal test suites |
|
||||||
- Help understand system behavior under various conditions |
|
||||||
|
|
||||||
## Running Informational Tests |
|
||||||
|
|
||||||
Informational tests are **not** run by `make test`. To build and run them: |
|
||||||
|
|
||||||
```bash |
|
||||||
cd tests/informational |
|
||||||
make |
|
||||||
./time_step_test |
|
||||||
``` |
|
||||||
|
|
||||||
## Current Tests |
|
||||||
|
|
||||||
### `test_time_step_stability.cpp` |
|
||||||
|
|
||||||
Determines the maximum stable time step for the RK4 integration across different orbital regimes. |
|
||||||
|
|
||||||
**Test Bodies:** |
|
||||||
- **Mercury Orbiter** (MESSENGER-like): 200 km altitude, ~12 hour period |
|
||||||
- **Io** (Jupiter's moon): 421,700 km orbit, ~1.77 day period |
|
||||||
- **Moon** (Earth's moon): 384,400 km orbit, ~27.3 day period |
|
||||||
|
|
||||||
**Stability Criteria:** |
|
||||||
- Energy drift < 1% over 100 orbits |
|
||||||
- Distance drift < 5% |
|
||||||
- No SOI transitions (body doesn't change parent) |
|
||||||
|
|
||||||
**Configuration:** Uses `test_time_step_stability.toml` |
|
||||||
|
|
||||||
**Usage:** |
|
||||||
```bash |
|
||||||
# From project root directory: |
|
||||||
./tests/informational/test_time_step_stability |
|
||||||
|
|
||||||
# Or using the Makefile: |
|
||||||
cd tests/informational |
|
||||||
make run-timestep |
|
||||||
|
|
||||||
# Run specific test case |
|
||||||
./tests/informational/test_time_step_stability '[timestep][stability]' |
|
||||||
|
|
||||||
# Run with verbose output to see binary search progress |
|
||||||
./tests/informational/test_time_step_stability -s "[timestep][stability]" |
|
||||||
``` |
|
||||||
|
|
||||||
**Note:** Tests must be run from the project root directory because config file paths are relative to the root. |
|
||||||
|
|
||||||
**Expected Output:** |
|
||||||
The test performs a binary search to find the maximum stable time step, printing progress for each dt value tested. Results show the minimum stable dt across all tested bodies with recommendations. |
|
||||||
|
|
||||||
## Adding New Informational Tests |
|
||||||
|
|
||||||
1. Create a new `.cpp` file in this directory |
|
||||||
2. Add corresponding `.toml` config file if needed |
|
||||||
3. Add a target in the `Makefile`: |
|
||||||
```makefile |
|
||||||
your_test_name: your_test_name.o |
|
||||||
g++ $(BUILD_DIR)/*.o -o $@ -lCatch2Main -lCatch2 -lm |
|
||||||
``` |
|
||||||
4. Update this README with test description |
|
||||||
5. Follow the naming convention: `test_<purpose>.cpp` and `test_<purpose>.toml` |
|
||||||
@ -1,296 +0,0 @@ |
|||||||
#include <catch2/catch_test_macros.hpp> |
|
||||||
#include "../../src/physics.h" |
|
||||||
#include "../../src/simulation.h" |
|
||||||
#include "../../src/config_loader.h" |
|
||||||
#include "../../src/test_utilities.h" |
|
||||||
#include <cmath> |
|
||||||
#include <cstdio> |
|
||||||
|
|
||||||
struct TestBody { |
|
||||||
const char* name; |
|
||||||
int body_index; |
|
||||||
int parent_index; |
|
||||||
double expected_period_days; |
|
||||||
}; |
|
||||||
|
|
||||||
struct StabilityResult { |
|
||||||
const char* name; |
|
||||||
double max_stable_dt; |
|
||||||
double period_days; |
|
||||||
}; |
|
||||||
|
|
||||||
const double SECONDS_PER_DAY = 86400.0; |
|
||||||
const double MIN_DT = 30.0; |
|
||||||
const double MAX_DT = 600.0; |
|
||||||
const double ENERGY_TOLERANCE = 1.0; |
|
||||||
const int NUM_ORBITS = 100; |
|
||||||
|
|
||||||
double calculate_orbital_period(CelestialBody* body, CelestialBody* parent) { |
|
||||||
Vec3 relative_pos = vec3_sub(body->global_position, parent->global_position); |
|
||||||
double r = vec3_magnitude(relative_pos); |
|
||||||
Vec3 relative_vel = vec3_sub(body->global_velocity, parent->global_velocity); |
|
||||||
double v = vec3_magnitude(relative_vel); |
|
||||||
|
|
||||||
double specific_energy = (v * v) / 2.0 - G * parent->mass / r; |
|
||||||
double semi_major_axis = -G * parent->mass / (2.0 * specific_energy); |
|
||||||
double period_seconds = 2.0 * M_PI * sqrt(pow(semi_major_axis, 3.0) / (G * parent->mass)); |
|
||||||
|
|
||||||
return period_seconds; |
|
||||||
} |
|
||||||
|
|
||||||
bool is_dt_stable(SimulationState* sim, const TestBody& test_body, double dt, int num_orbits) { |
|
||||||
SimulationState* test_sim = create_simulation(sim->max_bodies, sim->max_craft, sim->max_maneuvers, dt); |
|
||||||
|
|
||||||
REQUIRE(load_system_config(test_sim, "tests/informational/test_time_step_stability.toml")); |
|
||||||
|
|
||||||
int body_index = test_body.body_index; |
|
||||||
int parent_index = test_body.parent_index; |
|
||||||
|
|
||||||
double initial_energy = calculate_system_total_energy(test_sim); |
|
||||||
|
|
||||||
Vec3 initial_pos_relative = vec3_sub( |
|
||||||
test_sim->bodies[body_index].global_position, |
|
||||||
test_sim->bodies[parent_index].global_position |
|
||||||
); |
|
||||||
double initial_distance = vec3_magnitude(initial_pos_relative); |
|
||||||
|
|
||||||
double period = calculate_orbital_period(&test_sim->bodies[body_index], &test_sim->bodies[parent_index]); |
|
||||||
double max_time = period * num_orbits; |
|
||||||
|
|
||||||
bool completed = true; |
|
||||||
while (test_sim->time < max_time) { |
|
||||||
update_simulation(test_sim); |
|
||||||
|
|
||||||
if (test_sim->bodies[body_index].parent_index != parent_index) { |
|
||||||
completed = false; |
|
||||||
break; |
|
||||||
} |
|
||||||
} |
|
||||||
|
|
||||||
double final_energy = calculate_system_total_energy(test_sim); |
|
||||||
double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0; |
|
||||||
|
|
||||||
Vec3 final_pos_relative = vec3_sub( |
|
||||||
test_sim->bodies[body_index].global_position, |
|
||||||
test_sim->bodies[parent_index].global_position |
|
||||||
); |
|
||||||
double final_distance = vec3_magnitude(final_pos_relative); |
|
||||||
double distance_drift_percent = fabs((final_distance - initial_distance) / initial_distance) * 100.0; |
|
||||||
|
|
||||||
bool stable = completed && (energy_drift_percent < ENERGY_TOLERANCE) && (distance_drift_percent < 5.0); |
|
||||||
|
|
||||||
destroy_simulation(test_sim); |
|
||||||
|
|
||||||
return stable; |
|
||||||
} |
|
||||||
|
|
||||||
double find_max_stable_dt(SimulationState* sim, const TestBody& test_body) { |
|
||||||
double low = MIN_DT; |
|
||||||
double high = MAX_DT; |
|
||||||
double max_stable = low; |
|
||||||
|
|
||||||
printf("Testing %s (period ~%.2f days):\n", test_body.name, test_body.expected_period_days); |
|
||||||
|
|
||||||
for (int iter = 0; iter < 10; iter++) { |
|
||||||
double mid = (low + high) / 2.0; |
|
||||||
bool stable = is_dt_stable(sim, test_body, mid, NUM_ORBITS); |
|
||||||
|
|
||||||
if (stable) { |
|
||||||
max_stable = mid; |
|
||||||
low = mid; |
|
||||||
printf(" dt=%.0fs: STABLE\n", mid); |
|
||||||
} else { |
|
||||||
high = mid; |
|
||||||
printf(" dt=%.0fs: UNSTABLE\n", mid); |
|
||||||
} |
|
||||||
|
|
||||||
if (high - low < 5.0) break; |
|
||||||
} |
|
||||||
|
|
||||||
printf(" Maximum stable dt: %.0f seconds\n\n", max_stable); |
|
||||||
return max_stable; |
|
||||||
} |
|
||||||
|
|
||||||
void print_summary(const StabilityResult* results, int num_results, double min_stable_dt, double default_dt) { |
|
||||||
printf("\n"); |
|
||||||
printf("===============================================================================\n"); |
|
||||||
printf(" TIME STEP STABILITY TEST RESULTS\n"); |
|
||||||
printf("===============================================================================\n\n"); |
|
||||||
|
|
||||||
printf("STABILITY CRITERIA:\n"); |
|
||||||
printf(" - Energy drift < %.1f%% over %d orbits\n", ENERGY_TOLERANCE, NUM_ORBITS); |
|
||||||
printf(" - Distance drift < 5.0%%\n"); |
|
||||||
printf(" - No SOI transitions (parent changes)\n\n"); |
|
||||||
|
|
||||||
printf("PER-BODY RESULTS:\n"); |
|
||||||
printf("+----------------------+----------------+------------------+----------------+\n"); |
|
||||||
printf("| Body | Period (days) | Max Stable dt (s) | Stability Status |\n"); |
|
||||||
printf("+----------------------+----------------+------------------+----------------+\n"); |
|
||||||
|
|
||||||
for (int i = 0; i < num_results; i++) { |
|
||||||
const StabilityResult& r = results[i]; |
|
||||||
double ratio = default_dt / r.max_stable_dt; |
|
||||||
const char* status = ratio < 0.5 ? "Very Stable" : ratio < 0.8 ? "Stable" : "Limited Margin"; |
|
||||||
printf("| %-20s | %14.2f | %16.0f | %-14s |\n", r.name, r.period_days, r.max_stable_dt, status); |
|
||||||
} |
|
||||||
|
|
||||||
printf("+----------------------+----------------+------------------+----------------+\n\n"); |
|
||||||
|
|
||||||
printf("OVERALL ANALYSIS:\n"); |
|
||||||
printf(" Minimum stable time step: %.0f seconds\n", min_stable_dt); |
|
||||||
printf(" Recommended safe dt: %.0f seconds (0.7x safety margin)\n", min_stable_dt * 0.7); |
|
||||||
printf(" Current default dt: %.0f seconds\n", default_dt); |
|
||||||
printf(" Current dt stability: %.0fx\n", default_dt / min_stable_dt); |
|
||||||
printf("\n"); |
|
||||||
|
|
||||||
if (default_dt < min_stable_dt * 0.7) { |
|
||||||
printf("STATUS: Current time step (60s) is VERY STABLE with good margin.\n"); |
|
||||||
printf(" Can be increased significantly if needed.\n\n"); |
|
||||||
} else if (default_dt < min_stable_dt) { |
|
||||||
printf("STATUS: Current time step (60s) is STABLE with adequate margin.\n"); |
|
||||||
printf(" Moderate increases possible.\n\n"); |
|
||||||
} else { |
|
||||||
printf("STATUS: Current time step (60s) is near stability limit.\n"); |
|
||||||
printf(" Consider reducing for safety.\n\n"); |
|
||||||
} |
|
||||||
|
|
||||||
printf("RECOMMENDATIONS:\n"); |
|
||||||
printf(" - For MESSENGER-like close orbits: Keep dt <= %.0f seconds\n", min_stable_dt); |
|
||||||
printf(" - For planetary missions: Current dt=60s is excellent\n"); |
|
||||||
printf(" - For Moon-scale orbits: Could use dt=120s+ safely\n"); |
|
||||||
printf("\n"); |
|
||||||
printf("===============================================================================\n\n"); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Time step stability - Mercury orbiter (MESSENGER-like)", "[timestep][stability]") { |
|
||||||
const double BASE_DT = 60.0; |
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, BASE_DT); |
|
||||||
|
|
||||||
TestBody mercury_orbiter = {"Mercury_Orbiter", 1, 0, 0.5}; |
|
||||||
double max_dt = find_max_stable_dt(sim, mercury_orbiter); |
|
||||||
|
|
||||||
INFO("Mercury orbiter maximum stable dt: " << max_dt << " seconds"); |
|
||||||
|
|
||||||
REQUIRE(max_dt >= MIN_DT); |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Time step stability - Io (Jupiter's moon)", "[timestep][stability]") { |
|
||||||
const double BASE_DT = 60.0; |
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, BASE_DT); |
|
||||||
|
|
||||||
TestBody io = {"Io", 3, 2, 1.77}; |
|
||||||
double max_dt = find_max_stable_dt(sim, io); |
|
||||||
|
|
||||||
INFO("Io maximum stable dt: " << max_dt << " seconds"); |
|
||||||
|
|
||||||
REQUIRE(max_dt >= MIN_DT); |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Time step stability - Moon (Earth's moon)", "[timestep][stability]") { |
|
||||||
const double BASE_DT = 60.0; |
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, BASE_DT); |
|
||||||
|
|
||||||
TestBody moon = {"Moon", 5, 4, 27.3}; |
|
||||||
double max_dt = find_max_stable_dt(sim, moon); |
|
||||||
|
|
||||||
INFO("Moon maximum stable dt: " << max_dt << " seconds"); |
|
||||||
|
|
||||||
REQUIRE(max_dt >= MIN_DT); |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Find minimum stable time step across all bodies", "[timestep][stability]") { |
|
||||||
const double BASE_DT = 60.0; |
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, BASE_DT); |
|
||||||
|
|
||||||
TestBody bodies[] = { |
|
||||||
{"Mercury_Orbiter", 1, 0, 0.5}, |
|
||||||
{"Io", 3, 2, 1.77}, |
|
||||||
{"Moon", 5, 4, 27.3} |
|
||||||
}; |
|
||||||
|
|
||||||
StabilityResult results[3]; |
|
||||||
double max_dt = MAX_DT; |
|
||||||
|
|
||||||
printf("\n=== Finding minimum stable dt across all bodies ===\n\n"); |
|
||||||
|
|
||||||
for (int i = 0; i < 3; i++) { |
|
||||||
results[i].name = bodies[i].name; |
|
||||||
results[i].period_days = bodies[i].expected_period_days; |
|
||||||
results[i].max_stable_dt = find_max_stable_dt(sim, bodies[i]); |
|
||||||
if (results[i].max_stable_dt < max_dt) { |
|
||||||
max_dt = results[i].max_stable_dt; |
|
||||||
} |
|
||||||
} |
|
||||||
|
|
||||||
printf("\n=== RESULTS ===\n"); |
|
||||||
printf("Minimum stable time step: %.0f seconds\n", max_dt); |
|
||||||
printf("Recommended safe time step: %.0f seconds (%.0fx safety margin)\n", max_dt * 0.7, 1.0/0.7); |
|
||||||
|
|
||||||
INFO("Minimum stable dt: " << max_dt << " seconds"); |
|
||||||
|
|
||||||
print_summary(results, 3, max_dt, BASE_DT); |
|
||||||
|
|
||||||
REQUIRE(max_dt >= MIN_DT); |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
|
|
||||||
TEST_CASE("Verify current default dt (60s) stability", "[timestep][stability]") { |
|
||||||
const double DT = 60.0; |
|
||||||
const int NUM_ORBITS = 10; |
|
||||||
|
|
||||||
SimulationState* sim = create_simulation(10, 0, 0, DT); |
|
||||||
REQUIRE(load_system_config(sim, "tests/informational/test_time_step_stability.toml")); |
|
||||||
|
|
||||||
struct BodyTest { |
|
||||||
int body_index; |
|
||||||
int parent_index; |
|
||||||
const char* name; |
|
||||||
}; |
|
||||||
|
|
||||||
BodyTest tests[] = { |
|
||||||
{1, 0, "Mercury_Orbiter"}, |
|
||||||
{3, 2, "Io"}, |
|
||||||
{5, 4, "Moon"} |
|
||||||
}; |
|
||||||
|
|
||||||
for (int t = 0; t < 3; t++) { |
|
||||||
int body_index = tests[t].body_index; |
|
||||||
int parent_index = tests[t].parent_index; |
|
||||||
const char* name = tests[t].name; |
|
||||||
|
|
||||||
double period = calculate_orbital_period(&sim->bodies[body_index], &sim->bodies[parent_index]); |
|
||||||
double max_time = period * NUM_ORBITS; |
|
||||||
|
|
||||||
double initial_energy = calculate_system_total_energy(sim); |
|
||||||
|
|
||||||
INFO("Testing " << name << " with dt=" << DT << "s for " << NUM_ORBITS << " orbits"); |
|
||||||
|
|
||||||
bool completed = true; |
|
||||||
while (sim->time < max_time) { |
|
||||||
update_simulation(sim); |
|
||||||
|
|
||||||
if (sim->bodies[body_index].parent_index != parent_index) { |
|
||||||
completed = false; |
|
||||||
break; |
|
||||||
} |
|
||||||
} |
|
||||||
|
|
||||||
double final_energy = calculate_system_total_energy(sim); |
|
||||||
double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0; |
|
||||||
|
|
||||||
INFO(name << " completed: " << (completed ? "yes" : "no")); |
|
||||||
INFO(name << " energy drift: " << energy_drift_percent << "%"); |
|
||||||
|
|
||||||
REQUIRE(completed); |
|
||||||
REQUIRE(energy_drift_percent < ENERGY_TOLERANCE); |
|
||||||
} |
|
||||||
|
|
||||||
destroy_simulation(sim); |
|
||||||
} |
|
||||||
@ -1,80 +0,0 @@ |
|||||||
# Time Step Stability Test Configuration |
|
||||||
# Contains bodies with different orbital periods for testing RK4 stability |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Mercury" |
|
||||||
mass = 3.285e23 |
|
||||||
radius = 2.439e6 |
|
||||||
parent_index = -1 |
|
||||||
color = { r = 0.7, g = 0.7, b = 0.7 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 0.0, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
# Low Mercury orbiter at 200 km altitude (MESSENGER-like) |
|
||||||
# Period ~12 hours - most restrictive case |
|
||||||
[[bodies]] |
|
||||||
name = "Mercury_Orbiter" |
|
||||||
mass = 1000.0 |
|
||||||
radius = 1.0 |
|
||||||
parent_index = 0 |
|
||||||
color = { r = 1.0, g = 0.0, b = 1.0 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 2.639e6, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
[[bodies]] |
|
||||||
name = "Jupiter" |
|
||||||
mass = 1.898e27 |
|
||||||
radius = 6.9911e7 |
|
||||||
parent_index = -1 |
|
||||||
color = { r = 0.9, g = 0.7, b = 0.5 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 0.0, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
# Io - Jupiter's moon |
|
||||||
# Period ~1.77 days |
|
||||||
[[bodies]] |
|
||||||
name = "Io" |
|
||||||
mass = 8.93e22 |
|
||||||
radius = 1.822e6 |
|
||||||
parent_index = 2 |
|
||||||
color = { r = 0.9, g = 0.9, b = 0.3 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 4.217e8, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
|
|
||||||
[[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 |
|
||||||
} |
|
||||||
|
|
||||||
# Moon - Earth's moon |
|
||||||
# Period ~27.3 days |
|
||||||
[[bodies]] |
|
||||||
name = "Moon" |
|
||||||
mass = 7.342e22 |
|
||||||
radius = 1.737e6 |
|
||||||
parent_index = 4 |
|
||||||
color = { r = 0.7, g = 0.7, b = 0.7 } |
|
||||||
orbit = { |
|
||||||
semi_major_axis = 3.844e8, |
|
||||||
eccentricity = 0.0, |
|
||||||
true_anomaly = 0.0 |
|
||||||
} |
|
||||||
Loading…
Reference in new issue