Phase 4: Update initialization to use orbital mechanics module

docs/unified_orbital_elements_plan.md

@@ -245,32 +245,32 @@ struct Spacecraft {
    - For elliptical orbits (e < 1): `semi_major_axis > 0`
    - For hyperbolic orbits (e > 1): `semi_major_axis` negative or handle separately
 
-### Phase 4: Update Initialization
-
-10. **Replace `initialize_bodies()` with `initialize_orbital_objects()` (simulation.cpp)**
-    - For each body:
-      - If `parent_index >= 0`:
-        - Call `orbital_elements_to_cartesian(body->orbit, parent->mass, &local_pos, &local_vel)`
-        - Set `body->local_position = local_pos`, `body->local_velocity = local_vel`
-        - Compute global: `body->global_position = vec3_add(parent->global_position, local_pos)`
-        - Compute global: `body->global_velocity = vec3_add(parent->global_velocity, local_vel)`
-        - Calculate SOI: `body->soi_radius = calculate_soi_radius(body, parent)`
-      - Else (root body):
-        - Set all positions/velocities to zero
-        - Set `soi_radius = 1e15`
-
-    - For each spacecraft:
-      - Call `orbital_elements_to_cartesian(craft->orbit, parent->mass, &local_pos, &local_vel)`
-      - Set `craft->local_position = local_pos`, `craft->local_velocity = local_vel`
-      - Compute global: `craft->global_position = vec3_add(parent->global_position, local_pos)`
-      - Compute global: `craft->global_velocity = vec3_add(parent->global_velocity, local_vel)`
-
-11. **Remove old `calc_orbital_velocity()` function** (simulation.cpp)
-    - No longer needed, replaced by orbital mechanics module
-
-12. **Update all calls to `initialize_bodies()` → `initialize_orbital_objects()`**
-    - `config_loader.cpp`: Update call after config load
-    - Tests: Update if they call initialization directly
+### Phase 4: Update Initialization ✅ COMPLETE
+
+10. ✅ **Replace `initialize_bodies()` with `initialize_orbital_objects()` (simulation.cpp)**
+     - For each body:
+       - If `parent_index >= 0`:
+         - Call `orbital_elements_to_cartesian(body->orbit, parent->mass, &local_pos, &local_vel)`
+         - Set `body->local_position = local_pos`, `body->local_velocity = local_vel`
+         - Compute global: `body->global_position = vec3_add(parent->global_position, local_pos)`
+         - Compute global: `body->global_velocity = vec3_add(parent->global_velocity, local_vel)`
+         - Calculate SOI: `body->soi_radius = calculate_soi_radius(body, parent)`
+       - Else (root body):
+         - Set all positions/velocities to zero
+         - Set `soi_radius = 1e15`
+
+     - For each spacecraft:
+       - Call `orbital_elements_to_cartesian(craft->orbit, parent->mass, &local_pos, &local_vel)`
+       - Set `craft->local_position = local_pos`, `craft->local_velocity = local_vel`
+       - Compute global: `craft->global_position = vec3_add(parent->global_position, local_pos)`
+       - Compute global: `craft->global_velocity = vec3_add(parent->global_velocity, local_vel)`
+
+11. ✅ **Remove old `calc_orbital_velocity()` function** (simulation.cpp)
+     - No longer needed, replaced by orbital mechanics module
+
+12. ✅ **Update all calls to `initialize_bodies()` → `initialize_orbital_objects()`**
+     - `config_loader.cpp`: Update call after config load
+     - Tests: Update if they call initialization directly
 
 ### Phase 5: Update Test Configs
 
@@ -463,12 +463,11 @@ This is a **major breaking change** that affects all configs and code:
    - Help identify old vs new config formats (if we ever need backward compat)
 
 ## Success Criteria
-
 1. [ ] All test configs updated to new format
 2. [ ] Config loader successfully parses `orbit` tables
 3. [ ] `orbital_elements_to_cartesian()` correctly converts all orbit types
-4. [ ] `initialize_orbital_objects()` sets up bodies and spacecraft correctly
-5. [ ] All references to `position`/`velocity` renamed to `global_*`
+4. [x] `initialize_orbital_objects()` sets up bodies and spacecraft correctly
+5. [x] All references to `position`/`velocity` renamed to `global_*`
 6. [ ] Renderer displays orbits correctly
 7. [ ] Full test suite passes (`make test`)
 8. [ ] Solar system simulation runs correctly

src/config_loader.cpp

@@ -259,7 +259,7 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
 
     toml_free(result);
 
-    initialize_bodies(sim);
+    initialize_orbital_objects(sim);
 
     printf("Loaded %d bodies from %s\n", body_count, filepath);
     return true;

src/maneuver.cpp

@@ -63,7 +63,7 @@ void apply_impulsive_burn(Spacecraft* craft, BurnDirection direction, double del
 
 void apply_custom_burn(Spacecraft* craft, Vec3 delta_v_local) {
     craft->local_velocity = vec3_add(craft->local_velocity, delta_v_local);
-    craft->velocity = vec3_add(craft->velocity, delta_v_local);
+    craft->global_velocity = vec3_add(craft->global_velocity, delta_v_local);
 }
 
 double calculate_orbital_velocity(Spacecraft* craft, double parent_mass) {

src/simulation.cpp

@@ -1,6 +1,7 @@
 #include "simulation.h"
 #include "spacecraft.h"
 #include "maneuver.h"
+#include "orbital_mechanics.h"
 #include <cassert>
 #include <cstdlib>
 #include <cstring>
@@ -167,41 +168,6 @@ void update_simulation(SimulationState* 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->global_position, parent->global_position);
-    double distance = vec3_magnitude(r);
-    double e = body->orbit.eccentricity;
-    double a = body->orbit.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);
-
-    // FIXME: this whole section is just wrong and will break with inclined orbits
-    //  it's also a failure of our testing that we're not catching it
-    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->global_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
@@ -211,26 +177,46 @@ double calculate_soi_radius(CelestialBody* body, CelestialBody* parent) {
     return body->orbit.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) {
+// Initialize orbital objects from orbital elements
+// Converts orbital elements to local position/velocity and computes global coordinates
+void initialize_orbital_objects(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->global_velocity = calc_orbital_velocity(body, parent);
-            body->local_position = vec3_sub(body->global_position, parent->global_position);
-            body->local_velocity = vec3_sub(body->global_velocity, parent->global_velocity);
+
+            Vec3 local_pos, local_vel;
+            orbital_elements_to_cartesian(body->orbit, parent->mass, &local_pos, &local_vel);
+
+            body->local_position = local_pos;
+            body->local_velocity = local_vel;
+            body->global_position = vec3_add(parent->global_position, local_pos);
+            body->global_velocity = vec3_add(parent->global_velocity, local_vel);
             body->soi_radius = calculate_soi_radius(body, parent);
-        } else { // root body
+        } else {
+            body->local_position = {0.0, 0.0, 0.0};
+            body->local_velocity = {0.0, 0.0, 0.0};
+            body->global_position = {0.0, 0.0, 0.0};
             body->global_velocity = {0.0, 0.0, 0.0};
-            body->local_position = body->global_position;
-            body->local_velocity = body->global_velocity;
+            body->soi_radius = 1e15;
+        }
+    }
+
+    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];
+
+            Vec3 local_pos, local_vel;
+            orbital_elements_to_cartesian(craft->orbit, parent->mass, &local_pos, &local_vel);
 
-            // Root body (like Sun) has infinite SOI, use a large value
-            body->soi_radius = 1e15;    // 1000 AU in meters
+            craft->local_position = local_pos;
+            craft->local_velocity = local_vel;
+            craft->global_position = vec3_add(parent->global_position, local_pos);
+            craft->global_velocity = vec3_add(parent->global_velocity, local_vel);
         }
     }
 }
@@ -342,13 +328,13 @@ void compute_global_coordinates(SimulationState* sim) {
     }
 }
 
-OrbitalMetrics calculate_orbital_elements(CelestialBody* body, CelestialBody* primary,
+OrbitalAnalysis 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;
 
-    OrbitalMetrics elem;
+    OrbitalAnalysis elem;
     elem.time_days = current_time / SECONDS_PER_DAY;
 
     Vec3 r_vec = vec3_sub(body->global_position, primary->global_position);

src/simulation.h

@@ -71,11 +71,12 @@ 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);
+// Initialize orbital objects from orbital elements
+// Converts orbital elements to local position/velocity and computes global coordinates
+void initialize_orbital_objects(SimulationState* sim);
 
 // Orbital elements calculation (for output/analysis)
-struct OrbitalMetrics {
+struct OrbitalAnalysis {
     double time_days;
     double semi_major_axis_au;
     double eccentricity;
@@ -85,7 +86,7 @@ struct OrbitalMetrics {
     double velocity_magnitude;
 };
 
-OrbitalMetrics calculate_orbital_elements(CelestialBody* body, CelestialBody* primary,
+OrbitalAnalysis calculate_orbital_elements(CelestialBody* body, CelestialBody* primary,
                                            CelestialBody* optional_ref_body, double current_time);
 
 #endif

src/test_utilities.cpp

@@ -3,14 +3,14 @@
 #include <cmath>
 
 double calculate_kinetic_energy(CelestialBody* body) {
-    double v_squared = body->velocity.x * body->velocity.x +
-                       body->velocity.y * body->velocity.y +
-                       body->velocity.z * body->velocity.z;
+    double v_squared = body->global_velocity.x * body->global_velocity.x +
+                       body->global_velocity.y * body->global_velocity.y +
+                       body->global_velocity.z * body->global_velocity.z;
     return 0.5 * body->mass * v_squared;
 }
 
 double calculate_potential_energy_pair(CelestialBody* body1, CelestialBody* body2) {
-    double distance = vec3_distance(body1->position, body2->position);
+    double distance = vec3_distance(body1->global_position, body2->global_position);
     if (distance < 1.0) distance = 1.0;
     return -G * body1->mass * body2->mass / distance;
 }
@@ -33,9 +33,9 @@ double calculate_system_total_energy(SimulationState* sim) {
 OrbitalMetrics calculate_orbital_metrics(CelestialBody* body, CelestialBody* parent) {
     OrbitalMetrics metrics;
 
-    Vec3 relative_pos = vec3_sub(body->position, parent->position);
+    Vec3 relative_pos = vec3_sub(body->global_position, parent->global_position);
     metrics.orbital_radius = vec3_magnitude(relative_pos);
-    metrics.velocity_magnitude = vec3_magnitude(body->velocity);
+    metrics.velocity_magnitude = vec3_magnitude(body->global_velocity);
     metrics.angular_position = atan2(relative_pos.y, relative_pos.x);
 
     metrics.kinetic_energy = calculate_kinetic_energy(body);
@@ -80,7 +80,7 @@ void reset_orbit_tracker(OrbitTracker* tracker) {
 void update_orbit_tracker(OrbitTracker* tracker, CelestialBody* body, CelestialBody* parent, double current_time) {
     if (tracker->orbit_completed) return;
 
-    Vec3 relative_pos = vec3_sub(body->position, parent->position);
+    Vec3 relative_pos = vec3_sub(body->global_position, parent->global_position);
     double current_angle = atan2(relative_pos.y, relative_pos.x);
 
     if (tracker->quadrant_transitions == 0) {