Fix SOI calculation and orbit direction bugs

- 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 <noreply@anthropic.com>

src/bodies.cpp

@@ -1,6 +1,7 @@
 #include "bodies.h"
 #include <cstdlib>
 #include <cstring>
+#include <cstdio>
 #include <cmath>
 
 // Create a new simulation

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,
+                z_axis.y * r.z - z_axis.z * r.y,
-                r.z * z_axis.x - r.x * z_axis.z,
+                z_axis.z * r.x - z_axis.x * r.z,
-                r.x * z_axis.y - r.y * z_axis.x
+                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.x = x_axis.y * r.z - x_axis.z * r.y;
-                vel_dir.y = r.z * x_axis.x - r.x * x_axis.z;
+                vel_dir.y = x_axis.z * r.x - x_axis.x * r.z;
-                vel_dir.z = r.x * x_axis.y - r.y * x_axis.x;
+                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,
+                z_axis.y * r.z - z_axis.z * r.y,
-                r.z * z_axis.x - r.x * z_axis.z,
+                z_axis.z * r.x - z_axis.x * r.z,
-                r.x * z_axis.y - r.y * z_axis.x
+                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.x = x_axis.y * r.z - x_axis.z * r.y;
-                vel_dir.y = r.z * x_axis.x - r.x * x_axis.z;
+                vel_dir.y = x_axis.z * r.x - x_axis.x * r.z;
-                vel_dir.z = r.x * x_axis.y - r.y * x_axis.x;
+                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);
         }
     }
 }

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) {

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);