Refactor velocity calculation and consolidate initialization

- Rename compute_orbital_velocity_from_vis_viva to calc_orbital_velocity
- Change velocity function to return Vec3 instead of modifying body directly
- Further condense initialize_bodies into single loop with inlined coordinate init
- Remove redundant calculate_initial_velocities, calculate_soi_radii, and initialize_local_coordinates functions
- All initialization now handled efficiently in one loop

src/simulation.cpp

@@ -165,8 +165,9 @@ void update_simulation(SimulationState* sim) {
     sim->time += sim->dt;
 }
 
-static void compute_orbital_velocity_from_vis_viva(CelestialBody* body,
-                                                     CelestialBody* parent) {
+// 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;
@@ -193,86 +194,39 @@ static void compute_orbital_velocity_from_vis_viva(CelestialBody* body,
     }
 
     vel_dir = vec3_normalize(vel_dir);
-    body->velocity = vec3_scale(vel_dir, speed);
-    body->velocity = vec3_add(body->velocity, parent->velocity);
+    Vec3 velocity = vec3_scale(vel_dir, speed);
+    return vec3_add(velocity, parent->velocity);
 }
 
-void calculate_initial_velocities(SimulationState* sim) {
-    for (int i = 0; i < sim->body_count; i++) {
-        CelestialBody* body = &sim->bodies[i];
-
-        if (body->parent_index == -1) {
-            body->velocity = {0.0, 0.0, 0.0};
-        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
-            compute_orbital_velocity_from_vis_viva(body, parent);
-        }
-    }
+// 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
 }
 
-// Calculate SOI radii for all bodies
-// r_soi = a * (m/M)^(2/5) where a = semi-major axis, m = body mass, M = parent mass
 // 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;
 
-        // Initialize velocities
-        if (body->parent_index == -1) {
-            body->velocity = {0.0, 0.0, 0.0};
-        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
-            compute_orbital_velocity_from_vis_viva(body, parent);
-        }
-
-        // Calculate SOI radii
-        if (body->parent_index == -1) {
-            body->soi_radius = 1e15;
-        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
-            double mass_ratio = body->mass / parent->mass;
-            body->soi_radius = body->semi_major_axis * pow(mass_ratio, 0.4);
-        }
-
-        // Initialize local coordinates
-        if (body->parent_index == -1) {
-            body->local_position = body->position;
-            body->local_velocity = body->velocity;
-        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
+        // 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);
-        }
-    }
-}
-
-// Calculate SOI radii for all bodies (kept for individual testing)
-// r_soi = a * (m/M)^(2/5) where a = semi-major axis, m = body mass, M = parent mass
-void calculate_soi_radii(SimulationState* sim) {
-    for (int i = 0; i < sim->body_count; i++) {
-        CelestialBody* body = &sim->bodies[i];
-
-        if (body->parent_index == -1) {
-            body->soi_radius = 1e15;
-        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
-            double mass_ratio = body->mass / parent->mass;
-            body->soi_radius = body->semi_major_axis * pow(mass_ratio, 0.4);
-        }
-    }
-}
-
-void initialize_local_coordinates(SimulationState* sim) {
-    for (int i = 0; i < sim->body_count; i++) {
-        CelestialBody* body = &sim->bodies[i];
-
-        if (body->parent_index == -1) {
+            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;
-        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
-            CelestialBody* parent = &sim->bodies[body->parent_index];
-            body->local_position = vec3_sub(body->position, parent->position);
-            body->local_velocity = vec3_sub(body->velocity, parent->velocity);
+
+            // Root body (like Sun) has infinite SOI, use a large value
+            body->soi_radius = 1e15;    // 1000 AU in meters
         }
     }
 }

src/simulation.h

@@ -39,15 +39,8 @@ int add_body_to_simulation(SimulationState* sim, CelestialBody* body);
 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);
+double calculate_soi_radius(CelestialBody* body, CelestialBody* parent);
 
-// Velocity initialization
-void calculate_initial_velocities(SimulationState* sim);
-
-// SOI helpers
-void calculate_soi_radii(SimulationState* sim);
-
-// Coordinate frame management
-void initialize_local_coordinates(SimulationState* sim);
 void compute_global_coordinates(SimulationState* sim);
 
 // Combined initialization - sets velocities, SOI radii, and local coordinates in single loop