# Newton-Raphson Analytical Propagation - Implementation Plan

## Overview
Plan to replace RK4 numerical integration with Newton-Raphson analytical propagation for significantly larger simulation timesteps while maintaining accuracy.

## Motivation
### Current Limitations with RK4
- Time step constrained to seconds/minutes for stability
- Mercury orbiter (MESSENGER-like) limits stability to ~270s max dt
- Default dt=60s (only 22% of stability limit)
- Numerical drift accumulates over long simulations

### Benefits of Analytical Propagation
- Time steps of hours/days with perfect 2-body accuracy
- No numerical drift (exact solution to Kepler's problem)
- Newton-Raphson converges in 3-5 iterations (very fast)
- Enables much faster simulation of long-duration missions

## Proposed Solution

### Hybrid Approach
Use analytical propagation for orbital motion, numerical integration during burns:

1. **Normal operation (99% of time)**
   - Newton-Raphson solves Kepler's equation for true anomaly at time t
   - Direct conversion from orbital elements to state vectors
   - Perfect energy conservation

2. **During burns (<1% of time)**
   - Switch to numerical integration (RK4) for flexible timestep
   - Apply thrust acceleration combined with gravity
   - After burn, convert state vectors back to orbital elements
   - Resume analytical propagation

## Architecture

### Data Structure Changes

#### Spacecraft Structure (enhancements)
```cpp
struct Spacecraft {
    // Existing fields
    char name[64];
    double mass;
    int parent_index;
    OrbitalElements orbit;
    Vec3 global_position;
    Vec3 global_velocity;
    Vec3 local_position;
    Vec3 local_velocity;

    // New fields for analytical propagation
    bool in_active_burn;              // Currently executing finite-duration burn

    // Burn state
    double burn_start_time;
    double burn_duration;
    double delta_v_remaining;
    Vec3 burn_acceleration;           // Constant thrust acceleration vector
};
```

#### SimulationState Structure (enhancements)
```cpp
struct SimulationState {
    // Existing fields...

    // New propagation control
    double analytical_dt;             // Time step for analytical propagation (hours/days)
    double burn_dt;                   // Time step during burns (seconds)
};
```

## Implementation Phases

### Phase 1: Core Mathematical Functions
**Status**: Not started
**Estimated effort**: 4-6 hours

#### 1.1 Cartesian to Orbital Elements Conversion
**Function signature**:
```cpp
OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, double parent_mass);
```

**Algorithm**:
1. Calculate specific angular momentum: `h = r × v`
2. Calculate eccentricity vector: `e = ((v² - μ/r)r - (r·v)v) / μ`
3. Calculate eccentricity magnitude: `e = |e|`
4. Calculate semi-major axis:
   - `a = -μ / (2ε)` for elliptical orbits (ε < 0)
   - `a = μ / (2ε)` for hyperbolic orbits (ε > 0)
5. Calculate true anomaly from `r·e = r·e·cos(ν)`
6. Calculate inclination from `h_z = h·cos(i)`
7. Calculate longitude of ascending node from node vector `n = K × h`
8. Calculate argument of periapsis from `e·n = e·cos(ω)`

**Edge cases**:
- Circular orbits (e ≈ 0): Set true anomaly to 0
- Equatorial orbits (i ≈ 0): Set Ω = 0, ω = λ (true longitude)
- Hyperbolic orbits: Handle negative semi-major axis

#### 1.2 Newton-Raphson Solver for Kepler's Equation
**Function signature**:
```cpp
double solve_kepler_equation(double mean_anomaly, double eccentricity);
```

**Algorithm**:
```
Initial guess: E₀ = M + e·sin(M) + (e²/2)·sin(2M)
Iteration: Eₙ₊₁ = Eₙ - (Eₙ - e·sin(Eₙ) - M) / (1 - e·cos(Eₙ))
Convergence: |Eₙ₊₁ - Eₙ| < 1e-10 or max 50 iterations
```

**Initial Guess Formula**:
```cpp
inline double
getInitialTrialValue(double mean_anom, double ecc)
{
    return mean_anom + ecc * sin(mean_anom)
        + ((pow(ecc, 2) / 2) * sin(2 * mean_anom));
}
```

**Optimization**:
- Use series expansion initial guess for faster convergence
- Use hyperbolic Kepler equation for e > 1
- Cache convergence threshold based on precision needs

#### 1.3 Analytical Propagation Function
**Function signature**:
```cpp
void analytical_propagation_step(Spacecraft* craft, double time, double parent_mass);
```

**Algorithm**:
1. Calculate mean motion: `n = √(μ/a³)`
2. Calculate mean anomaly at time t: `M = n·(t - t₀) + M₀`
3. Solve Kepler's equation for eccentric anomaly E (Newton-Raphson)
4. Convert to true anomaly: `tan(ν/2) = √((1+e)/(1-e))·tan(E/2)`
5. Calculate radius: `r = a(1 - e²) / (1 + e·cos(ν))`
6. Calculate position in orbital plane (perifocal frame)
7. Apply 3D rotation matrices (same as existing `orbital_elements_to_cartesian`)
8. Calculate velocity from vis-viva equation or orbital velocity equations

### Phase 2: Hybrid Integration System
**Status**: Not started
**Estimated effort**: 6-8 hours

#### 2.1 Propagation Mode Selection
**Function signature**:
```cpp
void update_spacecraft_analytical(SimulationState* sim, Spacecraft* craft);
```

**Logic**:
```cpp
if (craft->in_active_burn) {
    // Use numerical integration during burn
    update_during_burn(sim, craft);
} else {
    // Use analytical propagation (default)
    analytical_propagation_step(craft, sim->time, parent_mass);
}
```

#### 2.2 Burn Execution with Numerical Integration
**Function signature**:
```cpp
void update_during_burn(SimulationState* sim, Spacecraft* craft);
```

**Algorithm**:
```
while (burn_in_progress):
    chunk_dt = min(sim->burn_dt, remaining_burn_time, time_until_soi_transition)

    // Get current state from orbital elements
    orbital_elements_to_cartesian(craft->orbit, parent_mass, &r, &v);

    // Combined acceleration: gravity + thrust
    Vec3 gravity = calculate_gravity(r, parent_mass);
    Vec3 total_accel = vec3_add(gravity, craft->burn_acceleration);

    // Numerical integration over chunk
    rk4_step_with_external_force(&r, &v, chunk_dt, craft->mass, parent_mass, total_accel);

    // Update orbital elements after chunk
    craft->orbit = cartesian_to_orbital_elements(r, v, parent_mass);

    // Update burn state
    craft->delta_v_remaining -= vec3_magnitude(craft->burn_acceleration) * chunk_dt;
    sim->time += chunk_dt;

    // Check burn completion
    if (craft->delta_v_remaining <= 0):
        craft->in_active_burn = false
        craft->orbit = cartesian_to_orbital_elements(r, v, parent_mass)
```

#### 2.3 RK4 with External Force (Enhancement)
**Function signature**:
```cpp
void rk4_step_with_external_force(Vec3* position, Vec3* velocity, double dt,
                                  double body_mass, double parent_mass,
                                  Vec3 external_acceleration);
```

**Algorithm**:
Same as existing `rk4_step()` but add external acceleration to each k_vel evaluation.

### Phase 3: SOI Transition Handling
**Status**: Infrastructure exists, needs adaptation
**Estimated effort**: 8-12 hours

#### 3.1 Orbital Element Transformation Across SOI Boundaries
**Function signature**:
```cpp
OrbitalElements transform_orbital_elements_across_soi(
    OrbitalElements old_elements,
    Vec3 position_global,
    Vec3 velocity_global,
    CelestialBody* new_parent,
    CelestialBody* old_parent
);
```

**Algorithm options**:

**Option A: Direct Conversion (Simpler)**
1. Convert old orbital elements to state vectors in global frame (already have)
2. Convert state vectors to orbital elements relative to new parent (Phase 1.1)
3. Requires position/velocity of both parents in global frame

**Option B: Lambert's Problem (More accurate)**
1. Solve Lambert's problem for trajectory between parents' positions
2. More complex but handles edge cases better
3. Useful for interplanetary transfers

**Recommended**: Start with Option A, implement Option B if needed

#### 3.2 Update Existing SOI Detection
**Modifications needed**:
```cpp
void update_soi_transitions(SimulationState* sim) {
    for (each spacecraft):
        if (craft crosses SOI boundary):
            // Transform orbital elements (analytical propagation)
            craft->orbit = transform_orbital_elements_across_soi(
                craft->orbit,
                craft->global_position,
                craft->global_velocity,
                new_parent,
                old_parent
            );
}
```

### Phase 4: Burn Command Interface
**Status**: Not started
**Estimated effort**: 4-6 hours

#### 4.1 Impulsive Burn Command
**Function signature**:
```cpp
void execute_impulsive_burn(Spacecraft* craft, Vec3 delta_v);
```

**Algorithm**:
```
1. Get current state from orbital elements
   orbital_elements_to_cartesian(craft->orbit, parent_mass, &r, &v);

2. Apply impulsive Δv
   v_new = v + delta_v

3. Convert back to orbital elements
   craft->orbit = cartesian_to_orbital_elements(r_new, v_new, parent_mass);
```

#### 4.2 Finite Duration Burn Command
**Function signature**:
```cpp
void start_continuous_burn(Spacecraft* craft, Vec3 thrust_acceleration, double duration);
```

**Algorithm**:
```
1. Set burn state
   craft->in_active_burn = true
   craft->burn_start_time = current_time
   craft->burn_duration = duration
   craft->burn_acceleration = thrust_acceleration
   craft->delta_v_remaining = |thrust_acceleration| * duration
```

### Phase 5: Testing and Validation
**Status**: Not started
**Estimated effort**: 8-12 hours

#### 5.1 Unit Tests
**Test cases**:
1. `cartesian_to_orbital_elements` conversion:
   - Circular orbits
   - Elliptical orbits
   - Parabolic orbits
   - Hyperbolic orbits
   - Equatorial orbits
   - Polar orbits
   - High inclination orbits

2. Newton-Raphson convergence:
   - Small eccentricities (e < 0.1)
   - Moderate eccentricities (0.1 < e < 0.5)
   - High eccentricities (e > 0.9)
   - Near-parabolic (e ≈ 1.0)
   - Hyperbolic (e > 1.0)

3. Analytical propagation accuracy:
   - Compare to RK4 for same orbits
   - Energy conservation over 1000 orbits
   - Period accuracy verification

#### 5.2 Integration Tests
**Test scenarios**:
1. Hohmann transfer:
   - Compare analytical vs. RK4 results
   - Verify orbital period match

2. Continuous thrust orbit raising:
   - Validate energy change
   - Check final orbit parameters

3. SOI transition with analytical propagation:
   - Earth-Moon transfer
   - Jupiter-Io transition

4. Long-duration simulation:
   - Multi-year Earth-Mars mission
   - Verify no numerical drift

#### 5.3 Performance Benchmarks
**Metrics**:
1. Time to simulate 1 Earth year with analytical vs. RK4
2. Newton-Raphson convergence rate (iterations vs. eccentricity)
3. Burn execution time (numerical phase)
4. Memory usage overhead

**Expected results**:
- 10-100x faster for large timesteps (hours/days)
- Negligible overhead for small timesteps
- Constant-time Newton-Raphson convergence

## Migration Strategy

### Phase A: Parallel Implementation (No Breaking Changes)
- Add new functions to `physics.h` and `physics.cpp`
- Keep existing `rk4_step()` unchanged (for burn integration)
- Both methods available simultaneously

### Phase B: Gradual Migration
- Enable analytical mode for test spacecraft
- Validate against existing RK4 results
- Update test configs to use analytical mode

### Phase C: Make Default
- After validation, make analytical propagation the default
- Keep RK4 available for burn integration and special cases (n-body perturbations)

## Technical Challenges

### Challenge 1: Numerical Precision with Large Timesteps
**Issue**: Floating-point errors may accumulate when jumping days/weeks
**Mitigation**:
- Use double precision (already using)
- Implement orbital element normalization after large jumps
- Consider splitting large timesteps into smaller chunks for precision

### Challenge 2: SOI Transition During Burn
**Issue**: What if burn crosses SOI boundary?
**Solutions**:
- Pause burn at SOI boundary, complete transition, resume burn
- Use combined acceleration during transition (numerical integration)
- Design burns to avoid SOI crossings (planning constraint)

### Challenge 3: Hyperbolic Trajectories
**Issue**: Hyperbolic Kepler equation different from elliptical
**Solution**:
- Implement hyperbolic Kepler solver: `H - e·sinh(H) = M`
- Detect orbit type from eccentricity
- Use appropriate solver based on orbit type

### Challenge 4: Eccentricity Near 1.0 (Parabolic)
**Issue**: Numerical instability at e ≈ 1.0
**Solution**:
- Treat parabolic as special case (semi-latus rectum)
- Use universal variable formulation for robustness
- Add tolerance band around e = 1.0

### Challenge 5: Continuous Thrust Optimization
**Issue**: Small burn chunks may be inefficient
**Solution**:
- Adaptive burn chunk sizing based on acceleration magnitude
- Larger chunks for low-thrust, smaller for high-thrust
- Cache intermediate calculations

## Performance Considerations

### Expected Performance Gains
| Scenario | RK4 dt | Analytical dt | Speedup |
|----------|--------|--------------|---------|
| Low Earth Orbit | 60s | 3600s (1 hour) | 60x |
| Geostationary Orbit | 60s | 3600s (1 hour) | 60x |
| Moon orbit | 60s | 86400s (1 day) | 1440x |
| Interplanetary | 60s | 172800s (2 days) | 2880x |

### Computational Cost Analysis
**Newton-Raphson per step**:
- 3-5 iterations
- Each iteration: trig functions, basic arithmetic
- Cost: ~100-200 FLOPs per step

**Comparison to RK4**:
- RK4: 4 force evaluations per step
- Each force evaluation: sqrt, division, vector operations
- Cost: ~50-80 FLOPs per force evaluation × 4 = ~200-320 FLOPs per step

**Conclusion**: Similar per-step computational cost, but analytical steps are 10-1000x larger

### Memory Overhead
- Minimal: Store orbital elements instead of position/velocity
- Already storing both in current implementation
- Negligible additional memory usage

## Dependencies
- None beyond current math library (cmath)
- Optional: Advanced orbital mechanics library for Lambert's problem (Phase 3.1 Option B)

## Risk Assessment

| Risk | Probability | Impact | Mitigation |
|------|-------------|--------|------------|
| Numerical instability at e ≈ 1.0 | Medium | High | Implement universal variable formulation |
| SOI transition errors | Low | High | Extensive testing with Moon/Phobos scenarios |
| Performance regression for small dt | Low | Low | Keep RK4 available, benchmark extensively |
| Burn integration accuracy | Medium | Medium | Adaptive timestep, validate against pure numerical |
| Complex implementation | High | Medium | Incremental phases, parallel implementation |

## Success Criteria

### Functional Requirements
- [ ] Newton-Raphson solves Kepler's equation for all eccentricity ranges
- [ ] Analytical propagation matches RK4 to within 1% for circular/elliptical orbits
- [ ] Impulsive burns correctly update orbital elements
- [ ] Continuous burns maintain numerical accuracy
- [ ] SOI transitions preserve orbital mechanics correctly

### Performance Requirements
- [ ] Analytical propagation is 10x faster than RK4 for dt > 1 hour
- [ ] Newton-Raphson converges in < 10 iterations for e < 0.99
- [ ] Memory overhead < 5% compared to RK4

### Quality Requirements
- [ ] Test coverage > 90% for new functions
- [ ] No regression in existing test suite
- [ ] Documentation updated for all new APIs

## References

### Algorithm References
1. "Fundamentals of Astrodynamics and Applications" - David Vallado
2. "Orbital Mechanics for Engineering Students" - Howard Curtis
3. "Methods of Orbit Determination" - Pedro Escobal

### Kepler's Equation Solvers
1. Newton-Raphson method with series expansion initial guess
2. Danby's method (higher convergence rate)
3. Universal variable formulation (handles all orbit types)

### Orbital Element Conversion
1. "Orbital Elements from State Vectors" - Vallado Chapter 2
2. "State Vectors from Orbital Elements" - existing implementation

## Future Enhancements

### Post-Implementation
1. Universal variable formulation (unifies elliptical/parabolic/hyperbolic)
2. Perturbations via Gauss's variational equations
3. Higher-order burn optimization (optimal control)
4. Real-time trajectory optimization
5. Monte Carlo uncertainty propagation

### Advanced Features
1. N-body perturbations with analytical corrections
2. Solar radiation pressure modeling
3. Atmospheric drag during low-thrust ascent
4. Multi-body gravity assists
5. Lunar descent powered flight

## Timeline Estimation

| Phase | Effort | Dependencies |
|------|--------|--------------|
| Phase 1: Core Math | 4-6 hours | None |
| Phase 2: Hybrid System | 6-8 hours | Phase 1 |
| Phase 3: SOI Handling | 8-12 hours | Phase 2 |
| Phase 4: Burn Interface | 4-6 hours | Phase 2 |
| Phase 5: Testing | 8-12 hours | Phases 1-4 |
| **Total** | **30-44 hours** | |

## Decision Points

### Before Starting Phase 1
- [ ] Confirm desired time step sizes (hours vs. days)
- [ ] Decide on hyperbolic/parabolic handling requirements
- [ ] Choose orbital element conversion algorithm (direct vs. Lambert)

### Before Starting Phase 3
- [ ] Validate Phase 2 burn execution accuracy
- [ ] Choose SOI transformation method (Option A vs. B)

### Before Phase 5
- [ ] Define performance benchmarks and acceptance criteria
- [ ] Identify critical test scenarios (Earth-Moon, Jupiter-Io, etc.)

## Open Questions

1. Should we implement universal variable formulation for robustness?
2. Do we need support for optimal control (continuous thrust optimization)?
3. What tolerance for Kepler's equation solver (1e-10 vs. 1e-12)?
4. Do we need to support non-Keplerian orbits (perturbed, n-body)?

---

**Document Status**: Planning - Not Implemented
**Last Updated**: Session: Time Step Stability Analysis
**Next Review**: When ready to begin implementation