Refactor: Break up update_simulation into helper functions

Extract large update_simulation function into focused helpers:
- update_bodies_physics() - Handle body RK4 and SOI transitions
- update_spacecraft_physics() - RK4 integration for spacecraft
- execute_pending_maneuvers() - Check triggers and execute burns
- compute_spacecraft_globals() - Global coords for spacecraft

Simplified update_simulation to call helpers sequentially:
1. Bodies physics
2. Bodies global coords
3. Spacecraft physics
4. Maneuver execution
5. Spacecraft global coords
6. Time step

Benefits:
- Single responsibility per function
- Cleaner, more readable code
- Easier to test individual components
- Net change: +96 lines (helper functions) -86 lines (removal)

All maneuver tests still passing.

docs/refactor_update_simulation.md

@@ -0,0 +1,62 @@
+# Refactoring Plan: Break Up update_simulation
+
+## Current State
+`update_simulation()` function is long with multiple responsibilities:
+1. Body RK4 updates and SOI transitions (~60 lines)
+2. Body global coordinate computation (~5 lines)
+3. Spacecraft RK4 updates (~15 lines)
+4. Maneuver execution checking (~20 lines)
+5. Spacecraft global coordinate computation (~15 lines)
+
+## Proposed Helper Functions
+
+### 1. `update_bodies_physics(SimulationState* sim)`
+- Handle body RK4 integration
+- Handle SOI transitions
+- ~60 lines of current body update logic
+
+### 2. `compute_body_globals(SimulationState* sim)`
+- Compute global coordinates for all bodies
+- Reuse existing function or inline it
+- ~5 lines
+
+### 3. `update_spacecraft_physics(SimulationState* sim)`
+- RK4 integration for spacecraft
+- ~15 lines
+
+### 4. `execute_pending_maneuvers(SimulationState* sim)`
+- Check maneuver triggers
+- Execute triggered maneuvers
+- ~20 lines
+
+### 5. `compute_spacecraft_globals(SimulationState* sim)`
+- Compute global coordinates for spacecraft
+- ~15 lines
+
+### 6. Simplified `update_simulation(SimulationState* sim)`
+```cpp
+void update_simulation(SimulationState* sim) {
+    update_bodies_physics(sim);
+    compute_body_globals(sim);
+    update_spacecraft_physics(sim);
+    execute_pending_maneuvers(sim);
+    compute_spacecraft_globals(sim);
+    sim->time += sim->dt;
+}
+```
+
+## Implementation Order
+1. Move body update logic to `update_bodies_physics()`
+2. Move spacecraft physics logic to `update_spacecraft_physics()`
+3. Create `execute_pending_maneuvers()` from existing maneuver loop
+4. Move spacecraft globals logic to `compute_spacecraft_globals()`
+5. Simplify `update_simulation()` to call helpers
+6. Update `simulation.h` with helper declarations
+7. Test to ensure behavior is identical
+
+## Benefits
+- Single responsibility for each function
+- Easier to test individual components
+- Cleaner `update_simulation()` that reads like a sequence
+- Easier to debug issues in specific areas
+- Functions can be reused or tested independently

src/simulation.cpp

@@ -159,6 +159,82 @@ void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_ax
 }
 
 void update_simulation(SimulationState* sim) {
+    update_bodies_physics(sim);
+    compute_global_coordinates(sim);
+    update_spacecraft_physics(sim);
+    execute_pending_maneuvers(sim);
+    compute_spacecraft_globals(sim);
+    sim->time += sim->dt;
+}
+
+// Calculate orbital velocity using vis-viva equation
+// Returns velocity vector for body relative to parent
+static Vec3 calc_orbital_velocity(CelestialBody* body, CelestialBody* parent) {
+    Vec3 r = vec3_sub(body->position, parent->position);
+    double distance = vec3_magnitude(r);
+    double e = body->eccentricity;
+    double a = body->semi_major_axis;
+    double v_squared;
+
+    if (fabs(e) < 0.0001) {
+        v_squared = G * parent->mass / a;
+    } else if (fabs(e - 1.0) < 0.0001) {
+        v_squared = 2.0 * G * parent->mass / distance;
+    } else {
+        v_squared = G * parent->mass * (2.0 / distance - 1.0 / a);
+    }
+
+    assert(v_squared >= 0);
+    double speed = (double) sqrt(v_squared);
+
+    Vec3 z_axis = {0.0, 0.0, 1.0};
+    Vec3 vel_dir = vec3_cross(z_axis, r);
+
+    if (vec3_magnitude(vel_dir) < 0.01) {
+        Vec3 x_axis = {1.0, 0.0, 0.0};
+        vel_dir = vec3_cross(z_axis, r);
+    }
+
+    vel_dir = vec3_normalize(vel_dir);
+    Vec3 velocity = vec3_scale(vel_dir, speed);
+    return vec3_add(velocity, parent->velocity);
+}
+
+// Calculate SOI radius for a single body
+// r_soi = a * (m/M)^(2/5) where a = semi-major axis, m = body mass, M = parent mass
+// Returns SOI radius in meters
+double calculate_soi_radius(CelestialBody* body, CelestialBody* parent) {
+    assert(body != nullptr && parent != nullptr);
+    double mass_ratio = body->mass / parent->mass;
+    return body->semi_major_axis * pow(mass_ratio, 0.4);  // 2/5 = 0.4
+}
+
+// Combined initialization - sets velocities, SOI radii, and local coordinates in single loop
+void initialize_bodies(SimulationState* sim) {
+    for (int i = 0; i < sim->body_count; i++) {
+        CelestialBody* body = &sim->bodies[i];
+        CelestialBody* parent = NULL;
+
+        // Set parent pointer if not root body
+        if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
+            parent = &sim->bodies[body->parent_index];
+            body->velocity = calc_orbital_velocity(body, parent);
+            body->local_position = vec3_sub(body->position, parent->position);
+            body->local_velocity = vec3_sub(body->velocity, parent->velocity);
+            body->soi_radius = calculate_soi_radius(body, parent);
+        } else { // root body
+            body->velocity = {0.0, 0.0, 0.0};
+            body->local_position = body->position;
+            body->local_velocity = body->velocity;
+
+            // Root body (like Sun) has infinite SOI, use a large value
+            body->soi_radius = 1e15;    // 1000 AU in meters
+        }
+    }
+}
+
+// Simulation update helper functions
+void update_bodies_physics(SimulationState* sim) {
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];
 
@@ -169,27 +245,22 @@ void update_simulation(SimulationState* sim) {
         int new_parent = find_dominant_body(sim, i);
 
         if (new_parent != body->parent_index) {
-            // Convert current local coordinates to global coordinates using old parent
             if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
                 CelestialBody* old_parent = &sim->bodies[body->parent_index];
                 body->position = vec3_add(body->local_position, old_parent->position);
                 body->velocity = vec3_add(body->local_velocity, old_parent->velocity);
             } else {
-                // old_parent is root (Sun): local = global
                 body->position = body->local_position;
                 body->velocity = body->local_velocity;
             }
 
-            // Update parent index
             body->parent_index = new_parent;
 
-            // Convert global coordinates to local coordinates using new parent
             if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
                 CelestialBody* new_parent_body = &sim->bodies[body->parent_index];
                 body->local_position = vec3_sub(body->position, new_parent_body->position);
                 body->local_velocity = vec3_sub(body->velocity, new_parent_body->velocity);
             } else {
-                // new_parent is root (Sun): global = local
                 body->local_position = body->position;
                 body->local_velocity = body->velocity;
             }
@@ -197,14 +268,14 @@ void update_simulation(SimulationState* sim) {
 
         if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
             CelestialBody* parent = &sim->bodies[body->parent_index];
+
             rk4_step(&body->local_position, &body->local_velocity,
                      sim->dt, body->mass, parent->mass);
         }
     }
+}
 
-    compute_global_coordinates(sim);
+void update_spacecraft_physics(SimulationState* sim) {
-    
-    // Update spacecraft
     for (int i = 0; i < sim->craft_count; i++) {
         Spacecraft* craft = &sim->spacecraft[i];
 
@@ -217,30 +288,32 @@ void update_simulation(SimulationState* sim) {
         rk4_step(&craft->local_position, &craft->local_velocity,
                  sim->dt, craft->mass, parent->mass);
     }
-    
+}
-    // Execute pending maneuvers (before computing globals)
+
+void execute_pending_maneuvers(SimulationState* sim) {
     for (int i = 0; i < sim->maneuver_count; i++) {
         Maneuver* maneuver = &sim->maneuvers[i];
-        
+
         if (maneuver->executed) {
             continue;
         }
-        
+
         if (maneuver->craft_index < 0 || maneuver->craft_index >= sim->craft_count) {
             continue;
         }
-        
+
         Spacecraft* craft = &sim->spacecraft[maneuver->craft_index];
-        
+
         if (check_maneuver_trigger(maneuver, craft, sim)) {
             execute_maneuver(maneuver, craft, sim->time);
         }
     }
-    
+}
-    // Compute spacecraft global coordinates (after maneuvers)
+
+void compute_spacecraft_globals(SimulationState* sim) {
     for (int i = 0; i < sim->craft_count; i++) {
         Spacecraft* craft = &sim->spacecraft[i];
-        
+
         if (craft->parent_index >= 0 && craft->parent_index < sim->body_count) {
             CelestialBody* parent = &sim->bodies[craft->parent_index];
             craft->position = vec3_add(parent->position, craft->local_position);
@@ -250,75 +323,6 @@ void update_simulation(SimulationState* sim) {
             craft->velocity = craft->local_velocity;
         }
     }
-    
-    sim->time += sim->dt;
-}
-
-
-// Calculate orbital velocity using vis-viva equation
-// Returns velocity vector for body relative to parent
-static Vec3 calc_orbital_velocity(CelestialBody* body, CelestialBody* parent) {
-    Vec3 r = vec3_sub(body->position, parent->position);
-    double distance = vec3_magnitude(r);
-    double e = body->eccentricity;
-    double a = body->semi_major_axis;
-    double v_squared;
-
-    if (fabs(e) < 0.0001) {
-        v_squared = G * parent->mass / a;
-    } else if (fabs(e - 1.0) < 0.0001) {
-        v_squared = 2.0 * G * parent->mass / distance;
-    } else {
-        v_squared = G * parent->mass * (2.0 / distance - 1.0 / a);
-    }
-
-    assert(v_squared >= 0);
-    double speed = (double) sqrt(v_squared);
-
-    Vec3 z_axis = {0.0, 0.0, 1.0};
-    Vec3 vel_dir = vec3_cross(z_axis, r);
-
-    if (vec3_magnitude(vel_dir) < 0.01) {
-        Vec3 x_axis = {1.0, 0.0, 0.0};
-        vel_dir = vec3_cross(z_axis, r);
-    }
-
-    vel_dir = vec3_normalize(vel_dir);
-    Vec3 velocity = vec3_scale(vel_dir, speed);
-    return vec3_add(velocity, parent->velocity);
-}
-
-// Calculate SOI radius for a single body
-// r_soi = a * (m/M)^(2/5) where a = semi-major axis, m = body mass, M = parent mass
-// Returns SOI radius in meters
-double calculate_soi_radius(CelestialBody* body, CelestialBody* parent) {
-    assert(body != nullptr && parent != nullptr);
-    double mass_ratio = body->mass / parent->mass;
-    return body->semi_major_axis * pow(mass_ratio, 0.4);  // 2/5 = 0.4
-}
-
-// Combined initialization - sets velocities, SOI radii, and local coordinates in single loop
-void initialize_bodies(SimulationState* sim) {
-    for (int i = 0; i < sim->body_count; i++) {
-        CelestialBody* body = &sim->bodies[i];
-        CelestialBody* parent = NULL;
-
-        // Set parent pointer if not root body
-        if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            parent = &sim->bodies[body->parent_index];
-            body->velocity = calc_orbital_velocity(body, parent);
-            body->local_position = vec3_sub(body->position, parent->position);
-            body->local_velocity = vec3_sub(body->velocity, parent->velocity);
-            body->soi_radius = calculate_soi_radius(body, parent);
-        } else { // root body
-            body->velocity = {0.0, 0.0, 0.0};
-            body->local_position = body->position;
-            body->local_velocity = body->velocity;
-
-            // Root body (like Sun) has infinite SOI, use a large value
-            body->soi_radius = 1e15;    // 1000 AU in meters
-        }
-    }
 }
 
 void compute_global_coordinates(SimulationState* sim) {

src/simulation.h

@@ -57,6 +57,12 @@ double calculate_soi_radius(CelestialBody* body, CelestialBody* parent);
 
 void compute_global_coordinates(SimulationState* sim);
 
+// Simulation update helper functions
+void update_bodies_physics(SimulationState* sim);
+void update_spacecraft_physics(SimulationState* sim);
+void execute_pending_maneuvers(SimulationState* sim);
+void compute_spacecraft_globals(SimulationState* sim);
+
 // Combined initialization - sets velocities, SOI radii, and local coordinates in single loop
 void initialize_bodies(SimulationState* sim);