diff --git a/docs/hierarchical_frames_plan.md b/docs/hierarchical_frames_plan.md
new file mode 100644
index 0000000..0bfdea2
--- /dev/null
+++ b/docs/hierarchical_frames_plan.md
@@ -0,0 +1,285 @@
+# Hierarchical Coordinate Frames Implementation Plan
+
+## Goal
+Transform the simulation from global coordinate space to hierarchical local frames to:
+1. Improve numerical precision for moon orbits
+2. Isolate planetary perturbations from affecting moons
+3. Enable future patched conics implementation for satellites
+4. Support SOI transitions with proper coordinate transformations
+
+## Current System Analysis
+
+### Storage (Global Coordinates):
+- `CelestialBody.position` - absolute position in Sun-centered frame
+- `CelestialBody.velocity` - absolute velocity in Sun-centered frame
+- `CelestialBody.parent_index` - determines which body to calculate gravity from
+
+### Physics Integration:
+- `evaluate_acceleration()` - calculates gravity force from parent only (2-body approximation)
+- `rk4_step()` - integrates using global coordinates
+- SOI transitions change `parent_index`, but coordinates stay in global frame
+
+### Key Observation:
+Line 145 in simulation.cpp already composes velocities: `body->velocity = vec3_add(body->velocity, parent->velocity)`
+
+This indicates the system is **already partially thinking in local frames** during initialization, but integrates in global frame.
+
+## Proposed Architecture: "Local Frame Integration"
+
+### Data Structure: Dual Storage (Recommended)
+Store both local and global coordinates:
+```cpp
+struct CelestialBody {
+    char name[64];
+    double mass;
+    double radius;
+    
+    // NEW: Local coordinates (relative to parent)
+    Vec3 local_position;     // position relative to parent
+    Vec3 local_velocity;     // velocity relative to parent
+    
+    // Keep for rendering/backward compatibility
+    Vec3 position;           // global position (computed from local)
+    Vec3 velocity;           // global velocity (computed from local)
+    
+    double soi_radius;
+    int parent_index;
+    float color[3];
+    double eccentricity;
+    double semi_major_axis;
+};
+```
+
+**Advantages:**
+- Easy rendering (global positions readily available)
+- Easy SOI checks (global positions)
+- Clear separation of concerns
+- Memory negligible for ~14 bodies (48 bytes per body)
+
+### Physics Integration Flow
+
+```
+Current:
+  body (global) → RK4 in global coords → body (global)
+
+Proposed:
+  body (local) → RK4 in local coords → body (local) → compute global
+```
+
+### Update Order: Parent-First (Recommended)
+```cpp
+void update_simulation(SimulationState* sim) {
+    // 1. Update all root bodies (in their own frame = global)
+    for (int i = 0; i < sim->body_count; i++) {
+        if (sim->bodies[i].parent_index == -1) {
+            rk4_step_local(&sim->bodies[i], ctx, sim->dt);
+        }
+    }
+    
+    // 2. Compute global positions for roots
+    compute_global_coordinates_for_roots(sim);
+    
+    // 3. Update all child bodies (in parent's frame)
+    for (int i = 0; i < sim->body_count; i++) {
+        CelestialBody* body = &sim->bodies[i];
+        if (body->parent_index >= 0) {
+            // Check for SOI transition
+            int new_parent = find_dominant_body(sim, i);
+            if (new_parent != body->parent_index && new_parent != -1) {
+                transition_to_new_parent(body, body->parent_index, 
+                                        new_parent, sim);
+            }
+            
+            // Integrate in local frame
+            rk4_step_local(body, ctx, sim->dt);
+        }
+    }
+    
+    // 4. Compute global positions for all children
+    compute_global_coordinates_for_children(sim);
+    
+    sim->time += sim->dt;
+}
+```
+
+### SOI Transition with Frame Transform
+
+**Critical piece for satellites and patched conics:**
+
+```cpp
+void transition_to_new_parent(CelestialBody* body, int old_parent_idx, 
+                               int new_parent_idx, SimulationState* sim) {
+    // Current state is in old parent's frame
+    Vec3 old_local_pos = body->local_position;
+    Vec3 old_local_vel = body->local_velocity;
+    
+    // Get old and new parent
+    CelestialBody* old_parent = (old_parent_idx >= 0) 
+                                 ? &sim->bodies[old_parent_idx] : NULL;
+    CelestialBody* new_parent = &sim->bodies[new_parent_idx];
+    
+    // Compute global position (or use cached global coords)
+    Vec3 global_pos = old_parent ? vec3_add(old_local_pos, old_parent->position)
+                                  : old_local_pos;
+    Vec3 global_vel = old_parent ? vec3_add(old_local_vel, old_parent->velocity)
+                                  : old_local_vel;
+    
+    // Transform to new parent's frame
+    body->local_position = vec3_sub(global_pos, new_parent->position);
+    body->local_velocity = vec3_sub(global_vel, new_parent->velocity);
+    body->parent_index = new_parent_idx;
+}
+```
+
+## Implementation Phases
+
+### Phase 1: Foundation (No Behavior Change)
+**Goal:** Add local coordinate storage without changing physics
+
+**Tasks:**
+1. Add `local_position` and `local_velocity` to `CelestialBody`
+2. Add `initialize_local_coordinates()` - convert global→local on load
+3. Add `compute_global_coordinates()` - convert local→global after update
+4. Modify `update_simulation()` to call both functions
+5. **Verify:** All tests still pass (same behavior)
+
+**Files to modify:**
+- `src/simulation.h` (add fields)
+- `src/simulation.cpp` (add conversion functions)
+- `src/config_loader.cpp` (call initialize_local_coordinates)
+
+**Estimated complexity:** Low
+**Risk:** Very low (pure refactor, no logic change)
+
+**Expected outcome:**
+- ✅ Dual coordinate storage in place
+- ✅ No behavior change (all tests same status)
+- ✅ Foundation for local frame integration
+
+---
+
+### Phase 2: Local Frame Integration
+**Goal:** Actually integrate in local frames
+
+**Tasks:**
+1. Create `rk4_step_local()` - integrates using local coordinates
+2. Modify `evaluate_acceleration()` to work in local frame (parent at origin)
+3. Update `update_simulation()` to use local frame integration
+4. **Verify:** Moon tests improve (should see reduced drift)
+
+**Files to modify:**
+- `src/physics.cpp` (modify rk4_step, evaluate_acceleration)
+- `src/simulation.cpp` (update call sites)
+
+**Estimated complexity:** Medium
+**Risk:** Medium (core physics change)
+
+**Expected outcome:**
+- ✅ Moon drift issues **should be fixed** (improved numerical precision)
+- ✅ Test failures reduced (Moon, Io, Titan tests should pass)
+- ✅ Physics happens in local frames
+
+---
+
+### Phase 3: SOI Transition with Frame Transform
+**Goal:** Properly handle coordinate transformations during SOI crossings
+
+**Tasks:**
+1. Create `transition_to_new_parent()` function
+2. Modify SOI transition logic in `update_simulation()`
+3. Add tests for comet SOI transitions (Sun→Mars→Sun)
+4. **Verify:** Comet test still passes with smooth transitions
+
+**Files to modify:**
+- `src/simulation.cpp` (SOI transition logic)
+- `tests/test_comet_orbit.cpp` (verify transitions)
+
+**Estimated complexity:** Medium
+**Risk:** Medium (affects patched conics later)
+
+**Expected outcome:**
+- ✅ SOI transitions properly transform coordinates
+- ✅ Foundation for patched conics implementation
+- ✅ Comet transitions validated
+
+---
+
+### Phase 4: Parent-First Update Order
+**Goal:** Update hierarchy in correct order
+
+**Tasks:**
+1. Refactor `update_simulation()` to update roots first, then children
+2. Ensure parent global positions are current before children update
+3. **Verify:** No regression in tests
+
+**Files to modify:**
+- `src/simulation.cpp` (update_simulation)
+
+**Estimated complexity:** Low-Medium
+**Risk:** Low
+
+**Expected outcome:**
+- ✅ Hierarchical update order implemented
+- ✅ Parent positions current when updating children
+
+---
+
+### Phase 5: Validation & Optimization
+**Goal:** Ensure correctness and performance
+
+**Tasks:**
+1. Add test for frame transformations
+2. Profile performance (should be similar or better)
+3. Add documentation comments explaining coordinate systems
+4. Update implementation_plan.md
+
+**Files to modify:**
+- `tests/` (new frame transform tests)
+- `docs/implementation_plan.md`
+
+**Expected outcome:**
+- ✅ Fully validated hierarchical coordinate system
+- ✅ Documentation complete
+- ✅ Ready for satellite/spacecraft simulation
+
+## Design Decisions
+
+### Config Format
+Keep global positions in config (backward compatible). Convert to local coordinates on load via `initialize_local_coordinates()`.
+
+### Storage Strategy
+Dual storage (both local and global) for performance and simplicity.
+
+### Multi-level Hierarchies
+Not implementing at this time. Maximum 2 levels: Sun→Planet→Moon.
+Design allows future extension to Sun→Planet→Moon→Satellite if needed.
+
+### Timestep Strategy
+Single global timestep for entire simulation during initial implementation.
+Per-level timesteps deferred for future optimization.
+
+## Risk Assessment
+
+**Low Risk:**
+- Phase 1 (pure refactor, no logic change)
+- Phase 4 (update order change)
+
+**Medium Risk:**
+- Phase 2 (core physics change - but testable)
+- Phase 3 (frame transforms - complex but well-defined)
+
+**Mitigation:**
+- Implement phases incrementally with manual review after each phase
+- Keep old code commented for comparison
+- Add validation tests at each phase
+- Can roll back if tests regress
+
+## Expected Final Outcomes
+
+After all phases complete:
+- ✅ Moon orbital stability vastly improved (test failures fixed)
+- ✅ Numerical precision improved for nested orbits
+- ✅ SOI transitions with proper coordinate frame transformations
+- ✅ Foundation for patched conics and satellite simulation
+- ✅ Parent-first hierarchical update order
+- ✅ Fully documented coordinate system architecture
diff --git a/src/config_loader.cpp b/src/config_loader.cpp
index a34e678..3ce8930 100644
--- a/src/config_loader.cpp
+++ b/src/config_loader.cpp
@@ -133,6 +133,7 @@ bool load_system_config(SimulationState* sim, const char* filepath) {
 
     calculate_initial_velocities(sim);
     calculate_soi_radii(sim);
+    initialize_local_coordinates(sim);
 
     printf("Loaded %d bodies from %s\n", body_count, filepath);
     return true;
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 5c3aa7c..4480dd3 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -110,6 +110,8 @@ void update_simulation(SimulationState* sim) {
         }
     }
 
+    initialize_local_coordinates(sim);
+
     sim->time += sim->dt;
 }
 
@@ -184,6 +186,36 @@ void calculate_soi_radii(SimulationState* sim) {
     }
 }
 
+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->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);
+        }
+    }
+}
+
+void compute_global_coordinates(SimulationState* sim) {
+    for (int i = 0; i < sim->body_count; i++) {
+        CelestialBody* body = &sim->bodies[i];
+
+        if (body->parent_index == -1) {
+            body->position = body->local_position;
+            body->velocity = body->local_velocity;
+        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
+            CelestialBody* parent = &sim->bodies[body->parent_index];
+            body->position = vec3_add(body->local_position, parent->position);
+            body->velocity = vec3_add(body->local_velocity, parent->velocity);
+        }
+    }
+}
+
 OrbitalElements calculate_orbital_elements(CelestialBody* body, CelestialBody* primary,
                                            CelestialBody* optional_ref_body, double current_time) {
     const double AU = 1.496e11;
diff --git a/src/simulation.h b/src/simulation.h
index 1ecee79..c5255c1 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -8,8 +8,10 @@ struct CelestialBody {
     char name[64];
     double mass;              // kg
     double radius;            // meters
-    Vec3 position;            // meters from origin
-    Vec3 velocity;            // m/s
+    Vec3 local_position;      // position relative to parent (meters)
+    Vec3 local_velocity;      // velocity relative to parent (m/s)
+    Vec3 position;            // global position (meters from origin)
+    Vec3 velocity;            // global velocity (m/s)
     double soi_radius;        // sphere of influence radius (meters)
     int parent_index;         // index of gravitational parent (-1 for root body like Sun)
     float color[3];           // RGB color for rendering
@@ -41,6 +43,10 @@ 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);
+
 // Orbital elements calculation
 struct OrbitalElements {
     double time_days;