Phase 2: Local frame integration (Earth Moon test now passing!)

Modified physics integration to operate in local coordinate frames.
Parent bodies are treated as origin in child's reference frame.

Changes:
- rk4_step() now integrates local_position/local_velocity
- evaluate_acceleration() calculates gravity with parent at origin
- update_simulation() computes global coords FROM local after integration
- Added compute_global_coordinates() after root and child updates

Results:
- Earth Moon orbital stability test: NOW PASSING! (was failing)
- Moon drift reduced significantly (improved numerical precision)
- Test improvements: 3 failures → 2 failures (7/9 passing)
- Still failing: Io (Jupiter) and Titan (Saturn) - likely need tuning

Implements hierarchical_frames_plan.md Phase 2

src/physics.cpp

@@ -72,8 +72,8 @@ Vec3 calculate_acceleration(Vec3 force, double mass) {
 
 Vec3 evaluate_acceleration(Vec3 pos, Vec3 vel, AccelerationContext* ctx) {
     CelestialBody temp_body = *ctx->current_body;
-    temp_body.position = pos;
-    temp_body.velocity = vel;
+    temp_body.local_position = pos;
+    temp_body.local_velocity = vel;
 
     Vec3 total_force = {0.0, 0.0, 0.0};
 
@@ -89,7 +89,15 @@ Vec3 evaluate_acceleration(Vec3 pos, Vec3 vel, AccelerationContext* ctx) {
     } else {
         if (temp_body.parent_index >= 0 && temp_body.parent_index < ctx->sim->body_count) {
             CelestialBody* parent = &ctx->sim->bodies[temp_body.parent_index];
-            total_force = calculate_gravity_force(&temp_body, parent);
+            
+            double distance = vec3_magnitude(pos);
+            if (distance < 1.0) {
+                distance = 1.0;
+            }
+            
+            double force_magnitude = G * temp_body.mass * parent->mass / (distance * distance);
+            Vec3 direction = vec3_normalize(vec3_scale(pos, -1.0));
+            total_force = vec3_scale(direction, force_magnitude);
         }
     }
 
@@ -100,8 +108,8 @@ void rk4_step(CelestialBody* body, AccelerationContext* ctx, double dt) {
     Vec3 k1_vel, k2_vel, k3_vel, k4_vel;
     Vec3 k1_pos, k2_pos, k3_pos, k4_pos;
 
-    Vec3 pos0 = body->position;
-    Vec3 vel0 = body->velocity;
+    Vec3 pos0 = body->local_position;
+    Vec3 vel0 = body->local_velocity;
 
     k1_vel = evaluate_acceleration(pos0, vel0, ctx);
     k1_pos = vel0;
@@ -126,6 +134,6 @@ void rk4_step(CelestialBody* body, AccelerationContext* ctx, double dt) {
     Vec3 k_pos_sum = vec3_add(vec3_add(k1_pos, vec3_scale(k2_pos, 2.0)),
                               vec3_add(vec3_scale(k3_pos, 2.0), k4_pos));
 
-    body->velocity = vec3_add(vel0, vec3_scale(k_vel_sum, dt / 6.0));
-    body->position = vec3_add(pos0, vec3_scale(k_pos_sum, dt / 6.0));
+    body->local_velocity = vec3_add(vel0, vec3_scale(k_vel_sum, dt / 6.0));
+    body->local_position = vec3_add(pos0, vec3_scale(k_pos_sum, dt / 6.0));
 }

src/simulation.cpp

@@ -88,6 +88,8 @@ void update_simulation(SimulationState* sim) {
         }
     }
 
+    compute_global_coordinates(sim);
+
     for (int i = 0; i < sim->body_count; i++) {
         CelestialBody* body = &sim->bodies[i];
 
@@ -110,7 +112,7 @@ void update_simulation(SimulationState* sim) {
         }
     }
 
-    initialize_local_coordinates(sim);
+    compute_global_coordinates(sim);
 
     sim->time += sim->dt;
 }