From 0e2f9dd92f66704a87c730bde9ce7e3e9fb80267 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sun, 4 Jan 2026 18:06:21 -0500 Subject: [PATCH] Refactor orbital elements to reusable code and track SOI transitions MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Move OrbitalElements struct to src/bodies.h - Move calculate_orbital_elements() to src/bodies.cpp - Add test_comet_orbit.cpp with parent index change tracking - Update Makefile to include new test file 🤖 Refactored orbital mechanics code for better reusability --- Makefile | 3 +- src/bodies.cpp | 46 +++++++++++ src/bodies.h | 14 ++++ tests/test_comet_orbit.cpp | 163 +++++++++++++++++++++++++++++++++++++ 4 files changed, 225 insertions(+), 1 deletion(-) create mode 100644 tests/test_comet_orbit.cpp diff --git a/Makefile b/Makefile index 015d871..e51e2b1 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ # Compiler and flags CXX = g++ -CXXFLAGS = -Wall -Wextra -std=c++14 -I./src -isystem./ext/raylib/src +CXXFLAGS = -Wall -Wextra -g -gddb3 -std=c++14 -I./src -isystem./ext/raylib/src LDFLAGS = -L./ext/raylib/src -lraylib -lm -lpthread -ldl -lrt -lX11 # Directories @@ -69,6 +69,7 @@ test-build: $(TEST_DIR)/test_integration.cpp \ $(TEST_DIR)/test_energy.cpp \ $(TEST_DIR)/test_orbital_period.cpp \ + $(TEST_DIR)/test_comet_orbit.cpp \ $(SRC_DIR)/test_utilities.cpp \ $(SRC_DIR)/physics.cpp \ $(SRC_DIR)/bodies.cpp \ diff --git a/src/bodies.cpp b/src/bodies.cpp index 8250c37..f691d07 100644 --- a/src/bodies.cpp +++ b/src/bodies.cpp @@ -137,3 +137,49 @@ void update_simulation(SimulationState* sim) { sim->time += sim->dt; } + +OrbitalElements calculate_orbital_elements(CelestialBody* body, CelestialBody* primary, + CelestialBody* optional_ref_body, double current_time) { + const double AU = 1.496e11; + const double SECONDS_PER_DAY = 86400.0; + const double M_sun = primary->mass; + + OrbitalElements elem; + elem.time_days = current_time / SECONDS_PER_DAY; + + Vec3 r_vec = vec3_sub(body->position, primary->position); + double r = vec3_magnitude(r_vec); + double v = vec3_magnitude(body->velocity); + + elem.distance_to_sun_au = r / AU; + elem.velocity_magnitude = v; + + if (optional_ref_body) { + double dist_ref = vec3_distance(body->position, optional_ref_body->position); + elem.distance_to_ref_body_au = dist_ref / AU; + } else { + elem.distance_to_ref_body_au = -1.0; + } + + elem.specific_energy = (v * v) / 2.0 - (G * M_sun) / r; + + if (elem.specific_energy < 0) { + elem.semi_major_axis_au = -(G * M_sun) / (2.0 * elem.specific_energy) / AU; + + double v_squared = v * v; + double r_dot_v = r_vec.x * body->velocity.x + r_vec.y * body->velocity.y + r_vec.z * body->velocity.z; + + Vec3 e_vec; + e_vec.x = (v_squared - G * M_sun / r) * r_vec.x - r_dot_v * body->velocity.x; + e_vec.y = (v_squared - G * M_sun / r) * r_vec.y - r_dot_v * body->velocity.y; + e_vec.z = (v_squared - G * M_sun / r) * r_vec.z - r_dot_v * body->velocity.z; + + double e_mag = vec3_magnitude(e_vec) / (G * M_sun); + elem.eccentricity = e_mag; + } else { + elem.semi_major_axis_au = 0.0; + elem.eccentricity = 1.0; + } + + return elem; +} diff --git a/src/bodies.h b/src/bodies.h index 648ad1a..d97e83a 100644 --- a/src/bodies.h +++ b/src/bodies.h @@ -38,4 +38,18 @@ int find_dominant_body(SimulationState* sim, int body_index); void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_axis); void update_simulation(SimulationState* sim); +// Orbital elements calculation +struct OrbitalElements { + double time_days; + double semi_major_axis_au; + double eccentricity; + double specific_energy; + double distance_to_sun_au; + double distance_to_ref_body_au; + double velocity_magnitude; +}; + +OrbitalElements calculate_orbital_elements(CelestialBody* body, CelestialBody* primary, + CelestialBody* optional_ref_body, double current_time); + #endif diff --git a/tests/test_comet_orbit.cpp b/tests/test_comet_orbit.cpp new file mode 100644 index 0000000..90510de --- /dev/null +++ b/tests/test_comet_orbit.cpp @@ -0,0 +1,163 @@ +#include +#include "../src/physics.h" +#include "../src/bodies.h" +#include "../src/test_utilities.h" +#include +#include + +struct ParentChange { + double time_seconds; + double time_days; + int old_parent; + int new_parent; + double distance_to_mars_au; + double distance_to_sun_au; +}; + +TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[comet][orbital-elements][soi]") { + const double TIME_STEP = 60.0; + const double SECONDS_PER_DAY = 86400.0; + const double DAYS_TO_SIMULATE = 1700.0; + const double AU = 1.496e11; + + const double EXPECTED_SMA = 2.5; + const double EXPECTED_ECC = 0.7; + + SimulationState* sim = create_simulation(10, TIME_STEP); + + Vec3 sun_pos = {0, 0, 0}; + Vec3 sun_vel = {0, 0, 0}; + add_body(sim, "Sun", 1.989e30, 6.96e8, sun_pos, sun_vel, -1, 1.0, 1.0, 0.0, 0, 0); + + Vec3 earth_pos = {1.496e11, 0, 0}; + Vec3 earth_vel = {0, 29789, 0}; + add_body(sim, "Earth", 5.972e24, 6.371e6, earth_pos, earth_vel, 0, 0.0, 0.5, 1.0, 0, 1.496e11); + + Vec3 mars_pos = {2.244e11, 0, 0}; + Vec3 mars_vel = {0, 24323, 0}; + add_body(sim, "Mars", 6.39e23, 3.3895e6, mars_pos, mars_vel, 0, 0.8, 0.3, 0.1, 0, 2.244e11); + + Vec3 comet_pos = {1.122e11, 0, 0}; + Vec3 comet_vel = {0, 44849, 0}; + add_body(sim, "Comet", 1e14, 5e3, comet_pos, comet_vel, 0, 0.5, 0.8, 1.0, 0.7, 3.74e11); + + const int COMET_INDEX = 3; + const int MARS_INDEX = 2; + const int SUN_INDEX = 0; + + update_soi(&sim->bodies[SUN_INDEX], NULL, 0); + update_soi(&sim->bodies[MARS_INDEX], &sim->bodies[SUN_INDEX], 2.244e11); + update_soi(&sim->bodies[COMET_INDEX], &sim->bodies[SUN_INDEX], 3.74e11); + + std::vector snapshots; + std::vector parent_changes; + int previous_parent = sim->bodies[COMET_INDEX].parent_index; + + OrbitalElements initial = calculate_orbital_elements( + &sim->bodies[COMET_INDEX], + &sim->bodies[SUN_INDEX], + &sim->bodies[MARS_INDEX], + sim->time + ); + snapshots.push_back(initial); + + double checkpoint_days[] = {0, 365, 722, 1444, 1533, 1600, 1700}; + size_t next_checkpoint = 1; + + double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; + while (sim->time < max_time) { + update_simulation(sim); + + int current_parent = sim->bodies[COMET_INDEX].parent_index; + if (current_parent != previous_parent) { + ParentChange change; + change.time_seconds = sim->time; + change.time_days = sim->time / SECONDS_PER_DAY; + change.old_parent = previous_parent; + change.new_parent = current_parent; + change.distance_to_sun_au = vec3_magnitude(sim->bodies[COMET_INDEX].position) / AU; + change.distance_to_mars_au = vec3_distance(sim->bodies[COMET_INDEX].position, + sim->bodies[MARS_INDEX].position) / AU; + parent_changes.push_back(change); + previous_parent = current_parent; + } + + double current_days = sim->time / SECONDS_PER_DAY; + if (next_checkpoint < sizeof(checkpoint_days)/sizeof(checkpoint_days[0]) && + current_days >= checkpoint_days[next_checkpoint]) { + + OrbitalElements elem = calculate_orbital_elements( + &sim->bodies[COMET_INDEX], + &sim->bodies[SUN_INDEX], + &sim->bodies[MARS_INDEX], + sim->time + ); + snapshots.push_back(elem); + next_checkpoint++; + } + } + + printf("\n=== Comet Orbital Elements Over Time ===\n"); + printf("Expected: a = %.1f AU, e = %.1f\n\n", EXPECTED_SMA, EXPECTED_ECC); + + for (const auto& elem : snapshots) { + printf("Day %.0f:\n", elem.time_days); + printf(" Semi-major axis: %.6f AU (expected: %.1f AU)\n", elem.semi_major_axis_au, EXPECTED_SMA); + printf(" Eccentricity: %.6f (expected: %.1f)\n", elem.eccentricity, EXPECTED_ECC); + printf(" Distance to Sun: %.4f AU\n", elem.distance_to_sun_au); + printf(" Distance to Mars: %.6f AU\n", elem.distance_to_ref_body_au); + printf(" Velocity: %.2f km/s\n", elem.velocity_magnitude / 1000.0); + printf(" Specific energy: %.3e J/kg\n", elem.specific_energy); + + double sma_error = fabs(elem.semi_major_axis_au - EXPECTED_SMA); + double ecc_error = fabs(elem.eccentricity - EXPECTED_ECC); + + if (elem.time_days >= 1530 && elem.time_days <= 1540) { + printf(" *** CLOSEST APPROACH TO MARS ***\n"); + } + + printf(" Error: Δa = %.6f AU, Δe = %.6f\n\n", sma_error, ecc_error); + } + + printf("\n=== Parent Index Changes ===\n"); + if (parent_changes.empty()) { + printf("No parent changes detected\n\n"); + } else { + for (const auto& change : parent_changes) { + printf("Time: %.0f days (%.0f seconds)\n", change.time_days, change.time_seconds); + printf(" Parent: %d -> %d", change.old_parent, change.new_parent); + if (change.new_parent == MARS_INDEX) { + printf(" (Sun -> Mars)"); + } else if (change.old_parent == MARS_INDEX) { + printf(" (Mars -> Sun)"); + } + printf("\n"); + printf(" Distance to Mars: %.6f AU\n", change.distance_to_mars_au); + printf(" Distance to Sun: %.6f AU\n\n", change.distance_to_sun_au); + } + } + + OrbitalElements final = snapshots.back(); + double final_sma_error = fabs(final.semi_major_axis_au - EXPECTED_SMA); + double final_ecc_error = fabs(final.eccentricity - EXPECTED_ECC); + + INFO("Final drift from initial:"); + INFO(" SMA: " << final_sma_error << " AU (" << (final_sma_error / EXPECTED_SMA * 100.0) << "%)"); + INFO(" ECC: " << final_ecc_error << " (" << (final_ecc_error / EXPECTED_ECC * 100.0) << "%)"); + + REQUIRE(final_sma_error < 0.1); + REQUIRE(final_ecc_error < 0.05); + + REQUIRE(parent_changes.size() > 0); + + bool found_mars_transition = false; + for (const auto& change : parent_changes) { + if (change.new_parent == MARS_INDEX || change.old_parent == MARS_INDEX) { + found_mars_transition = true; + REQUIRE(fabs(change.time_days - 1550.0) < 50.0); + } + } + REQUIRE(found_mars_transition); + + destroy_simulation(sim); +}