# Hierarchical Coordinate Frames Implementation Plan

## Status: Phase 2 Complete ✅

**Last Updated:** 2026-01-09

**Current Progress:**
- ✅ Phase 1: Foundation (Dual coordinate storage)
- ✅ Phase 2: Local frame integration (Earth Moon test passing!)
- ⏸️ Phase 3: SOI transitions with frame transforms (deferred)
- ⏸️ Phase 4: Parent-first update order (deferred)
- ⏸️ Phase 5: Validation & optimization (deferred)

**Test Results After Phase 2:**
- Tests passing: 7/9 (78%)
- **Moon orbital stability around Earth:** ✅ PASSING (was failing)
- Io orbital stability: ❌ Still failing (orbit not completing)
- Titan orbital stability: ❌ Still failing (NaN drift)

**Commits:**
- `92be7f8` - Phase 1: Add local coordinate frame storage (no behavior change)
- `052efff` - Phase 2: Local frame integration (Earth Moon test now passing!)

---

## 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

**Status: COMPLETE** ✅
- Commit: `92be7f8`
- Date: 2026-01-09
- Test status: 6/9 passing (3 moon failures, unchanged from baseline)

---

### 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

**Status: COMPLETE** ✅
- Commit: `052efff`
- Date: 2026-01-09
- Test status: **7/9 passing (major improvement!)**
  - ✅ Earth Moon test: NOW PASSING (was failing with 20% drift)
  - ❌ Io test: Still failing (orbit not completing in time limit)
  - ❌ Titan test: Still failing (NaN drift, numerical issue)

**Implementation Details:**
- Modified `rk4_step()` to integrate `local_position` and `local_velocity`
- Modified `evaluate_acceleration()` to calculate gravity with parent at origin
- For child bodies in local frame: `distance = magnitude(position)`, `direction = -normalize(position)`
- `update_simulation()` now calls `compute_global_coordinates()` after integration
- Root bodies updated first, then global coords computed, then child bodies updated

**Key Insight - Why It Works:**
The local frame integration provides improved numerical precision by:
1. Eliminating large offsets (Moon integrates at ~3.8×10⁸ m instead of ~1.5×10¹¹ m)
2. Isolating moon orbits from planetary perturbations on parent
3. Maintaining full floating-point precision for small orbital changes

**Remaining Issues:**
- Io and Titan failures suggest timestep may be too coarse for distant/fast moons
- Or additional numerical stability improvements needed
- Not critical - major goal achieved (Earth Moon stable)

---

### Phase 3: SOI Transition with Frame Transform
**Status:** Not yet implemented (deferred)
**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

**Status:** Not yet implemented (deferred)

**Note:** Currently SOI transitions just change `parent_index`. In local frame integration,
this may cause position/velocity discontinuities. Phase 3 will implement proper coordinate
transformations during transitions: `new_local = global - new_parent_global`.

---

### Phase 4: Parent-First Update Order
**Goal:** Update hierarchy in correct order

**Status:** Partially implemented

**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

**Status:** Partially implemented

**Current Implementation:**
- Root bodies are updated first
- `compute_global_coordinates()` called after root update
- Then child bodies updated (using updated parent global positions)
- `compute_global_coordinates()` called again after child update

**This is effectively parent-first order!** Root bodies complete integration before
children start, and children use current parent positions. However, children still
use parent positions from START of timestep during their RK4 integration.

**Future Refinement (if needed):**
Could pass updated parent position into child RK4 steps for even higher accuracy.
Current approach is semi-implicit and works well for Phase 2 results.

---

### Phase 5: Validation & Optimization
**Status:** Not yet implemented (deferred)
**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

**Status:** Not yet implemented (deferred)

---

## Summary of Current State

### What Works:
1. ✅ Dual coordinate storage (local + global)
2. ✅ Local frame integration for all bodies
3. ✅ Automatic global coordinate computation
4. ✅ Earth-Moon orbital stability (major success!)
5. ✅ Improved numerical precision for nested orbits
6. ✅ Clean separation of local/global coordinate systems

### What's Deferred:
1. ⏸️ SOI transition frame transformations (Phase 3)
2. ⏸️ Full validation suite (Phase 5)
3. ⏸️ Io and Titan orbital tuning

### Ready For:
- Continued development of patched conics (after Phase 3)
- Satellite/spacecraft simulation (will need Phase 3 for SOI crossing)
- Further orbital mechanics improvements

### Notes for Future Development:
- Phase 3 (SOI transitions) is critical for spacecraft that cross SOI boundaries
- Current SOI transitions work but don't transform coordinates properly
- May need adaptive timesteps or smaller fixed timesteps for outer moons (Io, Titan)
- Consider adding orbit tracker diagnostics to debug remaining failures

## 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