# Technical Reference - Orbital Mechanics Simulation

## Overview
N-body orbital mechanics simulator using analytical propagation for precise Keplerian trajectories. Supports elliptical, parabolic, and hyperbolic orbits with SOI (Sphere of Influence) transitions, impulsive burns, and 3D visualization using Raylib.

## Architecture
C-style C++ implementation (structs/functions, no classes/templates) with modular design. The system is organized into these interconnected modules: Main (entry point), Simulation (core loop), Orbital Mechanics (Keplerian propagation), Physics (vector/matrix math), Maneuver (impulsive burns), Spacecraft (craft logic), Config Validator (validation), Test Utilities (testing helpers), Renderer (3D visualization), UI Renderer (2D panels), and Config Loader (TOML parsing). Main depends on simulation, renderer, and UI renderer. Simulation module coordinates physics calculations and depends on orbital mechanics, maneuver, spacecraft, and config loader modules. Renderer depends on config loader for configuration. See module sections below for detailed descriptions.

## Core Data Structures

### Vec3 (physics.h)
3D vector for position, velocity, and acceleration
```cpp
struct Vec3 {
    double x, y, z;
};
```

### Mat3 (physics.h)
3x3 row-major matrix for coordinate transformations
```cpp
struct Mat3 {
    double m00, m01, m02;
    double m10, m11, m12;
    double m20, m21, m22;
};
```

### OrbitalElements (orbital_mechanics.h)
Keplerian orbital elements using union for parabolic/hyperbolic distinction
```cpp
struct OrbitalElements {
    union {
        double semi_major_axis;       // elliptical (e<1) and hyperbolic (e>1)
        double semi_latus_rectum;     // parabolic (e≈1)
    };
    double eccentricity;
    double true_anomaly;
    double inclination;
    double longitude_of_ascending_node;
    double argument_of_periapsis;
};
```

### CelestialBody (simulation.h)
Planet/moon with local and global coordinate frames
```cpp
struct CelestialBody {
    char name[64];
    double mass;              // kg
    double radius;            // meters
    int parent_index;         // index of gravitational parent (-1 for root)
    float color[3];           // RGB color for rendering

    OrbitalElements orbit;    // Keplerian elements from config

    Vec3 global_position;     // meters from origin
    Vec3 global_velocity;     // m/s
    Vec3 local_position;      // meters from parent
    Vec3 local_velocity;      // m/s relative to parent

    double soi_radius;        // sphere of influence radius (meters)
};
```

### Spacecraft (spacecraft.h)
Spacecraft similar to CelestialBody but without SOI or radius
```cpp
struct Spacecraft {
    char name[64];
    double mass;
    int parent_index;

    OrbitalElements orbit;

    Vec3 global_position;
    Vec3 global_velocity;
    Vec3 local_position;
    Vec3 local_velocity;
};
```

### Maneuver (maneuver.h)
Impulsive burn with execution trigger
```cpp
struct Maneuver {
    char name[64];
    int craft_index;
    BurnDirection direction;
    double delta_v;
    TriggerType trigger_type;
    double trigger_value;
    bool executed;
    double executed_time;
};
```

### SimulationState (simulation.h)
Top-level simulation container
```cpp
struct SimulationState {
    CelestialBody* bodies;
    int body_count;
    int max_bodies;

    Spacecraft* spacecraft;
    int craft_count;
    int max_craft;

    Maneuver* maneuvers;
    int maneuver_count;
    int max_maneuvers;

    double time;
    double dt;
    char config_name[256];
};
```

## Coordinate Frame System

### Local Frame
- Position/velocity relative to parent body
- Used for orbital mechanics calculations
- Primary coordinate system for propagation

### Global Frame
- Position/velocity from simulation origin (Sun at 0,0,0)
- Computed each frame: `global = parent.global + local`
- Used for SOI distance calculations and rendering

### Benefits
- Precision: Local coordinates use full double precision for small orbits (LEO, moon)
- Clarity: Separates orbit physics from system-wide positions
- Efficiency: SOI checks use global distances, propagation uses local
- Flexibility: Easy to add/remove bodies by updating parent_index

### Conversion Functions
- Local to global: `vec3_add(parent.global, local)`
- Global to local: `vec3_sub(global, parent.global)`

## Key Modules

### 1. Physics Module (physics.h/cpp)
Vector and matrix math utilities for orbital mechanics.

**Constants:**
- Gravitational constant: `G` in physics.h

**Vector Functions:**
- `vec3_add`, `vec3_sub`, `vec3_cross`, `vec3_scale`
- `vec3_magnitude`, `vec3_distance`, `vec3_normalize`, `vec3_dot`

**Matrix Functions:**
- `mat3_identity`, `mat3_multiply`, `mat3_multiply_vec3`
- `mat3_rotation_x`, `mat3_rotation_z`
- `mat3_rotation_orbital(omega, i, Omega)` - z-x-z Euler angles

**Physics Functions:**
- `calculate_acceleration(force, mass)`
- `evaluate_acceleration(relative_pos, body_mass, parent_mass)`
- `rk4_step()` - Available but NOT used in simulation (analytical used instead)

### 2. Orbital Mechanics Module (orbital_mechanics.h/cpp)
Keplerian orbit propagation and element conversion.

**Constants:**
- Parabolic tolerance: `PARABOLIC_TOLERANCE` in orbital_mechanics.h
- Kepler solver tolerance: `KEPLER_TOLERANCE` in orbital_mechanics.cpp

**Element Conversion:**
- `orbital_elements_to_cartesian(elements, parent_mass, out_pos, out_vel)` (orbital_mechanics.cpp:6-44)
  - Converts Keplerian elements to position/velocity vectors
  - Uses z-x-z Euler rotation: R_z(Ω) · R_x(i) · R_z(ω)

- `cartesian_to_orbital_elements(position, velocity, parent_mass)` (orbital_mechanics.cpp:186-299)
  - Reconstructs elements from state vectors
  - Handles near-parabolic (|ε| < 1e-10), near-circular (e < 1e-10) cases
  - Returns semi_latus_rectum for parabolic orbits

**Kepler Equation Solvers:**
- `solve_kepler_elliptical(mean_anomaly, eccentricity)`
  - Newton-Raphson: E - e·sin(E) = M

- `solve_kepler_hyperbolic(mean_anomaly, eccentricity)`
  - Newton-Raphson: H - e·sinh(H) = M

- `solve_barker_equation(mean_anomaly)`
  - Cubic solution: D + D³/3 = M where D = tan(ν/2)

**Anomaly Conversions:**
- `mean_anomaly_to_true_anomaly(mean_anomaly, eccentricity)`
  - Unified dispatcher for elliptical/hyperbolic

- `eccentric_to_true_anomaly(eccentric_anomaly, eccentricity)`
  - Tangent half-angle formula: tan(ν/2) = √((1+e)/(1-e)) · tan(E/2)
  - Near-parabolic special case for stability

- `hyperbolic_to_true_anomaly(hyperbolic_anomaly, eccentricity)`
  - tan(ν/2) = √((e+1)/(e-1)) · tanh(H/2)

**Primary Propagation Function:**
- `propagate_orbital_elements(elements, dt, parent_mass)` (orbital_mechanics.cpp:301-375)
  - **Critical: Main propagation method used by simulation**
  - Dispatches to parabolic/elliptical/hyperbolic based on eccentricity
  - Parabolic: Uses Barker's equation with mean motion n = √(μ/p³)
  - Elliptical: Kepler's equation with n = √(μ/a³)
  - Hyperbolic: Hyperbolic Kepler equation with n = √(μ/(-a)³)

### 3. Simulation Module (simulation.h/cpp)
Main simulation loop, SOI management, coordinate updates.

**SOI Calculation:**
- `calculate_soi_radius(body, parent)`
  - Hill sphere: r_soi = a · (m/M)^(2/5)

- `update_soi(body, parent, semi_major_axis)`
  - Updates body->soi_radius, sets 1e15 for root bodies

**SOI Transition Logic:**
- `find_dominant_body(sim, body_index)` (simulation.cpp:105-148)
  - For non-root bodies: checks if still within parent's SOI
  - For root bodies: finds closest body within SOI
  - Returns new parent index or 0 (Sun)

**Main Update Loop:**
- `update_simulation(sim)`
  - Sequence: bodies_physics → global_coords → spacecraft_physics → maneuvers → spacecraft_globals → time++

**Body Physics:**
- `update_bodies_physics(sim)`
  - Checks for SOI transitions each frame
  - On transition: reconstructs orbital elements for new parent
  - **Velocity deviation detection:** rebuilds elements if |v_local - v_expected| > 1e-6 (simulation.cpp:267-270)
  - Propagates using `propagate_orbital_elements()`
  - Updates local position/velocity from new elements

**Spacecraft Physics:**
- `update_spacecraft_physics(sim)`
  - Same velocity deviation detection as bodies (simulation.cpp:291-294)
  - Propagates and updates coordinates

**Coordinate Management:**
- `compute_global_coordinates(sim)`
  - Updates global positions from parent.global + local

- `compute_spacecraft_globals(sim)`
  - Same for spacecraft

**Initialization:**
- `initialize_orbital_objects(sim)`
  - Converts orbital elements to local coordinates for all bodies/craft
  - Calculates SOI radii for all bodies

**Dynamic Management:**
- `create_simulation(max_bodies, max_craft, max_maneuvers, dt)`
- `destroy_simulation(sim)`
- `add_body_to_simulation(sim, body)`
- `add_spacecraft(sim, craft)`

### 4. Maneuver Module (maneuver.h/cpp)
Impulsive burn execution with various burn directions.

**Burn Directions:**
- `BURN_PROGRADE` - Along velocity vector
- `BURN_RETROGRADE` - Opposite velocity
- `BURN_NORMAL` - Along angular momentum (orbit normal)
- `BURN_ANTINORMAL` - Opposite angular momentum
- `BURN_RADIAL_IN` - Toward parent
- `BURN_RADIAL_OUT` - Away from parent
- `BURN_CUSTOM` - User-specified vector

**Trigger Types:**
- `TRIGGER_TIME` - Execute at simulation time >= trigger_value
- `TRIGGER_TRUE_ANOMALY` - Execute when true anomaly within 0.01 rad of trigger_value

**Direction Calculation:**
- `get_burn_direction_vector(direction, local_pos, local_vel)`
  - Dispatches to specific direction calculators

**Burn Application:**
- `apply_impulsive_burn(craft, direction, delta_v)`
  - Instant velocity change: v_new = v_old + direction_vector · delta_v

- `apply_custom_burn(craft, delta_v_local)`
  - Applies arbitrary delta-v vector

**Trigger Checking:**
- `check_maneuver_trigger(maneuver, craft, sim)`
  - Time trigger: sim->time >= trigger_value
  - True anomaly trigger: computes ν from state, checks angular distance < 0.01 rad

**Execution:**
- `execute_maneuver(maneuver, craft, sim, current_time)` (maneuver.cpp:166-176)
  - Applies burn
  - **Reconstructs orbital elements** from new velocity
  - Marks executed, records time

### 5. Config Loader Module (config_loader.h/cpp)
TOML parsing for simulation configuration.

**Main Function:**
- `load_system_config(sim, filepath)` - Parses TOML file, loads bodies, spacecraft, maneuvers

**Config Format:**
- TOML format with [[bodies]], [[spacecraft]], [[maneuvers]] arrays
- See `docs/config_format.md` for complete documentation including all fields, validation rules, and examples

### 6. Config Validator Module (config_validator.h/cpp)
System-level validation to prevent simulation errors.

**Validation Functions:**
- `run_all_config_validations(sim)` - Master validation function
- See `docs/config_format.md` for all validation rules and constants

### 7. Renderer Module (renderer.h/cpp)
3D visualization using Raylib. See `docs/rendering.md` for detailed reference.

**Key Features:**
- Simulation XY plane → Render XZ plane transformation (90° rotation around X-axis)
- Linear distance/size scaling (1e-9 factor: 1 render unit = 1 billion meters)
- Relative rendering when body selected (children rendered around origin)
- Orbit path rendering (elliptical, parabolic, hyperbolic with appropriate segment counts)
- Wireframe spheres for bodies
- Child indicators for spacecraft/off-screen bodies
- Camera follow mode with orbit controls

### 8. UI Renderer Module (ui_renderer.h/cpp)
2D UI panels using raygui.

**Panels:**
- Info panel (bottom-left): Sim time, body count, FPS, controls, config name
- Objects list (top-left): All bodies and spacecraft with selection
- Info panel (top-right): Selected object details
- Maneuver list (below info panel): Pending/executed maneuvers

### 9. Test Utilities Module (test_utilities.h/cpp)
Energy calculations and orbit tracking for testing.

**OrbitalMetrics:**
- kinetic_energy, potential_energy, total_energy
- orbital_radius, velocity_magnitude
- angular_position

**OrbitTracker:**
- Tracks orbit completion with quadrant transitions
- 3D angle calculation using orbital elements
- `create_orbit_tracker_3d()` for inclined orbits

**Energy Functions:**
- `calculate_kinetic_energy(body)`
- `calculate_potential_energy_pair(body1, body2)`
- `calculate_system_total_energy(sim)`

**Orbit Tracking:**
- `create_orbit_tracker(body_index)`
- `update_orbit_tracker(tracker, body, parent, current_time)`
- `orbit_completed` flag, `time_at_completion` record

## Critical Implementation Details

### Analytical Propagation Algorithm

**Function:** `propagate_orbital_elements(elements, dt, parent_mass)` (orbital_mechanics.cpp:301-375)

**Elliptical Orbits (e < 1):**
1. Compute mean motion: n = √(μ/a³)
2. Convert true anomaly to eccentric anomaly: E = 2·atan(√((1-e)/(1+e)) · tan(ν/2))
3. Compute mean anomaly: M = E - e·sin(E)
4. Advance mean anomaly: M_new = M + n·dt
5. Solve Kepler equation for E_new using Newton-Raphson (max 50 iterations, tolerance 1e-10)
6. Convert back to true anomaly: ν_new = 2·atan(√((1+e)/(1-e)) · tan(E_new/2))
7. Return elements with updated true_anomaly

**Parabolic Orbits (|e-1| < PARABOLIC_TOLERANCE):**
1. Compute variable: D = tan(ν/2)
2. Compute Barker's mean anomaly: M = D + D³/3
3. Compute mean motion: n = √(μ/p³) where p = semi_latus_rectum
4. Advance: M_new = M + n·dt
5. Solve Barker's equation: D_new from cubic formula
6. Convert: ν_new = 2·atan(D_new)
7. Return elements with updated true_anomaly

**Hyperbolic Orbits (e > 1):**
1. Compute mean motion: n = √(μ/(-a)³) (note negative a)
2. Convert true anomaly to hyperbolic anomaly: H = 2·atanh(√((e-1)/(e+1)) · tan(ν/2))
3. Compute mean anomaly: M = e·sinh(H) - H
4. Advance: M_new = M + n·dt
5. Solve hyperbolic Kepler equation for H_new using Newton-Raphson (max 50 iterations, tolerance 1e-10)
6. Convert back: ν_new = 2·atan(√((e+1)/(e-1)) · tanh(H_new/2))
7. Return elements with updated true_anomaly

**Key Benefits:**
- Exact solution to Kepler's laws (no integration error)
- Energy conservation guaranteed
- Efficient for large time steps
- Handles all orbit types uniformly

### SOI Transition Logic

**Function:** `find_dominant_body(sim, body_index)` (simulation.cpp:105-148)

**Algorithm:**
1. If parent is not root (index != 0):
   - Check if distance to parent < parent.soi_radius
   - If yes: stay with current parent
   - If no: transition to root (Sun)

2. If parent is root (Sun):
   - Iterate through all other bodies
   - Find closest body where distance < body.soi_radius
   - Switch to that body (if found, else stay with root)

**Transition Handling:**
1. When `find_dominant_body()` returns new parent:
   - Compute current global position/velocity from old parent
   - Update `body->parent_index`
   - Compute new local coordinates relative to new parent
   - **Reconstruct orbital elements:** `cartesian_to_orbital_elements()`

**Orbital Element Reconstruction:**
- Converts local position/velocity to new Keplerian elements
- Handles edge cases (near-circular, near-parabolic)
- Ensures subsequent propagation is consistent with new frame

### Orbital Element Reconstruction After Burns

**Location:** `execute_maneuver()` (maneuver.cpp:166-176)

**Process:**
1. Apply impulsive burn to `craft->local_velocity`
2. **Reconstruct elements:** `cartesian_to_orbital_elements(local_pos, local_vel, parent_mass)`
3. New elements are used for subsequent analytical propagation

**Velocity Deviation Detection:**
- Before each propagation, check: `|v_local - v_expected_from_elements| > 1e-6` (simulation.cpp:267-270, 291-294)
- If exceeded: reconstruct elements from current state
- This catches numerical drift and ensures consistency

### 3D Orbital Orientation

**Rotation:** z-x-z Euler angles via `mat3_rotation_orbital(omega, i, Omega)`

**Sequence:**
1. Rotate by argument_of_periapsis (ω) around Z-axis
2. Rotate by inclination (i) around X-axis
3. Rotate by longitude_of_ascending_node (Ω) around Z-axis

**Matrix:** R_total = R_z(Ω) · R_x(i) · R_z(ω)

**Purpose:**
- Orient orbital plane in 3D space
- Align periapsis direction with eccentricity vector
- Define ascending node location

**Element Conversion:**
- `orbital_elements_to_cartesian()` applies rotation
- `cartesian_to_orbital_elements()` computes ω, i, Ω from angular momentum and eccentricity vectors

## Orbit Types

### Elliptical Orbits (e < 1)
- **Energy:** ε < 0 (negative specific orbital energy)
- **Semi-major axis:** a = -μ/(2ε) (positive)
- **Motion:** Bound, periodic
- **Propagation:** Uses elliptical Kepler equation: E - e·sin(E) = M
- **Examples:** Planets, moons, satellites

### Parabolic Orbits (|e-1| < PARABOLIC_TOLERANCE)
- **Energy:** ε ≈ 0 (zero specific orbital energy, escape velocity)
- **Parameter:** Semi-latus rectum p = h²/μ
- **Motion:** Escape trajectory, asymptotic velocity → 0 at infinity
- **Propagation:** Uses Barker's equation: D + D³/3 = M where D = tan(ν/2)
- **Examples:** Escape missions, comet-like trajectories

### Hyperbolic Orbits (e > 1)
- **Energy:** ε > 0 (positive specific orbital energy)
- **Semi-major axis:** a = -μ/(2ε) (negative)
- **Motion:** Unbound, excess velocity at infinity
- **Propagation:** Uses hyperbolic Kepler equation: H - e·sinh(H) = M
- **Examples:** Interplanetary trajectories, gravity assists

**Velocity Notes:**
- Elliptical: v < v_escape at periapsis
- Parabolic: v = v_escape at all points
- Hyperbolic: v > v_escape at all points
- All satisfy vis-viva equation: v² = μ(2/r - 1/a)

## Data Flow

### Initialization Sequence
The simulation initializes in this sequence: create_simulation() is called first, then load_system_config() parses the TOML file and loads bodies, spacecraft, and maneuvers. Next, run_all_config_validations() performs system-level validation. Then initialize_orbital_objects() converts orbital elements to local position/velocity for all bodies and spacecraft, computes global coordinates, and calculates SOI radii. Finally, the main simulation loop begins.

### Main Simulation Loop
The main simulation loop executes in this order: update_bodies_physics(), compute_global_coordinates(), update_spacecraft_physics(), execute_pending_maneuvers(), compute_spacecraft_globals(), then increments simulation time. Within update_bodies_physics(), for each body: check SOI via find_dominant_body, handle transitions by computing global coordinates from old parent, updating parent_index, computing new local coordinates, and reconstructing orbital elements. Then check velocity deviation with 1e-6 tolerance and reconstruct elements if needed. Propagate elements via propagate_orbital_elements() and update local position/velocity. compute_global_coordinates() updates all body global positions from parent.global + local. update_spacecraft_physics() performs the same velocity deviation check and propagation for spacecraft. execute_pending_maneuvers() checks each unexecuted maneuver for time or true anomaly triggers; if triggered, applies burn to local_velocity, reconstructs orbital elements, and marks executed. compute_spacecraft_globals() updates all spacecraft global positions.

### SOI Mechanics
SOI transitions are detected by calling find_dominant_body() before each physics update. If the parent changes, the body's global coordinates are computed in the old frame, the parent_index is updated, new local coordinates are computed, and orbital elements are reconstructed. Propagation then uses the new local frame.

## Location Hints

### Constants
- Gravitational constant: `G` in physics.h
- Parabolic tolerance: `PARABOLIC_TOLERANCE` in orbital_mechanics.h
- Kepler solver tolerance: `KEPLER_TOLERANCE` in orbital_mechanics.cpp
- SOI mass ratio: `MIN_MASS_RATIO` in config_validator.h

### Config Files
- Example: `configs/solar_system.toml`
- Format: TOML with [[bodies]], [[spacecraft]], [[maneuvers]] arrays

### Key Function Line Numbers
- `propagate_orbital_elements()`: orbital_mechanics.cpp:301-375
- `find_dominant_body()`: simulation.cpp:105-148
- `cartesian_to_orbital_elements()`: orbital_mechanics.cpp:186-299
- `orbital_elements_to_cartesian()`: orbital_mechanics.cpp:6-44
- `execute_maneuver()`: maneuver.cpp:166-176
- Velocity deviation check: simulation.cpp:267-270 (bodies), 291-294 (spacecraft)

### Build System
**Makefile targets:**
- `make` or `make all` - Build simulation
- `make raylib` - Build raylib dependency
- `make run` - Build and run simulation
- `make test` - Build and run automated tests
- `make test-build` - Build test executable only
- `make clean` - Remove build artifacts
- `make rebuild` - Clean and rebuild

**Dependencies:**
- g++ (C++14 standard)
- raylib (git submodule in ext/raylib/)
- raygui (header-only, in ext/raygui/)
- tomlc17 (C TOML parser, in ext/tomlc17/)
- Catch2 (test framework, for test builds only)

**Build output:**
- `orbit_sim` - Main simulation executable
- `orbit_test` - Test suite executable
- Object files in `build/` directory

**Testing:**
```bash
# Run all tests
make test

# Run specific test config
./orbit_test '[config_name]'

# Test with extra debug output (shows INFO messages)
./orbit_test -s '[config_name]'
```

**Test framework:**
- Catch2 for testing
- Use `WithinAbs(expected, tolerance)` for floating-point comparisons (NOT `Approx()`)
- Required header: `<catch2/matchers/catch_matchers_floating_point.hpp>`