From 89931c357166246c57c5531ddf8a2164f129594c Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Mon, 5 Jan 2026 12:25:24 -0500 Subject: [PATCH] Fix SOI calculation and orbit direction bugs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Fix SOI calculation to use semi_major_axis instead of current distance This caused incorrect SOI radii for eccentric orbits - Fix orbit direction: change cross product from r×z to z×r for counter-clockwise orbits - Add parent body change tracking to headless mode for debugging - Add debug output enabling in test 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.5 --- src/bodies.cpp | 1 + src/config_loader.cpp | 39 +++++++++++++++++++------------------- src/main.cpp | 18 ++++++++++++++++++ tests/test_comet_orbit.cpp | 13 ++++++++++--- 4 files changed, 49 insertions(+), 22 deletions(-) diff --git a/src/bodies.cpp b/src/bodies.cpp index f691d07..2446cc0 100644 --- a/src/bodies.cpp +++ b/src/bodies.cpp @@ -1,6 +1,7 @@ #include "bodies.h" #include #include +#include #include // Create a new simulation diff --git a/src/config_loader.cpp b/src/config_loader.cpp index f754d26..a16ee2c 100644 --- a/src/config_loader.cpp +++ b/src/config_loader.cpp @@ -132,21 +132,21 @@ void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params) { // v = sqrt(G * M_other / r) double speed = sqrt(G * other_mass / distance); - // Create velocity perpendicular to position vector (same logic as below) + // Create velocity perpendicular to position vector (counter-clockwise from +Z) Vec3 z_axis = {0.0, 0.0, 1.0}; Vec3 vel_dir = { - r.y * z_axis.z - r.z * z_axis.y, - r.z * z_axis.x - r.x * z_axis.z, - r.x * z_axis.y - r.y * z_axis.x + z_axis.y * r.z - z_axis.z * r.y, + z_axis.z * r.x - z_axis.x * r.z, + z_axis.x * r.y - z_axis.y * r.x }; // If r is parallel to z-axis, use x-axis instead double cross_mag = vec3_magnitude(vel_dir); if (cross_mag < 0.01) { Vec3 x_axis = {1.0, 0.0, 0.0}; - vel_dir.x = r.y * x_axis.z - r.z * x_axis.y; - vel_dir.y = r.z * x_axis.x - r.x * x_axis.z; - vel_dir.z = r.x * x_axis.y - r.y * x_axis.x; + vel_dir.x = x_axis.y * r.z - x_axis.z * r.y; + vel_dir.y = x_axis.z * r.x - x_axis.x * r.z; + vel_dir.z = x_axis.x * r.y - x_axis.y * r.x; } // Normalize and scale by orbital speed @@ -197,23 +197,23 @@ void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params) { // Create velocity perpendicular to position vector // If position is mostly in XY plane, make velocity in XY plane - // Cross product of r with z-axis gives perpendicular vector in XY plane + // Cross product z_axis x r gives counter-clockwise orbit when viewed from +Z Vec3 z_axis = {0.0, 0.0, 1.0}; - // Calculate cross product: r x z_axis + // Calculate cross product: z_axis x r Vec3 vel_dir = { - r.y * z_axis.z - r.z * z_axis.y, - r.z * z_axis.x - r.x * z_axis.z, - r.x * z_axis.y - r.y * z_axis.x + z_axis.y * r.z - z_axis.z * r.y, + z_axis.z * r.x - z_axis.x * r.z, + z_axis.x * r.y - z_axis.y * r.x }; // If r is parallel to z-axis, use x-axis instead double cross_mag = vec3_magnitude(vel_dir); if (cross_mag < 0.01) { Vec3 x_axis = {1.0, 0.0, 0.0}; - vel_dir.x = r.y * x_axis.z - r.z * x_axis.y; - vel_dir.y = r.z * x_axis.x - r.x * x_axis.z; - vel_dir.z = r.x * x_axis.y - r.y * x_axis.x; + vel_dir.x = x_axis.y * r.z - x_axis.z * r.y; + vel_dir.y = x_axis.z * r.x - x_axis.x * r.z; + vel_dir.z = x_axis.x * r.y - x_axis.y * r.x; } // Normalize and scale by orbital speed @@ -228,6 +228,7 @@ void calculate_velocities(SimulationState* sim, OrbitParams* orbit_params) { // Calculate SOI radii for all bodies void calculate_soi_radii(SimulationState* sim) { + const double AU = 1.496e11; for (int i = 0; i < sim->body_count; i++) { CelestialBody* body = &sim->bodies[i]; @@ -237,11 +238,11 @@ void calculate_soi_radii(SimulationState* sim) { } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) { CelestialBody* parent = &sim->bodies[body->parent_index]; - // Calculate semi-major axis (distance to parent) - double semi_major_axis = vec3_distance(body->position, parent->position); - // Update SOI using Hill sphere approximation - update_soi(body, parent, semi_major_axis); + update_soi(body, parent, body->semi_major_axis); + + printf(" %s SOI radius: %.6e m (%.6f AU)\n", + body->name, body->soi_radius, body->soi_radius / AU); } } } diff --git a/src/main.cpp b/src/main.cpp index 2e3662c..5ffbf9e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -91,8 +91,10 @@ void run_headless_simulation(SimulationState* sim, const ProgramArgs* args) { double next_output_time = 0.0; CelestialBody* initial_state = (CelestialBody*)malloc(sim->body_count * sizeof(CelestialBody)); + int* previous_parents = (int*)malloc(sim->body_count * sizeof(int)); for (int i = 0; i < sim->body_count; i++) { initial_state[i] = sim->bodies[i]; + previous_parents[i] = sim->bodies[i].parent_index; } printf("\n=== Initial State ===\n"); @@ -113,6 +115,21 @@ void run_headless_simulation(SimulationState* sim, const ProgramArgs* args) { while (sim->time < total_time) { update_simulation(sim); + for (int i = 0; i < sim->body_count; i++) { + int current_parent = sim->bodies[i].parent_index; + if (current_parent != previous_parents[i]) { + printf("\n*** PARENT CHANGE at Day %.1f ***\n", sim->time / SECONDS_PER_DAY); + printf(" Body: %s\n", sim->bodies[i].name); + printf(" Old parent: %d", previous_parents[i]); + if (previous_parents[i] >= 0) printf(" (%s)", sim->bodies[previous_parents[i]].name); + printf("\n"); + printf(" New parent: %d", current_parent); + if (current_parent >= 0) printf(" (%s)", sim->bodies[current_parent].name); + printf("\n\n"); + previous_parents[i] = current_parent; + } + } + if (sim->time >= next_output_time) { printf("Day %.1f: ", sim->time / SECONDS_PER_DAY); for (int i = 0; i < sim->body_count && i < 3; i++) { @@ -161,6 +178,7 @@ void run_headless_simulation(SimulationState* sim, const ProgramArgs* args) { } free(initial_state); + free(previous_parents); } void run_gui_simulation(SimulationState* sim) { diff --git a/tests/test_comet_orbit.cpp b/tests/test_comet_orbit.cpp index 90510de..d59a303 100644 --- a/tests/test_comet_orbit.cpp +++ b/tests/test_comet_orbit.cpp @@ -17,7 +17,7 @@ struct ParentChange { 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 DAYS_TO_SIMULATE = 5000.0; const double AU = 1.496e11; const double EXPECTED_SMA = 2.5; @@ -61,7 +61,7 @@ TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[ ); snapshots.push_back(initial); - double checkpoint_days[] = {0, 365, 722, 1444, 1533, 1600, 1700}; + double checkpoint_days[] = {0, 365, 722, 1444, 1533, 1600, 1700, 2000, 5000}; size_t next_checkpoint = 1; double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; @@ -86,6 +86,11 @@ TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[ if (next_checkpoint < sizeof(checkpoint_days)/sizeof(checkpoint_days[0]) && current_days >= checkpoint_days[next_checkpoint]) { + printf("\n day: %f, sim->time: %f, comet parent index: %d", + current_days, + sim->time, + sim->bodies[COMET_INDEX].parent_index); + OrbitalElements elem = calculate_orbital_elements( &sim->bodies[COMET_INDEX], &sim->bodies[SUN_INDEX], @@ -97,6 +102,7 @@ TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[ } } +#if 1 printf("\n=== Comet Orbital Elements Over Time ===\n"); printf("Expected: a = %.1f AU, e = %.1f\n\n", EXPECTED_SMA, EXPECTED_ECC); @@ -120,6 +126,7 @@ TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[ } printf("\n=== Parent Index Changes ===\n"); + if (parent_changes.empty()) { printf("No parent changes detected\n\n"); } else { @@ -136,7 +143,7 @@ TEST_CASE("Comet orbital elements and SOI transitions during Mars encounter", "[ printf(" Distance to Sun: %.6f AU\n\n", change.distance_to_sun_au); } } - +#endif 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);