Refactor orbital elements to reusable code and track SOI transitions

- 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
6 months ago · 0e2f9dd92f
4 changed files with 225 additions and 1 deletions
--- a/3
+++ b/3
@ -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 \
--- 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;
+}
--- 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
--- a/tests/test_comet_orbit.cpp
+++ b/tests/test_comet_orbit.cpp
@ -0,0 +1,163 @@
+#include <catch2/catch_test_macros.hpp>
+#include "../src/physics.h"
+#include "../src/bodies.h"
+#include "../src/test_utilities.h"
+#include <cmath>
+#include <vector>
+
+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<OrbitalElements> snapshots;
+    std::vector<ParentChange> 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);
+}