# Patched Conics and SOI Transition Implementation Plan

**Date:** January 14, 2026
**Status:** Implementation Ready (Decisions Made)
**Branch:** patched-conics

**Decisions Made:**
1. ✅ Hysteresis: Adaptive approach (Option B)
2. ✅ Integration timing: Current approach (Option A), TODO for future
3. ✅ Test priorities: Create all 3 test cases first
4. ✅ Adaptive timesteps: Deferred to later work (Phase 5)

**Next Step:** Begin Phase 1 - Fix Root Body Transitions

## Overview

This plan implements support for patched conics trajectory simulation, enabling satellites to transition between multiple spheres of influence (SOI) in complex orbital scenarios:
- Planet → Star → Planet transfers
- Star → Planet → Moon rendezvous
- Multi-leg interplanetary missions

## Current State Analysis

### ✅ What's Already Working

1. **SOI Transitions Are Already Implemented**
   - Lines 103-121 in `simulation.cpp` show coordinate transformation logic
   - Converts local→global using old parent, then global→local using new parent
   - This is essentially Phase 3 implementation from hierarchical_frames_plan.md

2. **Local Frame Integration** (Phase 2) ✅
   - All bodies integrate in local coordinates
   - Global coordinates computed after each timestep
   - Improved numerical precision for nested orbits

3. **Parent-First Update Order** (Phase 4) ✅
   - Root bodies skipped in loop (fixed at origin)
   - Child bodies integrate using parent coordinates
   - Hierarchical update order implemented

### 🔴 Critical Issues for Patched Conics

**Issue 1: Cannot Transition to Root Bodies**
```cpp
if (new_parent != body->parent_index && new_parent != -1) {
    // ... transition logic
}
```
The `new_parent != -1` condition prevents switching to Sun (parent_index = -1). This breaks the scenario: Planet→Sun→Moon rendezvous is impossible.

**Issue 2: Hysteresis Barrier**
The 0.5x hysteresis factor in `find_dominant_body()` (line 71) creates one-way barriers:
- Planet→Sun: Possible (easy entry)
- Sun→Planet: **Impossible** (can't exit due to hysteresis)

**Issue 3: Integration After Transition**
Transition happens before integration in the same timestep, using coordinates from the end of the previous timestep. This may cause velocity discontinuities.

### ⚠️ Potential Issues

**Issue 4: Numerical Precision**
Satellites crossing between star/planet/moon scales will see position magnitude changes of 10⁸ to 10¹¹ meters, potentially losing precision.

**Issue 5: Fixed Timestep**
60s timestep may be too coarse for fast orbital phases (moon capture) and too slow for deep-space phases.

## Implementation Phases

### Phase 1: Fix Root Body Transitions (Critical)

**Goal:** Allow satellites to switch to/from Sun (parent_index = -1)

**Changes:**
1. Remove `new_parent != -1` check in `simulation.cpp` line 104
2. Add special handling for root body transitions:
   ```cpp
   if (new_parent != body->parent_index) {
       // Convert local → global using old parent
       if (body->parent_index >= 0) {
           // old_parent is a real body
           CelestialBody* old_parent = &sim->bodies[body->parent_index];
           body->position = vec3_add(body->local_position, old_parent->position);
           body->velocity = vec3_add(body->local_velocity, old_parent->velocity);
       } else {
           // old_parent is root (Sun): local = global
           body->position = body->local_position;
           body->velocity = body->local_velocity;
       }

       body->parent_index = new_parent;

       // Convert global → local using new parent
       if (new_parent >= 0) {
           // new_parent is a real body
           CelestialBody* new_parent_body = &sim->bodies[new_parent];
           body->local_position = vec3_sub(body->position, new_parent_body->position);
           body->local_velocity = vec3_sub(body->velocity, new_parent_body->velocity);
       } else {
           // new_parent is root (Sun): global = local
           body->local_position = body->position;
           body->local_velocity = body->velocity;
       }
   }
   ```

3. **Update `find_dominant_body()`** to properly handle -1 returns

**Files to modify:**
- `src/simulation.cpp` (lines 103-121)
- `src/simulation.h` (no changes needed)

**Estimated complexity:** Low
**Risk:** Medium (affects core transition logic)

**Expected outcome:**
- ✅ Satellites can transition to/from Sun
- ✅ Enables Planet→Sun→Planet transfers
- ✅ Enables Star→Planet→Moon rendezvous

**Tests to add:**
- Test satellite transitioning from Earth to Sun
- Test satellite transitioning from Sun to Mars
- Test full round-trip: Earth→Sun→Earth

---

### Phase 2: Adaptive Hysteresis (Critical for Round-Trips)

**Current Problem:**
The 0.5x hysteresis factor prevents oscillation but creates one-way barriers:
```cpp
if (distance < dist_to_current * 0.5) {
    dominant = i;
}
```

**Decision:** Adaptive Hysteresis (Option B) - **DECIDED ✅**

- Keep hysteresis but only apply when already in SOI
- Allow free switching when outside current SOI
- **Pros:** Prevents oscillation while enabling round-trips
- **Cons:** More complex logic

**Option B Implementation:**
```cpp
if (can_switch && i != dominant) {
    if (dominant == -1) {
        dominant = i;
    } else {
        CelestialBody* current_parent = &sim->bodies[dominant];
        double dist_to_current = vec3_distance(body->position, current_parent->position);

        if (outside_current_soi) {
            // Outside current SOI: switch to closest body (no hysteresis)
            if (distance < dist_to_current) {
                dominant = i;
            }
        } else {
            // Inside current SOI: apply hysteresis to prevent oscillations
            if (distance < dist_to_current * 0.5) {
                dominant = i;
            }
        }
    }
}
```

**Files to modify:**
- `src/simulation.cpp` (line 70-75)

**Estimated complexity:** Low-Medium
**Risk:** Medium (may affect transition behavior)

**Expected outcome:**
- ✅ Enables round-trip transitions (Earth→Sun→Earth)
- ✅ Maintains stability by preventing oscillation when inside SOI
- ✅ Allows free switching when outside current SOI

**Tests to add:**
- Validate Satellite→Planet→Satellite round-trip
- Validate Earth→Sun→Mars→Sun→Earth full round-trip

---

### Phase 3: Refactor Transition to Separate Function

**Goal:** Cleaner code, easier to test, better separation of concerns

**Current state:** Transition logic is inline in `update_simulation()` (lines 105-120)

**Proposed refactoring:**

**Add to `simulation.h`:**
```cpp
void transition_to_new_parent(SimulationState* sim, CelestialBody* body,
                           int old_parent_idx, int new_parent_idx);
```

**Add to `simulation.cpp`:**
```cpp
void transition_to_new_parent(SimulationState* sim, CelestialBody* body,
                           int old_parent_idx, int new_parent_idx) {
    // Current state is in old parent's frame
    Vec3 old_local_pos = body->local_position;
    Vec3 old_local_vel = body->local_velocity;

    // Convert to global coordinates using old parent
    if (old_parent_idx >= 0 && old_parent_idx < sim->body_count) {
        CelestialBody* old_parent = &sim->bodies[old_parent_idx];
        body->position = vec3_add(old_local_pos, old_parent->position);
        body->velocity = vec3_add(old_local_vel, old_parent->velocity);
    } else {
        // old_parent is root (Sun): local = global
        body->position = old_local_pos;
        body->velocity = old_local_vel;
    }

    // Update parent index
    body->parent_index = new_parent_idx;

    // Convert to local coordinates using new parent
    if (new_parent_idx >= 0 && new_parent_idx < sim->body_count) {
        CelestialBody* new_parent_body = &sim->bodies[new_parent_idx];
        body->local_position = vec3_sub(body->position, new_parent_body->position);
        body->local_velocity = vec3_sub(body->velocity, new_parent_body->velocity);
    } else {
        // new_parent is root (Sun): global = local
        body->local_position = body->position;
        body->local_velocity = body->velocity;
    }
}
```

**Update `update_simulation()`:**
```cpp
int new_parent = find_dominant_body(sim, i);
if (new_parent != body->parent_index) {
    transition_to_new_parent(sim, body, body->parent_index, new_parent);
}
```

**Files to modify:**
- `src/simulation.h` (add function declaration)
- `src/simulation.cpp` (extract and refactor transition logic)

**Estimated complexity:** Low
**Risk:** Low (pure refactor, no behavior change)

**Expected outcome:**
- ✅ Cleaner code with better separation of concerns
- ✅ Easier to unit test transition logic
- ✅ Follows Phase 3 plan from hierarchical_frames_plan.md

**Tests to add:**
- Unit tests for `transition_to_new_parent()` with all scenarios:
  - Body→Body transition
  - Body→Root transition
  - Root→Body transition
  - Root→Root transition (edge case)

---

### Phase 4: Multi-Body Transition Test Configs

**Goal:** Create realistic test scenarios for patched conics

**Decision:** Create all three test cases first, expect failures for unimplemented features - **DECIDED ✅**

**Implementation approach:**
1. Create all three test configurations
2. Run tests to establish baseline failures
3. Tests validate features as they're implemented
4. Comprehensive coverage from the start

#### Test Config 1: Satellite Rendezvous
**File:** `tests/configs/satellite_rendezvous.toml`

```toml
[[bodies]]
name = "Sun"
mass = 1.989e30
radius = 6.96e8
position = { x = 0.0, y = 0.0, z = 0.0 }
parent_index = -1
color = { r = 1.0, g = 1.0, b = 0.0 }
eccentricity = 0.0
semi_major_axis = 0.0

[[bodies]]
name = "Earth"
mass = 5.972e24
radius = 6.371e6
position = { x = 1.496e11, y = 0.0, z = 0.0 }
parent_index = 0
color = { r = 0.0, g = 0.5, b = 1.0 }
eccentricity = 0.0
semi_major_axis = 1.496e11

[[bodies]]
name = "Moon"
mass = 7.342e22
radius = 1.737e6
position = { x = 1.49984e11, y = 0.0, z = 0.0 }
parent_index = 1
color = { r = 0.7, g = 0.7, b = 0.7 }
eccentricity = 0.0
semi_major_axis = 3.844e8

[[bodies]]
name = "Satellite"
mass = 1.0e3
radius = 1.0e1
position = { x = 1.500e11, y = 1.0e8, z = 0.0 }
parent_index = 1
color = { r = 1.0, g = 0.0, b = 1.0 }
eccentricity = 0.2
semi_major_axis = 4.0e8
```

**Scenario:** Satellite launches from Earth, transfers to Moon, returns

#### Test Config 2: Interplanetary Transfer
**File:** `tests/configs/interplanetary_transfer.toml`

```toml
[[bodies]]
name = "Sun"
mass = 1.989e30
radius = 6.96e8
position = { x = 0.0, y = 0.0, z = 0.0 }
parent_index = -1
color = { r = 1.0, g = 1.0, b = 0.0 }
eccentricity = 0.0
semi_major_axis = 0.0

[[bodies]]
name = "Earth"
mass = 5.972e24
radius = 6.371e6
position = { x = 1.496e11, y = 0.0, z = 0.0 }
parent_index = 0
color = { r = 0.0, g = 0.5, b = 1.0 }
eccentricity = 0.0
semi_major_axis = 1.496e11

[[bodies]]
name = "Mars"
mass = 6.39e23
radius = 3.3895e6
position = { x = 2.279e11, y = 0.0, z = 0.0 }
parent_index = 0
color = { r = 0.8, g = 0.3, b = 0.1 }
eccentricity = 0.0
semi_major_axis = 2.279e11

[[bodies]]
name = "Probe"
mass = 1.0e3
radius = 1.0e1
position = { x = 1.496e11, y = 0.0, z = 0.0 }
parent_index = 1
color = { r = 0.0, g = 1.0, b = 0.0 }
eccentricity = 0.5
semi_major_axis = 1.888e11
```

**Scenario:** Probe: Earth→Sun→Mars

#### Test Config 3: Moon Capture
**File:** `tests/configs/moon_capture.toml`

```toml
[[bodies]]
name = "Sun"
mass = 1.989e30
radius = 6.96e8
position = { x = 0.0, y = 0.0, z = 0.0 }
parent_index = -1
color = { r = 1.0, g = 1.0, b = 0.0 }
eccentricity = 0.0
semi_major_axis = 0.0

[[bodies]]
name = "Jupiter"
mass = 1.898e27
radius = 6.9911e7
position = { x = 7.785e11, y = 0.0, z = 0.0 }
parent_index = 0
color = { r = 0.9, g = 0.7, b = 0.5 }
eccentricity = 0.0
semi_major_axis = 7.785e11

[[bodies]]
name = "Ganymede"
mass = 1.48e23
radius = 2.634e6
position = { x = 7.796e11, y = 0.0, z = 0.0 }
parent_index = 1
color = { r = 0.6, g = 0.6, b = 0.5 }
eccentricity = 0.0
semi_major_axis = 1.070e9

[[bodies]]
name = "Comet"
mass = 1.0e14
radius = 5.0e3
position = { x = 1.0e12, y = 0.0, z = 0.0 }
parent_index = 0
color = { r = 0.5, g = 0.8, b = 1.0 }
eccentricity = 2.0
semi_major_axis = -5.0e11
```

**Scenario:** Comet: Sun→Jupiter→Ganymede→Sun

**Files to create:**
- `tests/configs/satellite_rendezvous.toml`
- `tests/configs/interplanetary_transfer.toml`
- `tests/configs/moon_capture.toml`

**Estimated complexity:** Low
**Risk:** Low (configs only, no code changes)

**Expected outcome:**
- ✅ Realistic test scenarios for patched conics
- ✅ Covers multi-leg missions
- ✅ Tests root body transitions

---

### Phase 5: Adaptive Timestepping (Performance + Accuracy)

**Goal:** Different timesteps for different orbital phases

**Decision:** DEFERRED to later work - **DECIDED ✅**

**Problem:** Fixed 60s timestep is:
- Too coarse for fast orbital phases (moon capture)
- Too slow for deep-space phases

**Proposed solution:** Adaptive timestep based on orbital period

**Implementation:**
```cpp
double calculate_adaptive_timestep(CelestialBody* body, CelestialBody* parent) {
    if (parent == NULL || body->semi_major_axis <= 0.0) {
        return 60.0;  // Default timestep
    }

    // Calculate orbital period using Kepler's third law
    double T =2.0 * M_PI * sqrt(pow(body->semi_major_axis, 3) / (G * parent->mass));

    // Use 1/1000 of orbital period as timestep
    double adaptive_dt = T / 1000.0;

    // Clamp to reasonable bounds
    adaptive_dt = fmax(adaptive_dt, 10.0);   // Minimum 10s
    adaptive_dt = fmin(adaptive_dt, 600.0);  // Maximum 600s

    return adaptive_dt;
}
```

**Changes required:**
1. Add per-body timesteps to `SimulationState`
2. Update `update_simulation()` to use adaptive timesteps
3. Add synchronization mechanism (multiple timesteps)

**Complexity:** High
**Risk:** High (may affect energy conservation, requires careful testing)

**Expected outcome:**
- ✅ Better accuracy for fast orbits (moon capture)
- ✅ Faster simulation for deep-space phases
- ✅ Energy conserved across transitions

**Tests to add:**
- Verify energy drift with adaptive timesteps
- Verify orbital period accuracy with adaptive timesteps
- Test stability across SOI transitions

---

### Phase 6: Enhanced Debugging & Visualization

**Goal:** Better tools for debugging transitions and planning missions

**Features to add:**

#### 6.1 Transition Logging
```cpp
void log_transition(SimulationState* sim, CelestialBody* body,
                  int old_parent, int new_parent, double time) {
    printf("[%.2f days] %s: parent %d → %d\n",
           time / 86400.0, body->name, old_parent, new_parent);
}
```

#### 6.2 Visual SOI Boundaries
- Add wireframe spheres in renderer for SOI radii
- Toggle with keyboard key
- Use different colors for different bodies

#### 6.3 Trajectory Prediction
- Predict future trajectory for given time span
- Show predicted SOI crossings
- Assist with mission planning

**Files to modify:**
- `src/simulation.cpp` (logging)
- `src/renderer.cpp` (SOI visualization)
- `src/renderer.cpp` (trajectory prediction - new function)

**Estimated complexity:** Medium
**Risk:** Low (mostly UI/visualization)

**Expected outcome:**
- ✅ Better debugging tools for transitions
- ✅ Visual SOI boundaries
- ✅ Trajectory prediction for mission planning

---

## Decisions Made

### Decision 1: Hysteresis Strategy ✅
**Chosen:** Option B - Adaptive hysteresis approach

**Rationale:**
- Prevents oscillation while enabling round-trips
- Maintains stability during normal operation
- Allows free switching when outside current SOI
- Best balance between stability and flexibility

**Implementation:** Will use adaptive hysteresis in Phase 2 (lines 70-75 in simulation.cpp)

---

### Decision 2: Integration Timing ✅
**Chosen:** Option A - Keep current approach (transition at start of timestep)

**Rationale:**
- Simple, works for most cases
- Start with proven approach
- Defer optimization to future work

**TODO:** Consider half-step transitions (Option B) if velocity discontinuities become problematic

---

### Decision 3: Test Priorities ✅
**Chosen:** Option C - Implement all three test cases first

**Rationale:**
- Create all test configurations upfront
- Expect failures for unimplemented features
- Use tests as validation throughout implementation
- Comprehensive coverage from the start

**Implementation:**
1. Create all three test configs (Phase 4)
2. Run tests to establish baseline failures
3. Implement features to make tests pass
4. Re-run tests after each phase

**Test scenarios:**
- Satellite Rendezvous (Earth→Moon→Earth)
- Interplanetary Transfer (Earth→Sun→Mars)
- Moon Capture (Sun→Jupiter→Ganymede→Sun)

---

### Decision 4: Adaptive Timesteps Priority ✅
**Chosen:** Option B - Defer to later work

**Rationale:**
- Focus on SOI transitions first
- Fixed 60s timestep works adequately for testing
- Optimize performance and accuracy after core features complete

**TODO:** Implement adaptive timesteps in future update (Phase 5)
- Will address: Too coarse for fast orbits, too slow for deep space
- Estimated complexity: High
- Timeline: 2-3 days

---

## Success Criteria

### Phase 1 Success
- [ ] Satellite can transition to/from Sun
- [ ] Tests pass for Earth→Sun→Mars
- [ ] No energy conservation issues during transitions

### Phase 2 Success
- [ ] Round-trip transitions work (Earth→Sun→Earth)
- [ ] No oscillation at SOI boundaries
- [ ] Tests validate stability

### Phase 3 Success
- [ ] Transition logic extracted to separate function
- [ ] Unit tests cover all transition scenarios
- [ ] Code is cleaner and more maintainable

### Phase 4 Success
- [ ] Three test configs created
- [ ] Baseline test failures documented
- [ ] Configs are realistic and well-documented

### Phase 5 Success (Deferred)
- [ ] Adaptive timesteps implemented
- [ ] Energy drift < 1% with adaptive timesteps
- [ ] Performance improved for deep-space phases

### Phase 6 Success (Optional)
- [ ] Transition logging works
- [ ] SOI boundaries visualized
- [ ] Trajectory prediction functional

---

## Timeline Estimate

- **Phase 1:** 1-2 days (fix root body transitions)
- **Phase 2:** 1 day (adaptive hysteresis - decision made ✅)
- **Phase 3:** 0.5 days (refactoring)
- **Phase 4:** 1 day (test configs - create all three first)
- **Phase 5:** 2-3 days (adaptive timesteps - DEFERRED ⏸️)
- **Phase 6:** 1-2 days (debugging/visualization - optional)

**Total for core features (Phases 1-4):** 4-5 days
**Total with all features:** 8-10 days (including deferred Phases 5-6)

**Implementation approach:**
1. Create all three test configs (Phase 4)
2. Implement core transition features (Phases 1-3)
3. Validate with tests after each phase
4. Defer Phases 5-6 to future work

---

## Dependencies

### Required
- ✅ Hierarchical coordinate frames (complete)
- ✅ Local frame integration (complete)
- ✅ Parent-first update order (complete)
- ⏸️ Root body transitions (Phase 1 - pending)

### Optional
- ⏸️ Adaptive timesteps (Phase 5 - optional)
- ⏸️ Debugging tools (Phase 6 - optional)

---

## Risks and Mitigations

### High Risk
- **Energy conservation during transitions**
  - Mitigation: Careful testing, energy drift checks
  - Backup: Keep old code for rollback

- **Numerical precision across scales**
  - Mitigation: Use local frames (already implemented)
  - Backup: Double precision (already used)

### Medium Risk
- **Oscillation at SOI boundaries**
  - Mitigation: Adaptive hysteresis (Phase 2)
  - Backup: Increase hysteresis factor if needed

- **Timestep too coarse for fast orbits**
  - Mitigation: Adaptive timesteps (Phase 5 - optional)
  - Backup: Reduce fixed timestep if needed

### Low Risk
- **Code complexity increases**
  - Mitigation: Good unit tests, refactoring (Phase 3)
  - Backup: Keep functions small and focused

---

## References

- `docs/hierarchical_frames_plan.md` - Phase 3: SOI Transition with Frame Transform
- `docs/implementation_plan.md` - SOI Transition Mechanics section
- `src/simulation.cpp` - Current SOI transition implementation (lines 103-121)
- `tests/test_soi_transition.cpp` - Current SOI transition tests