21 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;
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:
Gin physics.h
Vector Functions:
vec3_add,vec3_sub,vec3_cross,vec3_scalevec3_magnitude,vec3_distance,vec3_normalize,vec3_dot
Matrix Functions:
mat3_identity,mat3_multiply,mat3_multiply_vec3mat3_rotation_x,mat3_rotation_zmat3_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_TOLERANCEin orbital_mechanics.h - Kepler solver tolerance:
KEPLER_TOLERANCEin 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 vectorBURN_RETROGRADE- Opposite velocityBURN_NORMAL- Along angular momentum (orbit normal)BURN_ANTINORMAL- Opposite angular momentumBURN_RADIAL_IN- Toward parentBURN_RADIAL_OUT- Away from parentBURN_CUSTOM- User-specified vector
Trigger Types:
TRIGGER_TIME- Execute at simulation time >= trigger_valueTRIGGER_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.mdfor 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.mdfor 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_completedflag,time_at_completionrecord
Critical Implementation Details
Analytical Propagation Algorithm
Function: propagate_orbital_elements(elements, dt, parent_mass) (orbital_mechanics.cpp:301-375)
Elliptical Orbits (e < 1):
- Compute mean motion: n = √(μ/a³)
- Convert true anomaly to eccentric anomaly: E = 2·atan(√((1-e)/(1+e)) · tan(ν/2))
- Compute mean anomaly: M = E - e·sin(E)
- Advance mean anomaly: M_new = M + n·dt
- Solve Kepler equation for E_new using Newton-Raphson (max 50 iterations, tolerance 1e-10)
- Convert back to true anomaly: ν_new = 2·atan(√((1+e)/(1-e)) · tan(E_new/2))
- Return elements with updated true_anomaly
Parabolic Orbits (|e-1| < PARABOLIC_TOLERANCE):
- Compute variable: D = tan(ν/2)
- Compute Barker's mean anomaly: M = D + D³/3
- Compute mean motion: n = √(μ/p³) where p = semi_latus_rectum
- Advance: M_new = M + n·dt
- Solve Barker's equation: D_new from cubic formula
- Convert: ν_new = 2·atan(D_new)
- Return elements with updated true_anomaly
Hyperbolic Orbits (e > 1):
- Compute mean motion: n = √(μ/(-a)³) (note negative a)
- Convert true anomaly to hyperbolic anomaly: H = 2·atanh(√((e-1)/(e+1)) · tan(ν/2))
- Compute mean anomaly: M = e·sinh(H) - H
- Advance: M_new = M + n·dt
- Solve hyperbolic Kepler equation for H_new using Newton-Raphson (max 50 iterations, tolerance 1e-10)
- Convert back: ν_new = 2·atan(√((e+1)/(e-1)) · tanh(H_new/2))
- 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:
-
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)
-
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:
- 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:
- Apply impulsive burn to
craft->local_velocity - Reconstruct elements:
cartesian_to_orbital_elements(local_pos, local_vel, parent_mass) - 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:
- Rotate by argument_of_periapsis (ω) around Z-axis
- Rotate by inclination (i) around X-axis
- 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 rotationcartesian_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:
Gin physics.h - Parabolic tolerance:
PARABOLIC_TOLERANCEin orbital_mechanics.h - Kepler solver tolerance:
KEPLER_TOLERANCEin orbital_mechanics.cpp - SOI mass ratio:
MIN_MASS_RATIOin 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-375find_dominant_body(): simulation.cpp:105-148cartesian_to_orbital_elements(): orbital_mechanics.cpp:186-299orbital_elements_to_cartesian(): orbital_mechanics.cpp:6-44execute_maneuver(): maneuver.cpp:166-176- Velocity deviation check: simulation.cpp:267-270 (bodies), 291-294 (spacecraft)
Build System
Makefile targets:
makeormake all- Build simulationmake raylib- Build raylib dependencymake run- Build and run simulationmake test- Build and run automated testsmake test-build- Build test executable onlymake clean- Remove build artifactsmake 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 executableorbit_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 (NOTApprox()) - Required header:
<catch2/matchers/catch_matchers_floating_point.hpp>