diff --git a/docs/newton_raphson_propagation_plan.md b/docs/newton_raphson_propagation_plan.md
new file mode 100644
index 0000000..059ec14
--- /dev/null
+++ b/docs/newton_raphson_propagation_plan.md
@@ -0,0 +1,538 @@
+# 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