diff --git a/tests/test_energy.cpp b/tests/test_energy.cpp new file mode 100644 index 0000000..caa5830 --- /dev/null +++ b/tests/test_energy.cpp @@ -0,0 +1,60 @@ +#include +#include +#include "../src/physics.h" +#include "../src/simulation.h" +#include "../src/config_loader.h" +#include "../src/test_utilities.h" +#include + +using Catch::Matchers::WithinAbs; + +SCENARIO("Energy conservation in circular orbit", "[energy][sanity]") { + const double TIME_STEP = 60.0; + const double DAYS_TO_SIMULATE = 10.0; + const double SECONDS_PER_DAY = 86400.0; + + SimulationState* sim = create_simulation(2, 0, 0, TIME_STEP); + REQUIRE(load_system_config(sim, "tests/test_energy.toml")); + + const double initial_energy = calculate_system_total_energy(sim); + + double total_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; + while (sim->time < total_time) { + update_simulation(sim); + } + + const double final_energy = calculate_system_total_energy(sim); + const double energy_drift = fabs(final_energy - initial_energy) / fabs(initial_energy); + + INFO("Initial energy: " << initial_energy << " J"); + INFO("Final energy: " << final_energy << " J"); + INFO("Relative drift: " << energy_drift << " (fraction)"); + + REQUIRE_THAT(energy_drift, WithinAbs(0.0, 1e-12)); + + destroy_simulation(sim); +} + +SCENARIO("Orbit direction for zero inclination", "[direction][sanity]") { + const double TIME_STEP = 60.0; + const int STEPS = 1440; // 1 day at 60s steps + + 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 = vec3_sub(earth->global_position, sun->global_position); + const double theta_start = atan2(initial_rel.y, initial_rel.x); + + for (int i = 0; i < STEPS; i++) update_simulation(sim); + + Vec3 final_rel = vec3_sub(earth->global_position, sun->global_position); + const double delta = atan2(final_rel.y, final_rel.x) - theta_start; + + INFO("Delta: " << delta << " rad"); + REQUIRE_THAT(delta, WithinAbs(0.0172042841, 1e-6)); + + destroy_simulation(sim); +} diff --git a/tests/test_orbital_period.cpp b/tests/test_orbital_period.cpp index c399a88..ecc6ab9 100644 --- a/tests/test_orbital_period.cpp +++ b/tests/test_orbital_period.cpp @@ -16,80 +16,39 @@ SCENARIO("Orbital period measurement with analytical propagation", SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP); REQUIRE(load_system_config(sim, "tests/test_orbital_period.toml")); - SECTION("Earth completes one orbit in ~365 days") { - OrbitTracker* tracker = create_orbit_tracker(1); - CelestialBody* earth = &sim->bodies[1]; - CelestialBody* sun = &sim->bodies[0]; - - const double MAX_DAYS = 400.0; - const double max_time = MAX_DAYS * SECONDS_PER_DAY; - + auto propagate_and_check = [&](OrbitTracker* tracker, CelestialBody* body, + CelestialBody* parent, double max_time, + double expected_days) { while (sim->time < max_time && !tracker->orbit_completed) { update_simulation(sim); - update_orbit_tracker(tracker, earth, sun, sim->time); + update_orbit_tracker(tracker, body, parent, sim->time); } - REQUIRE(tracker->orbit_completed); - double measured_days = tracker->time_at_completion / SECONDS_PER_DAY; + const double measured_days = tracker->time_at_completion / SECONDS_PER_DAY; INFO("Measured period: " << measured_days << " days"); - - REQUIRE_THAT(measured_days, WithinAbs(365.2105, 0.1)); + REQUIRE_THAT(measured_days, WithinAbs(expected_days, 0.1)); destroy_orbit_tracker(tracker); + }; + + SECTION("Earth completes one orbit in ~365 days") { + CelestialBody* earth = &sim->bodies[1]; + CelestialBody* sun = &sim->bodies[0]; + OrbitTracker* tracker = create_orbit_tracker(1); + + const double max_time = 400.0 * SECONDS_PER_DAY; + propagate_and_check(tracker, earth, sun, max_time, 365.2105); } SECTION("Mars completes one orbit in ~671 days") { - OrbitTracker* tracker = create_orbit_tracker(2); CelestialBody* mars = &sim->bodies[2]; CelestialBody* sun = &sim->bodies[0]; + OrbitTracker* tracker = create_orbit_tracker(2); - const double MAX_DAYS = 750.0; - const double max_time = MAX_DAYS * SECONDS_PER_DAY; - - while (sim->time < max_time && !tracker->orbit_completed) { - update_simulation(sim); - update_orbit_tracker(tracker, mars, sun, sim->time); - } - - REQUIRE(tracker->orbit_completed); - - double measured_days = tracker->time_at_completion / SECONDS_PER_DAY; - INFO("Measured period: " << measured_days << " days"); - - REQUIRE_THAT(measured_days, WithinAbs(670.9345, 0.1)); - - destroy_orbit_tracker(tracker); + const double max_time = 750.0 * SECONDS_PER_DAY; + propagate_and_check(tracker, mars, sun, max_time, 670.9345); } } -SCENARIO("Orbit direction for zero inclination", - "[direction][sanity]") { - const double TIME_STEP = 60.0; - const double TEST_DURATION_DAYS = 1.0; - const int STEPS = (int)(TEST_DURATION_DAYS * 86400.0 / 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 = vec3_sub(earth->global_position, sun->global_position); - double theta_start = atan2(initial_rel.y, initial_rel.x); - - for (int i = 0; i < STEPS; i++) { - update_simulation(sim); - } - - Vec3 final_rel = vec3_sub(earth->global_position, sun->global_position); - double theta_final = atan2(final_rel.y, final_rel.x); - double delta = theta_final - theta_start; - - INFO("Initial angle: " << theta_start << " rad"); - INFO("Final angle: " << theta_final << " rad"); - INFO("Delta: " << delta << " rad"); - - REQUIRE_THAT(delta, WithinAbs(0.0172042841, 1e-6)); - REQUIRE(delta > 0.0); -}