claudes_game/docs/planning/newton_raphson_propagation_plan.md

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)

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)

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:

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:

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:

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:

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:

void update_spacecraft_analytical(SimulationState* sim, Spacecraft* craft);

Logic:

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:

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:

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:

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:

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:

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:

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