vibe coding an orbital mechanics simulation to try out claude code
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

22 KiB

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

struct Vec3 {
    double x, y, z;
};

Mat3 (physics.h)

3x3 row-major matrix for coordinate transformations

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

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

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

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

struct Maneuver {
    char name[64];
    int craft_index;
    BurnDirection direction;
    double delta_v;
    TriggerType trigger_type;
    double trigger_value;
    double scheduled_dt;      // time to propagate before burn (for exact position)
    bool executed;
    double executed_time;
};

SimulationState (simulation.h)

Top-level simulation container

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:219)

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

Exact Position Burn Execution

Purpose: True anomaly triggers must execute burns at the exact orbital position, not at the current position when crossing is detected.

Mechanism:

  1. check_maneuver_trigger() calculates scheduled_dt (time to reach target anomaly)
  2. If crossing will occur within current frame (scheduled_dt < sim->dt), trigger fires
  3. execute_pending_maneuvers() propagates spacecraft by scheduled_dt to exact position
  4. Burn executes at precise orbital location
  5. Remaining frame time (sim->dt - scheduled_dt) is propagated after burn
  6. Spacecraft marked as handled to skip redundant propagation in update_spacecraft_physics()

Wraparound Handling:

  • Special case for 2π→0 crossing at periapsis
  • When current_nu > 5.0 and future_nu < 1.0, wraparound crossing is detected
  • Prevents false "moving away" rejection near angle boundaries

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(), execute_pending_maneuvers(), update_spacecraft_physics(), 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. execute_pending_maneuvers() checks each unexecuted maneuver for time or true anomaly triggers. For true anomaly triggers: if crossing detected, sets scheduled_dt to time needed to reach target. When triggered, propagates spacecraft by scheduled_dt to exact position, executes burn, propagates remaining frame time, and marks spacecraft as handled to skip in update_spacecraft_physics(). update_spacecraft_physics() propagates spacecraft not already handled this frame. 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

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:219
  • 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:

# 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>