From 052efff5d1b70097ee1173e63d96d57b0dfe42b0 Mon Sep 17 00:00:00 2001
From: cinnaboot <cinnaboot@protonmail.com>
Date: Fri, 9 Jan 2026 10:01:11 -0500
Subject: [PATCH] Phase 2: Local frame integration (Earth Moon test now
 passing!)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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    | 22 +++++++++++++++-------
 src/simulation.cpp |  4 +++-
 2 files changed, 18 insertions(+), 8 deletions(-)

diff --git a/src/physics.cpp b/src/physics.cpp
index 0fd45e1..468548e 100644
--- a/src/physics.cpp
+++ b/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));
 }
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 4480dd3..9f152a4 100644
--- a/src/simulation.cpp
+++ b/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;
 }