Compare commits

..

No commits in common. 'eade44a245409a337ccdf0bfdff8f67c0d5aff10' and '0548d864c46dff9e1e8ac3a2587aa7f338c387aa' have entirely different histories.

  1. 3
      .gitignore
  2. 14
      AGENTS.md
  3. 14
      Makefile
  4. 24
      docs/TODO
  5. 256
      docs/planning/hohmann-rendezvous-quantization-fix.md
  6. 213
      docs/planning/propagation_refactor.md
  7. 366
      docs/planning/rendezvous_planning.md
  8. 82
      docs/session_summaries/2026-03-28-ui-fixes.md
  9. 782
      docs/technical_reference.md
  10. 204
      scripts/compute_rendezvous_params.py
  11. 72
      scripts/generate_interface_summaries.sh
  12. 38
      scripts/prompt.md
  13. 431
      scripts/simulate_rendezvous.py
  14. 732
      scripts/source_documentation_analysis.md
  15. 2
      src/config_loader.h
  16. 2
      src/config_validator.h
  17. 84
      src/maneuver.cpp
  18. 2
      src/maneuver.h
  19. 51
      src/orbital_mechanics.cpp
  20. 49
      src/orbital_objects.h
  21. 2
      src/renderer.h
  22. 206
      src/rendezvous.cpp
  23. 95
      src/rendezvous.h
  24. 69
      src/simulation.cpp
  25. 30
      src/simulation.h
  26. 5
      src/spacecraft.cpp
  27. 24
      src/spacecraft.h
  28. 52
      src/test_utilities.cpp
  29. 6
      src/test_utilities.h
  30. 16
      src/ui_renderer.cpp
  31. 2
      tests/test_cartesian_to_elements_advanced.cpp
  32. 2
      tests/test_hybrid_burns.cpp
  33. 2
      tests/test_hybrid_energy_conservation.cpp
  34. 9
      tests/test_maneuver_planning.cpp
  35. 2
      tests/test_maneuvers.cpp
  36. 13
      tests/test_omega_debug.cpp
  37. 2
      tests/test_orbit_rendering.cpp
  38. 2
      tests/test_periapsis_burn.cpp
  39. 613
      tests/test_rendezvous.cpp
  40. 70
      tests/test_rendezvous.toml

3
.gitignore vendored

@ -5,9 +5,6 @@ orbit_sim
orbit_test orbit_test
tests/informational/test_time_step_stability tests/informational/test_time_step_stability
# LLM summaries
src/*.summary.md
# Session notes and conversation logs # Session notes and conversation logs
docs/session_logs docs/session_logs

14
AGENTS.md

@ -43,15 +43,10 @@
## Common Commands ## Common Commands
- Build: make - Build: make
- Test (build + run all): make test - Test: make test
- Test (rebuild only): make test-build - Test (with extra debug info): ./orbit_test -s '[config_name]'
- Test (specific tag): ./build/orbit_test '[tag_name]'
- Test (with debug output): ./build/orbit_test -s '[tag_name]'
- needed to display 'INFO' statements for successful tests with Catch2 framework - needed to display 'INFO' statements for successful tests with Catch2 framework
- See docs/technical_reference.md for full build target reference - See README.md for full build instructions
## Simulation
The simulation binary is built at `./build/orbit_sim`. Agents should NOT attempt to run the graphical simulation. Focus on testing and code changes only.
## Testing Guidelines ## Testing Guidelines
- Always use `WithinAbs()` for floating-point comparisons - Always use `WithinAbs()` for floating-point comparisons
@ -59,6 +54,3 @@ The simulation binary is built at `./build/orbit_sim`. Agents should NOT attempt
- Required header: `<catch2/matchers/catch_matchers_floating_point.hpp>` - Required header: `<catch2/matchers/catch_matchers_floating_point.hpp>`
- Usage: `REQUIRE_THAT(value, WithinAbs(expected, absolute_margin))` - Usage: `REQUIRE_THAT(value, WithinAbs(expected, absolute_margin))`
## Documentation
- Always read in docs/technical_reference.md when starting a new session
- Always read in src/*.summary.md when starting a new session

14
Makefile

@ -10,11 +10,11 @@ LDFLAGS = -L./ext/raylib/src -lraylib -lm -lpthread -ldl -lrt -lX11
# Directories # Directories
SRC_DIR = src SRC_DIR = src
BUILD_DIR ?= build BUILD_DIR = build
RAYLIB_DIR = ext/raylib/src RAYLIB_DIR = ext/raylib/src
TARGET = ${BUILD_DIR}/orbit_sim TARGET = orbit_sim
TEST_DIR = tests TEST_DIR = tests
TEST_TARGET = ${BUILD_DIR}/orbit_test TEST_TARGET = orbit_test
# Source files (exclude old mission planning files) # Source files (exclude old mission planning files)
CPP_SOURCES := $(wildcard $(SRC_DIR)/*.cpp) CPP_SOURCES := $(wildcard $(SRC_DIR)/*.cpp)
@ -64,6 +64,10 @@ clean-all: clean
# Rebuild # Rebuild
rebuild: clean all rebuild: clean all
# Run the simulation
run: $(TARGET)
./$(TARGET)
$(TEST_OBJECTS): $(BUILD_DIR)/%.o : $(TEST_DIR)/%.cpp $(TEST_OBJECTS): $(BUILD_DIR)/%.o : $(TEST_DIR)/%.cpp
$(CXX) $(CXXFLAGS) -c $< -o $@ $(CXX) $(CXXFLAGS) -c $< -o $@
@ -77,11 +81,11 @@ test-build: $(BUILD_DIR) $(C_OBJECTS) $(CPP_OBJECTS) $(TEST_OBJECTS)
build/config_loader.o \ build/config_loader.o \
build/config_validator.o \ build/config_validator.o \
build/maneuver.o \ build/maneuver.o \
build/rendezvous.o \ build/spacecraft.o \
-o $(TEST_TARGET) -lCatch2Main -lCatch2 -lm -o $(TEST_TARGET) -lCatch2Main -lCatch2 -lm
# Run automated test suite # Run automated test suite
test: test-build test: test-build
./$(TEST_TARGET) ./$(TEST_TARGET)
.PHONY: all clean clean-all rebuild test test-build raylib ctags .PHONY: all clean clean-all rebuild run test test-build raylib ctags

24
docs/TODO

@ -6,24 +6,12 @@ If you see modifications to this file in git status, IGNORE them and do not comm
=== next steps === === next steps ===
- interplanetary transfers - debug maneuver at periapsis
- SOI boundary testing
- draw SOI boundry in graphical sim
- reorganize project into static library. easier for testing/llm doc generation,
and will be easier to try out new graphics/UI libraries
- can then better work on code cleaning issues, FIXMEs etc
- remove tests/informational/*
- remove RK4 integration implementation?
- refactor tests and check logic for edge cases
- add reset/load new config UI control - add reset/load new config UI control
- UI fixes - UI fixes
- SOI boundary testing
- add more documentation for the LLM to use 'make test-build' and run individual tests - draw SOI boundry in graphical sim
- write out/fix inline code docs with the documentation generation script, and review - interplanetary transfers
- can then work on a better method for commenting/detecting changes via script
- move summary docs to separate folder and combine into a big 'API reference' doc
=== Project Documentation === === Project Documentation ===
@ -37,13 +25,13 @@ If you see modifications to this file in git status, IGNORE them and do not comm
- actually, we should start out simpler. create and test maneuvers in planetary - actually, we should start out simpler. create and test maneuvers in planetary
frame before starting on interplanetary patched conics, oops frame before starting on interplanetary patched conics, oops
- will still need hohmann transfers, circularizing, plane change, etc. - will still need hohmann transfers, circularizing, plane change, etc.
- something is not working with hohmann transfers, maybe SOI/parent logic?
- orbital rendezveuz
- interplanetary/SOI changing transfers - interplanetary/SOI changing transfers
- remove bodies on non closed orbits after they are a certain distance from - remove bodies on non closed orbits after they are a certain distance from
the root body, or bary-center the root body, or bary-center
- should try to use the same code path for bodies and craft on 'update' functions - should try to use the same code path for bodies and craft on 'update' functions
as much as possible as much as possible
- combine code paths for trigger types in check_maneuver_trigger().
- see TODO in src/maneuver.cpp
=== Simulation Config Files === === Simulation Config Files ===
- a format with nested children would be nice - a format with nested children would be nice

256
docs/planning/hohmann-rendezvous-quantization-fix.md

@ -1,256 +0,0 @@
# Hohmann Transfer Rendezvous - Burn Timing Quantization Analysis
## Session Date
2026-04-19
## Problem Statement
The Hohmann transfer rendezvous simulation was failing with ~1.3 km separation despite correct phasing calculations. Investigation revealed two issues:
1. **Burn timing quantization**: Time-triggered burns execute at step boundaries, not exact trigger times
2. **cartesian_to_orbital_elements bug**: Coplanar orbit omega calculation was incorrect
## Root Cause Analysis
### 1. Burn Timing Quantization
**How it works:**
- `check_maneuver_trigger()` for `TRIGGER_TIME` uses simple comparison: `sim->time >= maneuver->trigger_value`
- When triggered, the burn executes at the **first step where `sim->time >= trigger_value`**
- No sub-step interpolation for time triggers
**Verified behavior** (from `test_maneuver_planning.cpp`):
| Trigger Time | Step Boundary | Actual Execution | Delay |
|-------------|---------------|-----------------|-------|
| t=305.0 | t=310.0 | t=310.0 | **5.0s** |
| t=300.0 | t=300.0 | t=300.0 | **0.0s** |
| t=62807.0 | t=62810.0 | t=62810.0 | **3.0s** |
**Impact on Hohmann transfer:**
- Arrival trigger: t=62804.47 (calculated precisely)
- Step boundaries: ..., 62800, 62810, ...
- Actual execution: t=62810 (5.53s late)
- Position drift: ~5.53s × ~7672 m/s ≈ **42 km** of orbital travel
### 2. cartesian_to_orbital_elements Bug
**Location:** `src/orbital_mechanics.cpp`, lines 300-320
**Bug:** For coplanar orbits (inclination < 0.01 rad), the function was setting `omega = 0.0` instead of computing the longitude of periapsis.
**Fix:**
```cpp
} else if (e > 1e-10) {
// Coplanar or near-circular: use longitude of periapsis
omega = atan2(e_vec.y, e_vec.x);
if (omega < 0.0) {
omega += 2.0 * M_PI;
}
} else {
omega = 0.0;
}
```
**Impact:** Without this fix, Hohmann separation goes from 8.75m → 3.22 million meters.
## DT Reduction Results
| TIME_STEP | Separation | Test Result |
|-----------|-----------|-------------|
| 10.0 s | 1,324 m | ❌ Failed (>100m) |
| 5.0 s | 458 m | ❌ Failed (>100m) |
| 2.0 s | 228 m | ❌ Failed (>100m) |
| 1.0 s | 55 m | ✅ Passed |
| 0.5 s | 69 m | ❌ Failed (>100m) |
| 0.1 s | 8.75 m | ✅ Passed |
**Note:** Values measured from `test_maneuver_planning.cpp` DT sweep test (2026-04-20).
The 0.5s value (69m) is higher than the original estimate (40m) due to the specific
trigger offset within the step boundary for this scenario. The 2.0s value (228m) is
also higher than the original estimate (150m). The 10.0s value (1324m) matches the
original measurement exactly.
**Key insight:** DT reduction dramatically improves accuracy:
- 24x improvement from 10s→1s
- 6x more from 1s→0.1s
## Test Results Summary
| Test Category | Before Fix | After Fix | Status |
|--------------|-----------|-----------|--------|
| rendezvous (8 cases) | 87 passed | **107 passed** | ✅ All pass |
| maneuver_planning (6 cases) | 3 passed | **6 passed** | ✅ All pass |
| omega (2 cases) | 1 failed | **2 passed** | ✅ All pass |
| **Total** | 156/160 pass | **154/154 pass** | All pass |
**Notes:**
- The old `rendezvous` (CW guidance) module was removed entirely, eliminating 3 pre-existing test failures
- `test_maneuver_timing.cpp` was merged into `test_maneuver_planning.cpp`
- `rendezvous_hohmann` was renamed to `rendezvous` (CW module removed, only Hohmann remains)
- All 154 remaining test cases pass (240,445 assertions)
## Current Code Path (After 2026-04-26 Sub-step Interpolation)
The maneuver trigger check and execution are merged into `update_spacecraft_physics()`. Both trigger types now use sub-step interpolation:
```cpp
// In update_spacecraft_physics(), per spacecraft:
check_maneuver_trigger(maneuver, craft, sim);
// → For TRIGGER_TIME: computes dt_to_burn = trigger_value - sim->time
// → For TRIGGER_TRUE_ANOMALY: computes dt_needed from mean anomaly delta
// → Both set scheduled_dt = dt_to_burn (0 to sim->dt)
if (maneuver_fired) {
craft->orbit = propagate_orbital_elements(craft->orbit, burn_dt, ...);
execute_maneuver(fired_maneuver, ...);
craft->orbit = propagate_orbital_elements(craft->orbit, remaining_dt, ...);
} else {
craft->orbit = propagate_orbital_elements(craft->orbit, sim->dt, ...);
}
```
**For `TRIGGER_TIME`:** `burn_dt` is the exact sub-step offset from step start to trigger time. The spacecraft propagates to the precise trigger position before the burn, then continues for `sim->dt - burn_dt`.
**For `TRIGGER_TRUE_ANOMALY`:** `burn_dt` is set to exact seconds-to-target by `check_maneuver_trigger()`, so the spacecraft propagates to the exact burn position before the burn executes.
**Edge case:** If `sim->time > trigger_value` (trigger passed in a previous step), `scheduled_dt` is clamped to 0 and the burn fires immediately at the current position.
## Burn Timing Quantization — RESOLVED (2026-04-26)
**Option A (Sub-step Interpolation)** is now implemented for both `TRIGGER_TIME` and `TRIGGER_TRUE_ANOMALY`.
**Implementation:**
- `check_maneuver_trigger()` computes `dt_to_burn = trigger_value - sim->time` for time triggers
- `update_spacecraft_physics()` propagates to exact burn time, executes burn, propagates remainder
- Quantization error is eliminated: burns execute at the precise trigger time
**Edge case:** When `sim->time > trigger_value` (trigger passed in a previous step), `scheduled_dt` is clamped to 0 and the burn fires immediately at the current position.
### Remaining Options (No Longer Needed)
### Option B: Snap Trigger Times to Step Boundaries
**Approach:** In `calculate_next_hohmann_wait_time()`, snap the calculated wait time to the nearest step boundary.
**Changes needed:**
1. In `calculate_next_hohmann_wait_time()`:
- After calculating wait time, snap to step boundary: `wait_time = ceil(wait_time / DT) * DT`
- This ensures the trigger aligns with a simulation step
**Pros:** Simple, minimal code changes
**Cons:** Introduces systematic timing error, may affect phasing accuracy
### Option C: Accept Quantization Error
**Approach:** Keep current behavior but set realistic thresholds based on DT.
**Changes needed:**
1. Calculate expected quantization error: `max_error = DT`
2. Set rendezvous threshold proportional to DT: `threshold = 100 * DT` (meters)
3. Document the limitation
**Pros:** Simplest, no code changes
**Cons:** Less accurate, threshold depends on DT choice
## Strategy for Testing with Larger Time Steps
### Goal
Understand the accuracy limitations of the simulation at realistic DT values (10s, 30s) to set appropriate rendezvous thresholds.
### Test Plan
#### Phase 1: Baseline at Current DT (0.1s)
- ✅ Already done: 8.75m separation at DT=0.1s
#### Phase 2: Systematic DT Sweep
Run the same Hohmann transfer test at increasing DT values:
| DT | Expected Steps | Expected Separation |
|----|---------------|-------------------|
| 0.1s | ~628,000 | ~8.75 m |
| 0.5s | ~125,600 | ~40 m (estimate) |
| 1.0s | ~62,800 | ~55 m |
| 2.0s | ~31,400 | ~100-200 m (estimate) |
| 5.0s | ~12,560 | ~500 m (estimate) |
| 10.0s | ~6,280 | ~1,324 m |
| 30.0s | ~2,093 | ~4,000 m (estimate) |
**Method:**
1. Create a new test file `tests/test_hohmann_dt_sweep.cpp`
2. Run the same Hohmann transfer scenario at each DT value
3. Record: final separation, radius error, relative velocity
4. Plot separation vs DT to determine the relationship
#### Phase 3: Quantization Impact Analysis
Test the effect of burn timing quantization specifically:
| Scenario | Trigger Offset | Expected Delay |
|----------|---------------|----------------|
| Exact boundary | 0s | 0s |
| 5s after boundary | 5s | 5s |
| 9s after boundary | 9s | 1s |
**Method:**
1. For each DT, run the Hohmann transfer multiple times with different trigger offsets
2. Measure the variation in final separation
3. Determine if quantization error dominates over integration error
#### Phase 4: Threshold Recommendation
Based on Phase 2 & 3 results, recommend:
- Maximum DT for rendezvous operations
- Separation threshold as a function of DT
- Whether sub-step interpolation is necessary
### Implementation Notes
- Use `calculate_next_hohmann_wait_time()` with `min_wait_time` to control trigger timing
- Keep all other parameters constant (initial conditions, maneuver DVs, etc.)
- Use `WithinAbs()` with increasing margins to find the threshold that passes at each DT
## Completed Work
### Files Modified
- `src/orbital_mechanics.cpp` - Fixed coplanar orbit omega calculation
- `src/rendezvous.cpp` (renamed from `rendezvous_hohmann.cpp`) - Added 3 new functions (validate, relative period, next wait time)
- `src/rendezvous.h` (renamed from `rendezvous_hohmann.h`) - Added function declarations
- `src/test_utilities.cpp` - Added `dump_simulation_state()` helper
- `src/test_utilities.h` - Added function declaration
- `tests/test_rendezvous.cpp` (renamed from `test_rendezvous_hohmann.cpp`) - Updated integration test with DT=0.1
- `tests/test_rendezvous.toml` (renamed from `test_rendezvous_hohmann.toml`) - Reverted to original values
- `tests/test_maneuver_planning.cpp` - Added 3 burn timing quantization tests (merged from test_maneuver_timing.cpp), plus 3 DT sweep tests (2026-04-20)
- `tests/test_omega_debug.cpp` - Updated to accept new coplanar omega behavior
- `Makefile` - Updated object file references
- `src/simulation.cpp` - Merged `execute_pending_maneuvers()` into `update_spacecraft_physics()` (2026-04-20)
- `src/simulation.h` - Removed `execute_pending_maneuvers()` declaration (2026-04-20)
### Files Removed
- `src/rendezvous.h` (old CW module) - replaced by Hohmann-only rendezvous.h
- `src/rendezvous.cpp` (old CW module) - replaced by Hohmann-only rendezvous.cpp
- `tests/test_rendezvous.cpp` (old CW tests) - replaced by Hohmann-only test_rendezvous.cpp
- `tests/test_rendezvous.toml` (old CW config) - replaced by Hohmann-only test_rendezvous.toml
## Remaining Work
### DT Sweep Tests (COMPLETED - 2026-04-20)
Measured in `test_maneuver_planning.cpp` via `DT sweep: Hohmann transfer separation vs time step`:
| DT | Expected Steps | Measured Separation |
|----|---------------|-------------------|
| 0.1s | ~628,000 | 8.75 m |
| 0.5s | ~125,600 | 69 m |
| 1.0s | ~62,800 | 55 m |
| 2.0s | ~31,400 | 228 m |
| 5.0s | ~12,560 | 458 m |
| 10.0s | ~6,280 | 1,324 m |
The separation scales roughly linearly with DT for larger values, consistent with
quantization error being the dominant factor (position error ≈ timing_error × velocity,
where timing_error is uniformly distributed in [0, DT]).
Additional tests:
- `DT sweep: quantization error is bounded by DT` — verifies error is always in [0, DT)
- `DT sweep: Hohmann arrival burn timing error` — measures exact timing error at each DT
### Threshold Recommendation
With sub-step interpolation implemented, burn timing quantization is eliminated. DT sweep results now reflect integration error rather than quantization error. Recommend:
- Maximum DT for rendezvous operations based on integration accuracy
- Separation threshold set by orbital dynamics, not quantization bounds
- Sub-step interpolation is now active for both trigger types

213
docs/planning/propagation_refactor.md

@ -1,213 +0,0 @@
# Propagation Call Chain Analysis
## Session Date
2026-04-20
## Objective
Audit all call sites of `propagate_orbital_elements()` for spacecraft and trace the `update_simulation()` call chain through `execute_pending_maneuvers()` and `update_spacecraft_physics()` to identify inefficiencies and confusing branching.
## Recent Changes
### 2026-04-20: Eliminated Redundant Propagation in `check_maneuver_trigger()`
Replaced the per-frame propagation probe in `check_maneuver_trigger()` (true-anomaly branch) with a direct analytical calculation using mean anomaly delta.
**What changed:**
- Removed `propagate_orbital_elements()` call from `check_maneuver_trigger()` — eliminates ~24k redundant Kepler solves per Hohmann transfer wait period
- Removed dead `angle_between()` static helper
- Replaced the `angle_between`, `future_diff`, and `wraparound_crossing` checks with a cleaner `dt_needed <= 0.0` guard
- Added TODO comment noting elliptical-orbits-only limitation
**Impact:**
- Same correctness, fewer function calls, simpler logic
- All 154 tests pass (240,445 assertions)
- For a Hohmann transfer with ~244,000s wait time and DT=10s: ~24,400 fewer `propagate_orbital_elements()` calls
**Remaining:** See Issue 1 below — still has a TODO for parabolic/hyperbolic orbit support.
### 2026-04-20: Merged Maneuver Execution into Single Propagation Pass
Consolidated `execute_pending_maneuvers()` into `update_spacecraft_physics()` so every spacecraft goes through exactly one propagation path.
**What changed:**
- Removed `spacecraft_handled_this_frame[256]` static array (magic number)
- Removed `reset_spacecraft_tracking()` function
- Removed `execute_pending_maneuvers()` function
- Check maneuver triggers before propagation in `update_spacecraft_physics`
- Propagate to burn time, execute burn, propagate remainder (same order as old code)
- Removed `execute_pending_maneuvers` declaration from `simulation.h`
**Impact:**
- Single propagation call per spacecraft per frame (was 2 separate code paths)
- No static array with hardcoded size
- Simpler call chain: no separate maneuver pass
- Enables sub-step interpolation for true-anomaly triggers
- Sets up foundation for interpolated time triggers
**Remaining:** Adds TODO in `update_spacecraft_physics` about testing interpolated burns.
---
## Call Chain: `update_simulation()` (Current)
```
update_simulation()
├── update_bodies_physics()
│ └── for each body:
│ ├── find_dominant_body()
│ ├── orbital_elements_to_cartesian() → expected_vel
│ ├── velocity drift check (every frame)
│ ├── cartesian_to_orbital_elements() (if drift > 1e-6)
│ └── propagate_orbital_elements()
├── compute_global_coordinates()
├── update_spacecraft_physics()
│ └── for each craft:
│ ├── orbital_elements_to_cartesian() → expected_vel
│ ├── velocity drift check
│ ├── cartesian_to_orbital_elements() (if drift)
│ ├── check_maneuver_trigger() ← inline, before propagation
│ │ ├── TRIGGER_TIME: sim->time >= trigger_value
│ │ └── TRIGGER_TRUE_ANOMALY:
│ │ └── analytical mean anomaly delta (no propagation)
│ └── if maneuver fired:
│ │ ├── propagate_orbital_elements(burn_dt)
│ │ ├── execute_maneuver()
│ │ │ ├── apply_impulsive_burn()
│ │ │ └── cartesian_to_orbital_elements()
│ │ └── propagate_orbital_elements(remaining_dt)
│ └── else:
│ └── propagate_orbital_elements(sim->dt)
├── compute_spacecraft_globals()
└── sim->time += sim->dt
```
## Issues Found
### Issue 1: Redundant Propagation for True-Anomaly Triggers
**Status:** RESOLVED — eliminated in favor of analytical mean anomaly delta calculation.
**Location:** `src/maneuver.cpp`, `check_maneuver_trigger()` (true-anomaly branch)
Previously, for every frame that a true-anomaly maneuver was pending, `propagate_orbital_elements()` was called as a "look-ahead" probe to determine if the target angle was approaching.
```cpp
// OLD (maneuver.cpp:147) — REMOVED
OrbitalElements future_elements = propagate_orbital_elements(craft->orbit, sim->dt, parent->mass);
```
**Resolution:** Replaced with direct analytical computation of `dt_needed` from mean anomaly delta. The analytical solution is more precise (no discretization error) and eliminates ~24k redundant Kepler solves per Hohmann transfer wait period at DT=10s.
**Remaining:** The current implementation only handles elliptical orbits. TODO comment added for parabolic (Barker's equation) and hyperbolic branches.
### Issue 2: Mixed Concerns in `execute_pending_maneuvers()`
**Status:** RESOLVED — merged into `update_spacecraft_physics()`.
The function previously performed two distinct responsibilities:
1. **Checking** trigger conditions (calls `check_maneuver_trigger()`)
2. **Executing** the burn (propagation → burn → propagation)
**Resolution:** Trigger checking now happens inline in `update_spacecraft_physics()` before the propagation decision. The code is cleaner: check → decide propagation amount → propagate → burn (if needed) → propagate remainder.
### Issue 3: Ambiguous `scheduled_dt` Semantics
**Status:** RESOLVED — both trigger types now set `scheduled_dt` to the exact sub-step offset.
| Trigger Type | `scheduled_dt` meaning | Set by |
|-------------|----------------------|--------|
| `TRIGGER_TIME` | `trigger_value - sim->time` (step start) | `check_maneuver_trigger()` |
| `TRIGGER_TRUE_ANOMALY` | Seconds until exact burn position | `check_maneuver_trigger()` |
In `update_spacecraft_physics()`:
```cpp
if (maneuver_fired) {
craft->orbit = propagate_orbital_elements(craft->orbit, burn_dt, ...);
execute_maneuver(fired_maneuver, ...);
craft->orbit = propagate_orbital_elements(craft->orbit, remaining_dt, ...);
}
```
Both trigger types now produce the same `scheduled_dt` semantics: seconds from step start until the burn executes. The only difference is how that value is computed (simple subtraction vs. mean anomaly propagation).
### Issue 4: Time-Triggered Burns Propagate from Wrong State
**Status:** RESOLVED — sub-step interpolation implemented (2026-04-26).
For time triggers, the sequence is now:
```
Frame N: sim->time = 310.0 (trigger_value = 305.0)
check: sim->time >= trigger_value → true
dt_to_burn = trigger_value - (sim->time - sim->dt) = 305 - 300 = 5.0
burn_dt = 5.0
remaining_dt = 5.0
propagate(5.0) → craft at trigger position
execute_maneuver() → burn applies at sim->time=305
propagate(5.0) → continue simulation
```
The spacecraft propagates to the exact trigger position before the burn, then continues for the remaining timestep. This eliminates the burn timing quantization problem documented in `docs/planning/hohmann-rendezvous-quantization-fix.md`.
**Edge case:** When `sim->time > trigger_value` (trigger passed in a previous step), `scheduled_dt` is clamped to 0 and the burn fires immediately at the current position.
### Issue 5: Hardcoded Array Size
**Status:** RESOLVED — removed `spacecraft_handled_this_frame[256]` static array.
**Location:** Previously `src/simulation.cpp`, static declaration
```cpp
// REMOVED
static bool spacecraft_handled_this_frame[256];
```
**Resolution:** The array is no longer needed since maneuver checking is inline in the propagation loop. No separate tracking mechanism required.
### Issue 6: Duplicated Propagation Logic
**Status:** RESOLVED — single propagation path.
The "normal" propagation path and the maneuver path are now the same code path:
```cpp
if (maneuver_fired) {
craft->orbit = propagate_orbital_elements(craft->orbit, burn_dt, ...);
execute_maneuver(...);
craft->orbit = propagate_orbital_elements(craft->orbit, remaining_dt, ...);
} else {
craft->orbit = propagate_orbital_elements(craft->orbit, sim->dt, ...);
}
```
When `burn_dt == 0` and `remaining_dt == sim->dt`, the maneuver path is: propagate(0) → burn → propagate(sim->dt). The normal path is: propagate(sim->dt). The propagation logic is no longer duplicated.
## Call Sites Summary
### Final State (3 call sites)
| Location | Function | Context | dt value |
|----------|----------|---------|----------|
| `simulation.cpp:287` | `update_bodies_physics()` | Normal body propagation | `sim->dt` |
| `simulation.cpp:315` | `update_spacecraft_physics()` | Normal craft propagation | `sim->dt` |
| `simulation.cpp:338` | `update_spacecraft_physics()` | Pre-burn sub-step propagation | `burn_dt` (0 to `sim->dt`) |
| `simulation.cpp:343` | `update_spacecraft_physics()` | Post-burn remaining propagation | `sim->dt - burn_dt` |
Remaining: 4 call sites (1 for bodies, 3 for spacecraft), 2 distinct contexts for spacecraft (normal propagation, sub-step execution).
## Remaining Work
### Parabolic/Hyperbolic Orbit Support (TODO in code)
The analytical `dt_needed` calculation in `check_maneuver_trigger()` only handles elliptical orbits. TODO comment added for:
- Parabolic: Barker's equation `D + D³/3 = M`, where `D = tan(ν/2)`
- Hyperbolic: `e·sinh(H) - H = M`, using hyperbolic anomaly `H`
### Code Path Unification (TODO in maneuver.cpp)
A TODO was added to combine the TRIGGER_TIME and TRIGGER_TRUE_ANOMALY branches in `check_maneuver_trigger()` into a single code path. Both now compute `scheduled_dt` the same way (sub-step offset), so the logic can be unified. The main complication is maneuver array ordering: when an earlier maneuver fires and changes the orbit, cached trigger times for later maneuvers become stale and must be recomputed.

366
docs/planning/rendezvous_planning.md

@ -1,366 +0,0 @@
# Orbital Rendezvous Planning - Implementation Plan
## Overview
This document outlines the implementation plan for adding orbital rendezvous planning capabilities to the simulation.
## Architecture
- C-style C++ (structs/functions, NO classes/templates)
- Follow existing patterns in src/
- Use existing orbital mechanics functions (orbital_mechanics.h/cpp)
- Integrate with maneuver system (maneuver.h/cpp)
## Implementation Plan
### Phase 1: Core - Lambert Solver, Relative Orbit Analysis
#### 1.1 Create Rendezvous Target Structure
**Location:** `maneuver.h`
```cpp
struct RendezvousTarget {
int target_craft_index;
OrbitalElements target_elements;
Vec3 target_position;
Vec3 target_velocity;
double max_encounter_distance;
double max_relative_velocity;
};
```
**Purpose:** Hold target state for rendezvous calculations
#### 1.2 Implement Lambert's Problem Solver
**Location:** `maneuver.cpp`
```cpp
// Solve two-point boundary value problem
// Given: r1, r2, time_of_flight
// Returns: v1 (departure velocity), v2 (arrival velocity)
bool solve_lambert(Vec3 r1, Vec3 r2, double time_of_flight,
double central_mass, Vec3* v1, Vec3* v2);
```
**Algorithm:** Universal variable formulation (Gooding's method)
- Handle elliptical, parabolic, and hyperbolic transfers
- Iterative solution with convergence tolerance
- Return success/failure status
#### 1.3 Implement Relative Orbit Analysis
**Location:** `maneuver.cpp`
```cpp
// Calculate relative orbit using Clohessy-Wiltshire (Hill's) equations
// Returns relative orbital elements in LVLH frame
void calculate_relative_orbit(Vec3 r_rel, Vec3 v_rel, double mu,
Vec3* h_rel, Vec3* e_rel);
// Calculate relative position/velocity in LVLH frame
void lvih_frame_transform(Vec3 r, Vec3 v, Vec3 r_parent, Vec3 v_parent,
Vec3* r_rel, Vec3* v_rel);
```
**Purpose:** Enable precise phasing calculations
---
### Phase 2: Planning - Phasing Maneuvers, Rendezvous Planning
#### 2.1 Implement Phasing Maneuver Planning
**Location:** `maneuver.cpp`
```cpp
// Calculate phasing orbit to adjust relative position along orbit
// Returns required delta-v and phasing orbit parameters
bool calculate_phasing_maneuver(Spacecraft* craft, Spacecraft* target,
double phasing_angle, double central_mass,
double* dv, double* phasing_period);
```
**Algorithm:**
- Compute current phase difference
- Calculate semi-major axis for phasing orbit
- Determine delta-v for phasing burn
- Compute phasing orbit period
#### 2.2 Implement Rendezvous Transfer Planning
**Location:** `maneuver.cpp`
```cpp
// Calculate optimal rendezvous transfer
// Returns departure burn parameters and transfer duration
bool calculate_rendezvous_transfer(Spacecraft* craft, Spacecraft* target,
double central_mass, double max_transfer_time,
double* departure_dv, double* transfer_time,
Vec3* departure_direction);
```
**Algorithm:**
- Use Lambert solver to find transfer orbit
- Optimize for minimum delta-v
- Calculate departure burn direction and magnitude
- Determine insertion burn requirements
---
### Phase 3: Integration - Complete calculate_rendezvous(), Configuration Support
#### 3.1 Add Rendezvous Maneuver Type
**Location:** `maneuver.h`
```cpp
enum RendezvousType {
RENDEZVOUS_PROGRADE,
RENDEZVOUS_RETROGRADE,
RENDEZVOUS_PHASING,
RENDEZVOUS_CUSTOM
};
struct Maneuver {
// ... existing fields ...
RendezvousType rendezvous_type;
int target_craft_index;
double max_encounter_distance;
double max_relative_velocity;
};
```
**Purpose:** Store rendezvous-specific parameters
#### 3.2 Add New Trigger Type
**Location:** `maneuver.h`
```cpp
enum TriggerType {
TRIGGER_TIME,
TRIGGER_TRUE_ANOMALY,
TRIGGER_RENDZVOUS_COMPLETE // New trigger for rendezvous completion
};
```
#### 3.3 Extend Config Loader
**Location:** `config_loader.cpp`
```cpp
// Parse rendezvous parameters from TOML
bool load_rendezvous_config(toml::table* root, SimulationState* sim);
```
**Config Format:**
```toml
[[rendezvous]]
craft_index = 0
target_craft_index = 1
max_encounter_distance = 1000.0 # meters
max_relative_velocity = 0.1 # m/s
```
---
### Phase 4: Execution - Multi-burn Handling, Rendezvous Detection
#### 4.1 Implement Rendezvous Execution Logic
**Location:** `maneuver.cpp`
```cpp
// Execute rendezvous maneuver with multi-burn sequence
void execute_rendezvous_maneuver(Maneuver* maneuver, Spacecraft* craft,
SimulationState* sim, double current_time);
// Handle mid-course corrections during transfer
void execute_mid_course_correction(Maneuver* maneuver, Spacecraft* craft,
SimulationState* sim);
```
**Features:**
- Multi-burn sequences (departure, mid-course, insertion)
- Update target orbital elements as simulation progresses
- Track rendezvous progress
#### 4.2 Implement Rendezvous Detection
**Location:** `maneuver.cpp`
```cpp
// Check if rendezvous is complete
bool check_rendezvous_complete(Spacecraft* craft, Spacecraft* target,
RendezvousTarget* target_params);
// Calculate encounter state
void compute_encounter_state(Spacecraft* craft, Spacecraft* target,
double central_mass, double* encounter_distance,
double* relative_velocity);
```
**Checks:**
- Distance within encounter parameters
- Relative velocity within bounds
- Proper phasing achieved
---
### Phase 5: Validation - Tests, Feasibility Checks
#### 5.1 Add Rendezvous Validation Functions
**Location:** `maneuver.cpp`
```cpp
// Validate rendezvous parameters
bool validate_rendezvous_parameters(Spacecraft* craft, Spacecraft* target,
RendezvousTarget* target_params,
double central_mass);
// Check if rendezvous is feasible
bool is_rendezvous_feasible(Spacecraft* craft, Spacecraft* target,
double central_mass, double max_dv_budget);
```
**Validation:**
- Delta-v budget feasibility
- Target reachability
- Encounter geometry validity
#### 5.2 Create Supporting Utility Functions
**Location:** `maneuver.cpp`
```cpp
// Calculate optimal transfer time
double calculate_transfer_time(Vec3 r1, Vec3 r2, double central_mass);
// Compute target state at encounter time
void compute_target_state_at_time(Spacecraft* target, double encounter_time,
SimulationState* sim,
Vec3* position, Vec3* velocity);
// Calculate insertion delta-v
double calculate_insertion_dv(Vec3 v_arrival, Vec3 v_target,
double max_relative_velocity);
```
#### 5.3 Create Unit Tests
**Location:** `tests/test_rendezvous.cpp`
**Test Cases:**
1. Lambert solver accuracy (known solutions)
2. Relative orbit calculations (Clohessy-Wiltshire verification)
3. Phasing maneuver calculations (phase angle accuracy)
4. Rendezvous feasibility checks (delta-v budget)
5. End-to-end rendezvous planning workflow
6. Rendezvous detection (distance/velocity thresholds)
7. Multi-burn sequence execution
8. Mid-course correction accuracy
**Testing Guidelines:**
- Use `WithinAbs()` for floating-point comparisons (NOT `Approx()`)
- Required header: `<catch2/matchers/catch_matchers_floating_point.hpp>`
- Test with various orbital configurations (elliptical, circular, inclined)
---
## Technical Considerations
### Lambert's Problem
- **Universal Variable Formulation:** More robust than classical methods
- **Convergence:** Newton-Raphson iteration with tolerance ~1e-10
- **Time of Flight:** Key variable affecting delta-v
- **Multiple Solutions:** Handle short-way and long-way transfers
### Relative Orbit Analysis
- **Clohessy-Wiltshire Equations:** Linearized relative motion in LVLH frame
- **Validity:** Small relative distances (< 10% of orbital radius)
- **Extensions:** Non-linear corrections for large separations
### Phasing Maneuvers
- **Semi-major Axis Selection:** Determines phasing orbit period
- **Phase Angle:** Angular separation along orbit
- **Delta-v:** Varies with phase angle and orbital parameters
### Multi-Body Considerations
- **Patched Conics:** For interplanetary rendezvous (advanced)
- **SOI Transitions:** Handle frame changes during transfer
- **Third-Body Perturbations:** For high-precision requirements
### Implementation Notes
- Use existing vector/orbital element functions for consistency
- Follow ZII (Zero Is Initialization) pattern: `MyStruct s = {0};`
- Minimal comments - code should be self-documenting
- No decorative comment blocks (===, ---, etc.)
- Small, focused functions
- Follow existing patterns in src/
---
## Priority Order
1. **Critical Path:** Lambert solver, relative orbit analysis
2. **Core Functionality:** Phasing maneuvers, rendezvous planning
3. **Integration:** Complete `calculate_rendezvous()`, configuration support
4. **Execution:** Multi-burn handling, rendezvous detection
5. **Validation:** Tests, feasibility checks
---
## Dependencies
### Existing Modules
- `physics.h/cpp` - Vector/matrix math
- `orbital_mechanics.h/cpp` - Orbital element conversions, propagation
- `simulation.h/cpp` - Simulation state management
- `spacecraft.h` - Spacecraft structure
- `maneuver.h/cpp` - Maneuver system (primary integration point)
- `config_loader.h/cpp` - TOML parsing
### New Files (Optional)
- `rendezvous.h` - Rendezvous-specific structures and declarations
- `rendezvous.cpp` - Rendezvous implementation (if code becomes large)
---
## Testing Strategy
### Unit Tests
- Lambert solver with known solutions
- Relative orbit calculations vs analytical solutions
- Phasing maneuver accuracy
- Feasibility checks
### Integration Tests
- End-to-end rendezvous planning
- Multi-burn sequence execution
- Rendezvous detection accuracy
### Validation Tests
- Energy conservation during transfers
- Orbital element consistency
- Long-term stability of rendezvous orbits
---
## Future Enhancements
### Advanced Features
- **Optimal Control:** Minimize fuel consumption
- **Autonomous Navigation:** Real-time rendezvous planning
- **Docking Procedures:** Final approach and docking sequence
- **Formation Flying:** Multiple spacecraft coordination
- **Gravity Assist:** Use planetary flybys for rendezvous
### Performance Optimizations
- **Adaptive Time Stepping:** Smaller steps during critical maneuvers
- **Parallel Computation:** Multi-threaded Lambert solver
- **GPU Acceleration:** For large-scale simulations
---
## References
### Orbital Mechanics
- Curtis, H. D. "Orbital Mechanics for Engineering Students"
- Battin, R. H. "An Introduction to the Mathematics and Methods of Astrodynamics"
- Vallado, D. A. "Fundamentals of Astrodynamics and Applications"
### Lambert's Problem
- Gooding, R. H. "A Procedure for the Solution of Lambert's Problem"
- Izzo, D. "Revisiting Lambert's Problem"
### Relative Orbit
- Clohessy, W. H., & Wiltshire, R. S. "Terminal Guidance System for Satellite Rendezvous"
- Hill, G. W. "Researches in the Lunar Theory"

82
docs/session_summaries/2026-03-28-ui-fixes.md

@ -1,82 +0,0 @@
# UI Fixes Session - 2026-03-28
## Changes Made
### 1. Fixed Special Characters (line 767)
- **Issue**: Unicode checkmark characters '✓' and '✗' may not render properly on all systems
- **Fix**: Replaced with ASCII brackets: `[OK]` and `[X]`
- **Location**: `render_orbital_preview()` function
### 2. Fixed Text Overlapping in Maneuver Create Tab (line 687)
- **Issue**: Orbital preview was rendered at `current_y + 10`, overlapping with the value_text_label at `current_y`
- **Fix**: Moved orbital preview to `current_y + 50` to add proper spacing
- **Location**: `render_maneuver_create_tab()` function
### 3. Fixed Dropdown Box Positioning in Hohmann Tab (lines 843, 859, 872, 885)
- **Issue**: All dropdown boxes were positioned at `current_y - 35`, which placed them 70 pixels above the labels
- **Fix**: Changed all dropdown positions to use `current_y` correctly
- **Location**: `render_maneuver_hohmann_tab()` function
- Spacecraft dropdown
- Transfer parent dropdown
- Current body dropdown
- Target body dropdown
### 4. Fixed Text Overlapping in Edit Maneuver Tab (line 1450)
- **Issue**: Same overlapping issue as in create tab - orbital preview at `current_y + 10`
- **Fix**: Moved orbital preview to `current_y + 50`
- **Location**: `render_maneuver_edit_tab()` function
### 5. Added Reset Simulation Button (lines 103-119)
- **New Feature**: Added `[R] Reset Simulation` button in info panel
- **Functionality**: Resets simulation time to 0 when pressed
- **Location**: `render_info()` function
### 6. Added Load Config Button (lines 120-127)
- **New Feature**: Added `[L] Load Config` button in info panel
- **Functionality**: Placeholder for reloading current config file
- **Location**: `render_info()` function
## Code Quality Notes
- All changes maintain existing code style and patterns
- No new warnings introduced (only pre-existing warnings remain)
- Changes are minimal and focused on specific issues
- Follows existing UI layout conventions
## Build Results
```bash
make clean && make
```
- Build successful with only pre-existing warnings
- No new compilation errors introduced
## Testing
To test the UI fixes:
```bash
make run
```
Then verify:
1. Special characters display correctly as `[OK]` and `[X]`
2. No text overlapping in maneuver creation/edit dialogs
3. Dropdown boxes are positioned correctly in Hohmann transfer tab
4. Reset and Load buttons appear in info panel
## Remaining TODO Items from docs/TODO
- [ ] debug maneuver at periapsis - physics issue, not UI
- [ ] add reset/load new config UI control - **COMPLETED**
- [ ] UI fixes - **COMPLETED**
- [ ] SOI boundary testing - visualization feature
- [ ] interplanetary transfers - physics feature
## Next Steps
1. Test the UI fixes with the simulation
2. Consider adding actual config file loading functionality
3. Evaluate SOI boundary visualization if physics is working correctly
4. Continue with interplanetary transfers or maneuver debugging as appropriate

782
docs/technical_reference.md

@ -1,44 +1,37 @@
# Orbital Mechanics Simulation - Conceptual Reference # Technical Reference - Orbital Mechanics Simulation
## Overview ## 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.
2-body orbital mechanics simulator using **analytical propagation** for precise Keplerian trajectories. Supports elliptical, parabolic, and hyperbolic orbits with dynamic Sphere of Influence (SOI) transitions, impulsive burns, and 3D visualization via Raylib.
## Architecture ## 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.
Modular C-style C++ (structs/functions, no classes). Module dependencies:
```
Main → Simulation → Orbital Mechanics → Physics
→ Maneuver
→ Rendezvous
→ Config Loader
```
Renderer and UI Renderer depend on Config Loader for display data. Test Utilities are standalone.
## Coordinate Frame System
**Local Frame**: Position/velocity relative to parent body. Primary system for orbital mechanics calculations. Maintains double-precision accuracy for small orbits (LEO, lunar).
**Global Frame**: Position/velocity from simulation origin (Sun at 0,0,0). Computed as `global = parent.global + local`. Used for SOI distance calculations and rendering.
**Conversion**: Simple vector addition/subtraction. All bodies and spacecraft store both local and global coordinates.
## Core Data Structures ## Core Data Structures
### Vec3 ### Vec3 (physics.h)
3D vector (x, y, z). Operations: add, sub, cross, scale, magnitude, distance, normalize, dot. 3D vector for position, velocity, and acceleration
```cpp
struct Vec3 {
double x, y, z;
};
```
### Mat3 ### Mat3 (physics.h)
3x3 row-major matrix. Used for orbital plane rotations (z-x-z Euler angles). 3x3 row-major matrix for coordinate transformations
```cpp
struct Mat3 {
double m00, m01, m02;
double m10, m11, m12;
double m20, m21, m22;
};
```
### OrbitalElements ### OrbitalElements (orbital_mechanics.h)
Keplerian elements with union for parabolic/hyperbolic distinction: Keplerian orbital elements using union for parabolic/hyperbolic distinction
```cpp ```cpp
struct OrbitalElements { struct OrbitalElements {
union { union {
double semi_major_axis; // elliptical (e<1), hyperbolic (e>1) double semi_major_axis; // elliptical (e<1) and hyperbolic (e>1)
double semi_latus_rectum; // parabolic (e≈1) double semi_latus_rectum; // parabolic (e≈1)
}; };
double eccentricity; double eccentricity;
@ -49,243 +42,534 @@ struct OrbitalElements {
}; };
``` ```
### CelestialBody ### CelestialBody (simulation.h)
Planet/moon with mass, radius, parent_index, color, orbital elements, local/global coordinates, SOI radius. Planet/moon with local and global coordinate frames
```cpp
### Spacecraft struct CelestialBody {
Similar to CelestialBody but without radius or SOI. Supports standalone (parent_index = -1) and orbiting spacecraft. char name[64];
double mass; // kg
### Maneuver double radius; // meters
Impulsive burn: name, craft_index, direction, delta_v, trigger_type, trigger_value, scheduled_dt, executed flag, executed_time. int parent_index; // index of gravitational parent (-1 for root)
float color[3]; // RGB color for rendering
### HohmannTransfer
Hohmann transfer calculation result: dv1 (first burn delta-v), dv2 (second burn delta-v), transfer_time (half-period of transfer ellipse), true_anomaly_2 (true anomaly at second burn, typically π).
### SimulationState
Top-level container: arrays of bodies/spacecraft/maneuvers (with counts and capacities), time, dt, config_name.
## Orbital Propagation Algorithm
Analytical propagation solves Kepler's equations exactly (no integration error, perfect energy conservation).
**Common Pattern**:
1. Compute mean motion: n = √(μ/|a|³)
2. Convert true anomaly to appropriate anomaly form
3. Advance mean anomaly: M_new = M + n·dt
4. Solve Kepler's equation for new anomaly
5. Convert back to true anomaly
6. Update elements
**Elliptical (e < 1)**: Newton-Raphson on E - e·sin(E) = M. Tolerance: 1e-10, max 50 iterations.
**Parabolic (|e-1| < PARABOLIC_TOLERANCE)**: Barker's equation D + D³/3 = M with closed-form solution using cubic roots. Mean motion: n = √(μ/p³) where p = semi_latus_rectum.
**Hyperbolic (e > 1)**: Newton-Raphson on H - e·sinh(H) = M. Mean motion: n = √(μ/(-a)³) using negative semi-major axis. Initial guess: H = log(2M/e) when e·sinh(M) > M, else H = M.
## Orbital Element Reconstruction
Elements are dynamically maintained and reconstructed from Cartesian state when:
- **SOI transitions**: Parent changes, local coordinates update
- **Burns**: Velocity changes impulsively
- **Velocity drift**: |v_local - v_expected| > 1e-6 m/s before propagation
Reconstruction uses `cartesian_to_orbital_elements()` which handles edge cases (near-circular e < 1e-10, near-parabolic). OrbitalElements orbit; // Keplerian elements from config
## Coordinate Transformations 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
Orbital plane orientation via z-x-z Euler rotation: double soi_radius; // sphere of influence radius (meters)
``` };
R_total = R_z(Ω) · R_x(i) · R_z(ω)
``` ```
Sequence: argument_of_periapsis (ω) → inclination (i) → longitude_of_ascending_node (Ω). ### Spacecraft (spacecraft.h)
Spacecraft similar to CelestialBody but without SOI or radius
**Forward**: `orbital_elements_to_cartesian()` applies rotation. ```cpp
**Reverse**: `cartesian_to_orbital_elements()` computes angles from angular momentum and eccentricity vectors. struct Spacecraft {
char name[64];
## Sphere of Influence Mechanics double mass;
int parent_index;
**SOI Radius**: Hill sphere approximation r_soi = a × (m/M)^(2/5).
**SOI Transitions**:
1. `find_dominant_body()` called before each physics update
2. For non-root bodies: stay if distance < parent.soi_radius, else switch to root
3. For root bodies: find closest body where distance < body.soi_radius
4. On parent change:
- Compute global position/velocity from old parent
- Update parent_index
- Compute new local coordinates
- Reconstruct orbital elements
**SOI in Loop**: Checked for all bodies and spacecraft before physics propagation.
## Maneuver System
### Burn Directions
- BURN_PROGRADE: normalized velocity vector
- BURN_RETROGRADE: opposite velocity
- BURN_NORMAL: normalized (position × velocity)
- BURN_ANTINORMAL: opposite normal
- BURN_RADIAL_IN: -normalized position
- BURN_RADIAL_OUT: normalized position
- BURN_CUSTOM: user-specified vector
### Exact Position Execution
Both trigger types use sub-step propagation to execute burns at the exact trigger time, eliminating quantization error:
1. **TRIGGER_TIME**: `check_maneuver_trigger()` computes `dt_to_burn = trigger_value - sim->time` (sim->time is step start). If the trigger falls within `[sim->time, sim->time + sim->dt]`, `scheduled_dt` is set to this offset.
2. **TRIGGER_TRUE_ANOMALY**: `check_maneuver_trigger()` converts current and target true anomaly to mean anomaly, computes delta-M, divides by mean motion to get `dt_needed`.
3. If `0 < scheduled_dt <= sim->dt`, trigger fires and `scheduled_dt` is set.
4. In `update_spacecraft_physics()`, for each spacecraft: check all pending maneuvers for that craft.
5. If a maneuver fires: propagate by `burn_dt` (scheduled_dt), execute burn, propagate remaining (`sim->dt - burn_dt`).
6. No separate maneuver execution step — all inline in the spacecraft propagation loop.
**TRIGGER_TIME**: `scheduled_dt` is 0 only when the trigger time has already passed (clamped to fire immediately). Otherwise it is the exact sub-step offset from step start.
**TRIGGER_TRUE_ANOMALY**: Sub-step timing via analytical mean anomaly calculation.
**Future TODO**: Parabolic (Barker's equation) and hyperbolic branches for `check_maneuver_trigger()`.
### Hohmann Transfer
`calculate_hohmann_transfer()` computes optimal two-burn transfer between two circular orbits using the vis-viva equation. Transfer time equals half the period of the transfer ellipse.
## Rendezvous Module
Handles orbital rendezvous planning and execution via Hohmann transfers and phasing maneuvers.
**Hohmann Transfer Planning**:
- `calculate_required_separation_for_hohmann()` - computes ideal angular separation between chaser and target at first burn
- `calculate_wait_time_for_hohmann()` - determines wait time before starting transfer (positive = wait, negative = late)
- `calculate_next_hohmann_wait_time()` - like above but always returns non-negative by advancing to next phasing opportunity
- `calculate_relative_orbit_period()` - time between consecutive phasing opportunities
**Maneuver Creation**:
- `create_hohmann_departure_maneuver()` - adds prograde burn to enter transfer orbit (supports true anomaly or immediate trigger)
- `create_hohmann_arrival_maneuver()` - adds circularization burn to match target orbit (triggered at target radius)
**Verification**:
- `verify_hohmann_transfer_orbit()` - checks if current orbit matches expected Hohmann transfer parameters
- `validate_hohmann_transfer_parameters()` - validates transfer parameters before calculation
- `hohmann_transfer_complete()` - checks if transfer time has elapsed and spacecraft is at target radius
## Simulation Loop
**Initialization**:
1. create_simulation()
2. load_system_config() - parse TOML
3. run_all_config_validations()
4. initialize_orbital_objects() - convert elements to Cartesian, compute globals, calculate SOI
5. Main loop begins
**Main Loop Order**:
1. update_bodies_physics() - SOI checks, drift detection, propagation
2. compute_global_coordinates()
3. update_spacecraft_physics() - maneuver checking, propagation, burns
4. compute_spacecraft_globals()
5. time += dt
**Spacecraft Physics Per-Frame** (`update_spacecraft_physics`):
- For each spacecraft:
1. Validate local velocity against expected Keplerian velocity; if vel_diff > 1e-6, recalculate orbital elements
2. Check all pending maneuvers for this craft — if a trigger fires:
a. Propagate by `burn_dt` (scheduled sub-step offset)
b. Execute the maneuver (apply delta-v)
c. Propagate remaining (`sim->dt - burn_dt`)
3. If no maneuver: propagate full `sim->dt`
**Body Physics Per-Frame** (`update_bodies_physics`):
- Check SOI via find_dominant_body()
- Handle transitions (compute global, update parent, compute local, reconstruct elements)
- Check velocity drift (> 1e-6 m/s) and reconstruct if needed
- Propagate via propagate_orbital_elements()
- Update local position/velocity
## Orbit Types
**Elliptical (e < 1)**: ε < 0, bound periodic motion. a = -μ/(2ε) > 0.
**Parabolic (|e-1| < tolerance)**: ε 0, escape trajectory. p = h²/μ.
**Hyperbolic (e > 1)**: ε > 0, unbound. a = -μ/(2ε) < 0.
All satisfy vis-viva: v² = μ(2/r - 1/a).
## Configuration
**TOML Format**: [[bodies]], [[spacecraft]], [[maneuvers]] arrays.
**Validation Rules**:
- Parent index ordering (parents before children)
- Orbital element consistency (e ≥ 0, parabolic p > 0, etc.)
- True anomaly ranges for hyperbolic orbits
- Mass ratios (MIN_MASS_RATIO = 1000.0 for large-radius bodies)
- SOI overlap (siblings must not overlap)
- Nested orbits (grandchild semi-major axis ≤ 5× parent SOI)
- Initial positions (no overlap with parent)
## Testing Utilities
**OrbitalMetrics**: kinetic_energy, potential_energy, total_energy, orbital_radius, velocity_magnitude, angular_position.
**OrbitTracker**: Tracks orbit completion via quadrant transitions and total rotation. 3D mode uses orbital elements for inclined orbits.
**Energy Functions**:
- KE = 0.5 × m × v²
- PE = -G × m1 × m2 / r (min distance clamped to 1.0)
- Total = ΣKE + ΣPE_pairs
## Visualization
**Renderer (Raylib)**:
- XY simulation plane → XZ render plane (90° rotation around X)
- Scale: 1e-9 (1 unit = 1 billion meters)
- Relative rendering when body selected
- Orbit paths with appropriate segment counts
- Wireframe spheres, child indicators
**UI Renderer (raygui)**:
- Bottom-left: Sim time, body count, FPS, controls, config
- Top-left: Objects list with selection
- Top-right: Selected object details
- Below info: Maneuver list
## Build System
**Targets**: make / make all (build + ctags), make raylib, make test, make test-build, make clean, make clean-all, make rebuild. OrbitalElements orbit;
**Dependencies**: g++ (C++14), raylib (git submodule), raygui (header-only), tomlc17, Catch2 (tests only). Vec3 global_position;
Vec3 global_velocity;
Vec3 local_position;
Vec3 local_velocity;
};
```
**Testing**: ### Maneuver (maneuver.h)
```bash Impulsive burn with execution trigger
make test # Build and run all tests ```cpp
make test-build # Rebuild test executable only 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;
};
``` ```
### Running Specific Tests ### 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];
};
```
The test binary is at `./build/orbit_test`. Catch2 v3.7.1 supports multiple filtering strategies: ## Coordinate Frame System
```bash ### Local Frame
./build/orbit_test '[analytical]' # By tag - Position/velocity relative to parent body
./build/orbit_test "*barker*" # By name (wildcard) - Used for orbital mechanics calculations
./build/orbit_test '*burn*' [hybrid] # Combined: name AND tag - Primary coordinate system for propagation
./build/orbit_test '[hohmann]~[energy]' # Exclude tag (~)
./build/orbit_test -s '[tag]' # With debug/INFO output
```
Tags are the most common filter. Multiple tags on the command line use AND logic (test must have all tags). The `~` operator excludes a tag. ### 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
### SCENARIO Sections ### 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
SCENARIO tests group related assertions under SECTION sub-tests. Each SECTION is an independent test case — setup code before all SECTIONs runs once per SECTION, so fixtures are recreated. This allows drilling into specific parts: ## 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
- 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: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:**
```bash ```bash
./build/orbit_test '[hohmann]' --section "First burn at perigee raises apogee" # Run all tests
``` make test
In `--list-tests` output, SCENARIO tests display with "Scenario: " prefix. When filtering by name, wildcards handle the prefix transparently. # Run specific test config
./orbit_test '[config_name]'
Use `WithinAbs(expected, tolerance)` for floating-point comparisons (NOT `Approx()`). # Test with extra debug output (shows INFO messages)
./orbit_test -s '[config_name]'
## Hybrid Documentation Strategy ```
This document provides high-level concepts and architecture. Per-module summary files (`src/*.summary.md`) contain detailed function documentation including signatures, parameters, formulas, and edge cases. Use both: this for the big picture, summaries for implementation details. **Test framework:**
- Catch2 for testing
- Use `WithinAbs(expected, tolerance)` for floating-point comparisons (NOT `Approx()`)
- Required header: `<catch2/matchers/catch_matchers_floating_point.hpp>`

204
scripts/compute_rendezvous_params.py

@ -1,204 +0,0 @@
#!/usr/bin/env python3
"""
Pre-compute Hohmann transfer rendezvous parameters for test validation.
Replicates the exact rendezvous module phasing logic from src/rendezvous.cpp.
Hardcoded from tests/test_rendezvous.toml no TOML parser needed.
Usage: python3 tests/compute_rendezvous_params.py
"""
import math
import sys
G = 6.67430e-11
# Central body
EARTH_MASS = 5.972e24
# Spacecraft orbits (from test_rendezvous.toml)
TARGET_R = 6.771e6 # 400 km altitude
TARGET_NU = 0.0
CHASER_R = 6.671e6 # 300 km altitude
CHASER_NU = 4.71238898038469 # 270 degrees
MU = G * EARTH_MASS
def calc_mean_motion(radius, mass):
"""n = sqrt(mu / a^3)"""
return math.sqrt(MU / (radius ** 3))
def hohmann_transfer_time(r1, r2, mass):
"""Half orbit of transfer ellipse."""
a_transfer = (r1 + r2) / 2.0
T_transfer = 2.0 * math.pi * math.sqrt(a_transfer ** 3 / MU)
return T_transfer / 2.0
def required_separation(r1, r2, mass):
"""
Required angular separation at first burn.
chaser_pos - target_pos = target_angle - pi
"""
transfer_time = hohmann_transfer_time(r1, r2, mass)
n2 = calc_mean_motion(r2, mass)
target_angle = n2 * transfer_time
return target_angle - math.pi
def normalize_angle_2pi(angle):
"""Normalize to [0, 2*pi)."""
while angle < 0.0:
angle += 2.0 * math.pi
while angle >= 2.0 * math.pi:
angle -= 2.0 * math.pi
return angle
def normalize_angle_pi(angle):
"""Normalize to [-pi, pi]."""
angle = normalize_angle_2pi(angle)
while angle > math.pi:
angle -= 2.0 * math.pi
while angle < -math.pi:
angle += 2.0 * math.pi
return angle
def calculate_wait_time_for_hohmann(r1, r2, angular_separation, mass):
"""
Wait time before Hohmann transfer.
Positive = wait, negative = transfer already late.
"""
required_sep = required_separation(r1, r2, mass)
n1 = calc_mean_motion(r1, mass)
n2 = calc_mean_motion(r2, mass)
rel_angular_vel = n1 - n2
current_sep = normalize_angle_pi(angular_separation)
required_sep = normalize_angle_pi(required_sep)
angle_to_close = required_sep - current_sep
return angle_to_close / rel_angular_vel
def calculate_relative_orbit_period(r1, r2, mass):
"""Time between consecutive phasing opportunities."""
n1 = calc_mean_motion(r1, mass)
n2 = calc_mean_motion(r2, mass)
rel_angular_vel = abs(n1 - n2)
return 2.0 * math.pi / rel_angular_vel
def calculate_next_hohmann_wait_time(r1, r2, angular_separation, mass, min_wait_time):
"""
Like calculate_wait_time_for_hohmann, but advances to next phasing
opportunity if wait_time < min_wait_time. Always returns non-negative.
"""
wait_time = calculate_wait_time_for_hohmann(r1, r2, angular_separation, mass)
rel_period = calculate_relative_orbit_period(r1, r2, mass)
while wait_time < min_wait_time:
wait_time += rel_period
return wait_time
def main():
print(f"Central body: Earth, mass = {EARTH_MASS:.6e} kg")
print(f"mu = {MU:.6e} m^3/s^2")
print(f"\n=== INITIAL ORBITAL ELEMENTS ===")
print(f"Chaser_Lower: r = {CHASER_R:.6e} m, nu = {CHASER_NU:.6f} rad ({math.degrees(CHASER_NU):.2f} deg)")
print(f"Target: r = {TARGET_R:.6e} m, nu = {TARGET_NU:.6f} rad ({math.degrees(TARGET_NU):.2f} deg)")
# Angular separation: chaser - target
angular_sep = CHASER_NU - TARGET_NU
angular_sep = normalize_angle_pi(angular_sep)
print(f"\nAngular separation (chaser - target): {angular_sep:.6f} rad ({math.degrees(angular_sep):.2f} deg)")
# Mean motions
n1 = calc_mean_motion(CHASER_R, EARTH_MASS)
n2 = calc_mean_motion(TARGET_R, EARTH_MASS)
print(f"\nMean motions:")
print(f" n1 (chaser): {n1:.10f} rad/s")
print(f" n2 (target): {n2:.10f} rad/s")
print(f" n1 - n2: {n1 - n2:.10f} rad/s")
# Orbital periods
p_chaser = 2.0 * math.pi / n1
p_target = 2.0 * math.pi / n2
print(f"\nOrbital periods:")
print(f" Chaser: {p_chaser:.2f} s ({p_chaser/3600:.2f} h)")
print(f" Target: {p_target:.2f} s ({p_target/3600:.2f} h)")
# Hohmann transfer
tt = hohmann_transfer_time(CHASER_R, TARGET_R, EARTH_MASS)
a_t = (CHASER_R + TARGET_R) / 2.0
print(f"\n=== HOHMANN TRANSFER ===")
print(f" Transfer semi-major axis: {a_t:.6e} m")
print(f" Transfer time: {tt:.6f} s ({tt/60:.2f} min)")
# Required separation
req_sep = required_separation(CHASER_R, TARGET_R, EARTH_MASS)
req_sep_norm = normalize_angle_pi(req_sep)
print(f"\n=== REQUIRED SEPARATION ===")
print(f" Raw: {req_sep:.6f} rad ({math.degrees(req_sep):.2f} deg)")
print(f" Norm: {req_sep_norm:.6f} rad ({math.degrees(req_sep_norm):.2f} deg)")
# Relative orbit period
rel_period = calculate_relative_orbit_period(CHASER_R, TARGET_R, EARTH_MASS)
print(f"\nRelative orbit period: {rel_period:.6f} s ({rel_period/3600:.2f} h)")
# Detailed phasing calculation
print(f"\n=== PHASING CALCULATION ===")
current_sep = normalize_angle_pi(angular_sep)
print(f" Current separation (normalized): {current_sep:.6f} rad ({math.degrees(current_sep):.2f} deg)")
print(f" Required separation (normalized): {req_sep_norm:.6f} rad ({math.degrees(req_sep_norm):.2f} deg)")
angle_to_close = req_sep_norm - current_sep
print(f" Angle to close: {angle_to_close:.6f} rad ({math.degrees(angle_to_close):.2f} deg)")
wait_time = calculate_wait_time_for_hohmann(CHASER_R, TARGET_R, angular_sep, EARTH_MASS)
print(f" Raw wait_time: {wait_time:.6f} s ({wait_time/3600:.2f} h)")
# Wait times for various DT values
dt_values = [0.1, 0.5, 1.0, 2.0, 5.0, 10.0]
print(f"\n=== WAIT TIME vs DT (via calculate_next_hohmann_wait_time) ===")
for dt in dt_values:
wt = calculate_next_hohmann_wait_time(CHASER_R, TARGET_R, angular_sep, EARTH_MASS, dt)
arrival = wt + tt
steps = int(arrival / dt) + 1
print(f" DT={dt:6.1f} s: wait={wt:12.2f} s arrival={arrival:12.2f} s steps~{steps}")
# Recommended values for TIME_STEP = 0.1
dt = 0.1
wt = calculate_next_hohmann_wait_time(CHASER_R, TARGET_R, angular_sep, EARTH_MASS, dt)
arrival = wt + tt
max_steps = int(arrival / dt) + 1000
print(f"\n=== RECOMMENDED FOR TEST (DT=0.1) ===")
print(f" wait_time: {wt:.2f} s")
print(f" arrival_time: {arrival:.2f} s")
print(f" expected_steps: {int(arrival / dt)}")
print(f" max_steps (with margin): {max_steps}")
print(f" safety_limit (1 yr): {3600.0 * 24.0 * 365.0:.2f} s")
print(f"\n Milestone step indices:")
print(f" just_before_departure: {int(wt / dt)}")
print(f" after_departure: {int(wt / dt) + 1}")
print(f" just_before_arrival: {int(arrival / dt)}")
# Verify against C++ test output
print(f"\n=== COMPARISON WITH C++ TEST OUTPUT ===")
print(f" Python wait_time: {wt:.2f} s")
print(f" C++ test wait_time: 60062.7 s")
print(f" Python arrival: {arrival:.2f} s")
print(f" C++ test arrival: 62804.5 s")
print(f" Match: {abs(wt - 60062.7) < 0.1 and abs(arrival - 62804.5) < 0.1}")
if __name__ == '__main__':
main()

72
scripts/generate_interface_summaries.sh

@ -1,72 +0,0 @@
#!/bin/bash
# Generate interface summaries for all modules in src/
# Creates temporary markdown files alongside source files for quick reference
set -e
SRC_DIR="/home/agent/dev/claudes_game/src"
TEMP_EXT=".summary.md"
PROMPT_FILE="scripts/prompt.md"
# Function to generate prompt for a specific module
generate_prompt() {
local module_name="$1"
sed "s|\${module_name}|$module_name|g" "$PROMPT_FILE" | sed "s|\${SRC_DIR}|$SRC_DIR|g"
}
# Function to process a single module
process_module() {
local module_name="$1"
local header_file="${SRC_DIR}/${module_name}.h"
local output_file="${SRC_DIR}/${module_name}${TEMP_EXT}"
if [ ! -f "$header_file" ]; then
echo "Warning: No header file found for ${module_name}, skipping"
return
fi
echo "----------"
echo "Processing ${module_name}..."
# Build command as array
CMD=()
CMD+=("pi" "--print" "--no-session" "--provider" "llama.cpp" "--model" "Qwen3.6-35B")
CMD+=("$(generate_prompt "$module_name")")
#echo "CMD would execute: ${CMD[*]}"
"${CMD[@]}"
echo "Generated: ${output_file}"
}
# List of modules to process (based on src/ directory)
# Modules to exclude from processing
EXCLUDED=("renderer" "ui_renderer" "config_validator")
# Filter function
is_excluded() {
local mod="$1"
for exc in "${EXCLUDED[@]}"; do
[[ "$mod" == "$exc" ]] && return 0
done
return 1
}
# Build filtered list
MODULES=()
for h in src/*.h; do
mod=$(basename "$h" .h)
is_excluded "$mod" || MODULES+=("$mod")
done
echo "=== Interface Summary Generator ==="
echo "Processing ${#MODULES[@]} modules..."
echo ""
for module in "${MODULES[@]}"; do
process_module "$module"
done
echo ""
echo "=== Complete ==="
echo "Generated ${TEMP_EXT} files alongside source files in ${SRC_DIR}"
echo "These files are temporary and should be added to .gitignore"

38
scripts/prompt.md

@ -1,38 +0,0 @@
Analyze the ${SRC_DIR}/${module_name}.h and ${SRC_DIR}/${module_name}.cpp files in the orbital mechanics simulation project.
Create a concise but comprehensive summary of all interface elements (functions, structs, enums, constants, defines)
that would be useful for a new agent session to quickly understand this module's API.
Output format:
## Path
/full/file/path
## Structs, Constants, Enums, Defines
```cpp
// Brief description of the struct
struct <StructName> {
<type> <field_name>; // Field purpose
};
// comments (if needed)
static const <type> <CONSTANT_NAME> = <value>; // or inline (if needed)
```
## Functions
### <FunctionGroupName>
```cpp
// Brief description of what the function does, (include details for complex algorithms)
<return_type> <function_name>(
<param_type> <param_name>, // brief description
<param_type> <param_name> // brief description
);
```
Keep summaries concise but do not oversimplify complex algorithms - include enough detail
to understand the logic and approach. Focus on public interfaces, not internal implementation details.
IMPORTANT: write output to one file alongside the module location: ${SRC_DIR}/${module_name}.summary.md

431
scripts/simulate_rendezvous.py

@ -1,431 +0,0 @@
#!/usr/bin/env python3
"""
Full analytical propagation simulation of the Hohmann rendezvous scenario.
Replicates the exact physics from src/orbital_mechanics.cpp and src/maneuver.cpp.
Step-by-step trace to find where the 11,578 km separation comes from.
Usage: python3 tests/simulate_rendezvous.py
"""
import math
import sys
G = 6.67430e-11
MU = G * 5.972e24 # Earth
# ---- Vector operations ----
def vadd(a, b): return (a[0]+b[0], a[1]+b[1], a[2]+b[2])
def vsub(a, b): return (a[0]-b[0], a[1]-b[1], a[2]-b[2])
def vscale(v, s): return (v[0]*s, v[1]*s, v[2]*s)
def vmag(v): return math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
def vdot(a, b): return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
def vcross(a, b): return (
a[1]*b[2] - a[2]*b[1],
a[2]*b[0] - a[0]*b[2],
a[0]*b[1] - a[1]*b[0]
)
def vnorm(v):
m = vmag(v)
if m < 1e-15: return (0, 0, 0)
return (v[0]/m, v[1]/m, v[2]/m)
def normalize_angle(angle):
while angle < 0.0: angle += 2*math.pi
while angle >= 2*math.pi: angle -= 2*math.pi
return angle
def normalize_angle_2pi(angle):
while angle < 0.0: angle += 2*math.pi
while angle >= 2*math.pi: angle -= 2*math.pi
return angle
def normalize_angle_pi(angle):
angle = normalize_angle_2pi(angle)
while angle > math.pi: angle -= 2*math.pi
while angle < -math.pi: angle += 2*math.pi
return angle
# ---- Kepler equation solvers (exact C++ logic) ----
def get_initial_trial_value(mean_anomaly, eccentricity):
return (mean_anomaly + eccentricity * math.sin(mean_anomaly)
+ ((eccentricity**2 / 2.0) * math.sin(2.0 * mean_anomaly)))
def solve_kepler_elliptical(mean_anomaly, eccentricity):
E = get_initial_trial_value(mean_anomaly, eccentricity)
E_prev = E + 2.0e-10
for _ in range(50):
if abs(E - E_prev) < 1e-10:
break
E_prev = E
sin_E = math.sin(E)
E = E - (E - eccentricity * sin_E - mean_anomaly) / (1.0 - eccentricity * math.cos(E))
return E
def eccentric_to_true_anomaly(eccentric_anomaly, eccentricity):
if abs(1.0 - eccentricity) < 0.01:
E = eccentric_anomaly
e = eccentricity
cos_E = math.cos(E)
sin_E = math.sin(E)
denom = 1.0 - e * cos_E
cos_nu = max(-1.0, min(1.0, (cos_E - e) / denom))
sin_nu = max(-1.0, min(1.0, sin_E * math.sqrt(1.0 - e*e) / denom))
return math.atan2(sin_nu, cos_nu)
tan_half_E = math.tan(eccentric_anomaly / 2.0)
tan_half_nu = math.sqrt((1.0 + eccentricity) / (1.0 - eccentricity)) * tan_half_E
return 2.0 * math.atan(tan_half_nu)
# ---- Propagation (exact C++ propagate_orbital_elements) ----
def propagate(elements, dt, parent_mass):
a = elements['a']
e = elements['e']
nu = elements['nu']
mu = MU # fixed for this sim
if e < 1.0:
n = math.sqrt(mu / a**3)
E = 2.0 * math.atan(math.sqrt((1.0 - e) / (1.0 + e)) * math.tan(nu / 2.0))
M = E - e * math.sin(E)
M = M + n * dt
E_new = get_initial_trial_value(M, e)
E_prev = E_new + 2.0e-10
for _ in range(50):
if abs(E_new - E_prev) < 1e-10:
break
E_prev = E_new
sin_E = math.sin(E_new)
E_new = E_new - (E_new - e * sin_E - M) / (1.0 - e * math.cos(E_new))
nu_new = 2.0 * math.atan(math.sqrt((1.0 + e) / (1.0 - e)) * math.tan(E_new / 2.0))
result = dict(elements)
result['nu'] = nu_new
return result
else:
# Hyperbolic (not needed for this test)
raise NotImplementedError("hyperbolic propagation not needed")
# ---- Cartesian from orbital elements ----
def orbital_to_cartesian(elements, parent_mass):
a = elements['a']
e = elements['e']
nu = elements['nu']
inc = elements['inc']
Omega = elements['Omega']
omega = elements['omega']
mu = MU
p = a * (1.0 - e*e)
r = p / (1.0 + e * math.cos(nu))
# Orbital plane position/velocity
x_orb = r * math.cos(nu)
y_orb = r * math.sin(nu)
vx_orb = -math.sqrt(mu / p) * math.sin(nu)
vy_orb = math.sqrt(mu / p) * (e + math.cos(nu))
# z-x-z rotation: Rz(Omega) * Rx(inc) * Rz(omega)
# Apply Rz(omega) first
cos_w = math.cos(omega)
sin_w = math.sin(omega)
x1 = x_orb * cos_w - y_orb * sin_w
y1 = x_orb * sin_w + y_orb * cos_w
# Then Rx(inc)
cos_i = math.cos(inc)
sin_i = math.sin(inc)
x2 = x1
y2 = y1 * cos_i
z2 = y1 * sin_i
# Then Rz(Omega)
cos_O = math.cos(Omega)
sin_O = math.sin(Omega)
pos = (x2 * cos_O - y2 * sin_O,
x2 * sin_O + y2 * cos_O,
z2)
# Same rotation for velocity
vx1 = vx_orb * cos_w - vy_orb * sin_w
vy1 = vx_orb * sin_w + vy_orb * cos_w
vx2 = vx1
vy2 = vy1 * cos_i
vz2 = vy1 * sin_i
vel = (vx2 * cos_O - vy2 * sin_O,
vx2 * sin_O + vy2 * cos_O,
vz2)
return pos, vel
# ---- Cartesian to orbital elements ----
def cartesian_to_elements(pos, vel, parent_mass):
mu = MU
r = vmag(pos)
v = vmag(vel)
# Specific orbital energy
specific_energy = -mu / r + v**2 / 2.0
# Semi-major axis
if abs(specific_energy) < 1e-10:
a = 1e10
else:
a = -mu / (2.0 * specific_energy)
# Angular momentum
h_vec = vcross(pos, vel)
h = vmag(h_vec)
# Eccentricity vector
r_dot_v = vdot(pos, vel)
e_vec = ((v**2 - mu/r) * pos[0] - r_dot_v * vel[0]) / mu, \
((v**2 - mu/r) * pos[1] - r_dot_v * vel[1]) / mu, \
((v**2 - mu/r) * pos[2] - r_dot_v * vel[2]) / mu
e = vmag(e_vec)
# True anomaly
if e < 1e-10:
nu = 0.0
else:
cos_nu = vdot(pos, e_vec) / (r * e)
cos_nu = max(-1.0, min(1.0, cos_nu))
if abs(cos_nu) > 1.0 - 1e-10:
h_cross_e = vcross(h_vec, e_vec)
denom = r * e * h
sin_nu = vdot(pos, h_cross_e) / denom if denom > 1e-10 else 0.0
else:
r_cross_h = vcross(pos, h_vec)
denom = r * e * h
sin_nu = vdot(r_cross_h, e_vec) / denom if denom > 1e-10 else 0.0
nu = math.atan2(sin_nu, cos_nu)
if nu == -math.pi:
nu = math.pi
nu = normalize_angle(nu)
# Inclination
if h > 1e-10:
i = math.acos(h_vec[2] / h)
else:
i = 0.0
# RAAN
n_vec = (0, 0, 1)
n = vcross(n_vec, h_vec)
n_mag = vmag(n)
if n_mag > 1e-10:
Omega = math.acos(n[0] / n_mag)
if n[1] < 0.0:
Omega = 2*math.pi - Omega
else:
Omega = 0.0
# Argument of periapsis
if e > 1e-10 and n_mag > 1e-10 and i > 0.01:
cos_omega = vdot(e_vec, n) / (e * n_mag)
n_cross_e = vcross(n, e_vec)
sin_omega = vdot(n_cross_e, h_vec) / (e * n_mag * h)
omega = math.atan2(sin_omega, cos_omega)
if omega < 0: omega += 2*math.pi
elif e > 1e-10:
omega = math.atan2(e_vec[1], e_vec[0])
if omega < 0: omega += 2*math.pi
else:
omega = 0.0
return {'a': a, 'e': e, 'nu': nu, 'inc': i, 'Omega': Omega, 'omega': omega}
# ---- Hohmann transfer calculations ----
def hohmann_transfer_time(r1, r2):
a_t = (r1 + r2) / 2.0
T = 2*math.pi * math.sqrt(a_t**3 / MU)
return T / 2.0
def required_separation(r1, r2):
tt = hohmann_transfer_time(r1, r2)
n2 = math.sqrt(MU / r2**3)
target_angle = n2 * tt
return target_angle - math.pi
def calc_mean_motion(radius):
return math.sqrt(MU / radius**3)
def calculate_wait_time_for_hohmann(r1, r2, angular_separation):
required_sep = required_separation(r1, r2)
n1 = calc_mean_motion(r1)
n2 = calc_mean_motion(r2)
rel_angular_vel = n1 - n2
current_sep = normalize_angle_pi(angular_separation)
required_sep = normalize_angle_pi(required_sep)
angle_to_close = required_sep - current_sep
return angle_to_close / rel_angular_vel
def relative_orbit_period(r1, r2):
n1 = calc_mean_motion(r1)
n2 = calc_mean_motion(r2)
return 2*math.pi / abs(n1 - n2)
def calculate_next_hohmann_wait_time(r1, r2, angular_sep, dt):
wait_time = calculate_wait_time_for_hohmann(r1, r2, angular_sep)
rel_period = relative_orbit_period(r1, r2)
while wait_time < dt:
wait_time += rel_period
return wait_time
# ---- Burn application ----
def apply_burn(pos, vel, direction, delta_v, parent_mass):
"""Apply impulsive burn in local orbital frame."""
# direction: 'prograde', 'retrograde', 'normal'
if direction == 'prograde':
d = vnorm(vel)
elif direction == 'retrograde':
d = vscale(vnorm(vel), -1)
elif direction == 'normal':
h = vcross(pos, vel)
d = vnorm(h)
else:
raise ValueError(f"Unknown direction: {direction}")
new_vel = vadd(vel, vscale(d, delta_v))
return pos, new_vel
# ---- Dump state helper ----
def dump_state(label, chaser, target, chaser_pos, chaser_vel, target_pos, target_vel, sim_time):
"""Print state at key simulation milestones, matching test_rendezvous.cpp dump_state."""
c_r = vmag(chaser_pos)
t_r = vmag(target_pos)
c_sep = vmag(vsub(chaser_pos, target_pos))
print(f"\n*** {label} (t={sim_time:.1f}s) ***")
print(f" Chaser: r={c_r:.0f} m, nu={chaser['nu']:.6f} rad ({math.degrees(chaser['nu']):.1f}°) "
f"a={chaser['a']:.0f} e={chaser['e']:.6f}")
print(f" pos={chaser_pos}, vel={chaser_vel}")
print(f" Target: r={t_r:.0f} m, nu={target['nu']:.6f} rad ({math.degrees(target['nu']):.1f}°) "
f"a={target['a']:.0f} e={target['e']:.6f}")
print(f" pos={target_pos}, vel={target_vel}")
print(f" Separation: {c_sep:.0f} m")
# ---- Full rendezvous scenario ----
def main():
# Initial conditions from test_rendezvous.toml
TARGET_R = 6.771e6
TARGET_NU = 0.0
CHASER_R = 6.671e6
CHASER_NU = 4.71238898038469 # 270 degrees
print("=== INITIAL STATE ===")
print(f"Chaser: r={CHASER_R:.1f} m, nu={math.degrees(CHASER_NU):.1f} deg")
print(f"Target: r={TARGET_R:.1f} m, nu={math.degrees(TARGET_NU):.1f} deg")
# Create orbital elements (coplanar, circular)
chaser = {'a': CHASER_R, 'e': 0.0, 'nu': CHASER_NU,
'inc': 0.0, 'Omega': 0.0, 'omega': 0.0}
target = {'a': TARGET_R, 'e': 0.0, 'nu': TARGET_NU,
'inc': 0.0, 'Omega': 0.0, 'omega': 0.0}
chaser_pos, chaser_vel = orbital_to_cartesian(chaser, 5.972e24)
target_pos, target_vel = orbital_to_cartesian(target, 5.972e24)
print(f"Chaser pos: {chaser_pos}, vel: {vmag(chaser_vel):.1f} m/s")
print(f"Target pos: {target_pos}, vel: {vmag(target_vel):.1f} m/s")
# Angular separation
angular_sep = chaser['nu'] - target['nu']
angular_sep = normalize_angle_pi(angular_sep)
print(f"\nAngular separation (chaser - target): {math.degrees(angular_sep):.1f} deg")
# Hohmann parameters
hohmann_tt = hohmann_transfer_time(CHASER_R, TARGET_R)
dv1 = math.sqrt(MU * (2/CHASER_R - 2/(CHASER_R + TARGET_R))) - math.sqrt(MU/CHASER_R)
dv2 = math.sqrt(MU/TARGET_R) - math.sqrt(MU * (2/TARGET_R - 2/(CHASER_R + TARGET_R)))
print(f"\nHohmann transfer: tt={hohmann_tt:.1f} s, dv1={dv1:.2f} m/s, dv2={dv2:.2f} m/s")
# Phasing
dt = 0.1
wait_time = calculate_next_hohmann_wait_time(CHASER_R, TARGET_R, angular_sep, dt)
arrival_time = wait_time + hohmann_tt
print(f"Wait time: {wait_time:.2f} s")
print(f"Arrival time: {arrival_time:.2f} s")
print(f"Steps: {int(arrival_time/dt)}")
# ---- Run simulation ----
print(f"\n=== SIMULATION (dt={dt}) ===")
sim_time = 0.0
steps = 0
chaser_executed = False
arrival_executed = False
while steps < int(arrival_time / dt) + 1000:
chaser, target, chaser_pos, chaser_vel, target_pos, target_vel = \
update_simulation(chaser, target, sim_time, dt, dv1, dv2, wait_time, arrival_time,
chaser_pos, chaser_vel, target_pos, target_vel,
chaser_executed, arrival_executed)
sim_time += dt
steps += 1
if steps == 1:
dump_state("T=0 (initial)", chaser, target, chaser_pos, chaser_vel, target_pos, target_vel, sim_time)
if steps == int(wait_time / dt):
dump_state("JUST BEFORE DEPARTURE", chaser, target, chaser_pos, chaser_vel, target_pos, target_vel, sim_time)
if steps == int(wait_time / dt) + 1:
dump_state("AFTER DEPARTURE BURN", chaser, target, chaser_pos, chaser_vel, target_pos, target_vel, sim_time)
if steps == int(arrival_time / dt) - 1:
dump_state("JUST BEFORE ARRIVAL BURN", chaser, target, chaser_pos, chaser_vel, target_pos, target_vel, sim_time)
if not chaser_executed and sim_time >= wait_time:
# Execute departure burn
print(f"\n *** DEPARTURE BURN at t={sim_time:.1f}s ***")
print(f" Before: pos={chaser_pos}, vel={chaser_vel}")
chaser_pos, chaser_vel = apply_burn(chaser_pos, chaser_vel, 'prograde', dv1, 5.972e24)
print(f" After: pos={chaser_pos}, vel={chaser_vel}")
chaser = cartesian_to_elements(chaser_pos, chaser_vel, 5.972e24)
chaser_executed = True
print(f" Chaser: r={vmag(chaser_pos):.0f} nu={math.degrees(chaser['nu']):.1f}° "
f"a={chaser['a']:.0f} e={chaser['e']:.6f}")
if not arrival_executed and sim_time >= arrival_time:
# Execute arrival burn
chaser_pos, chaser_vel = apply_burn(chaser_pos, chaser_vel, 'prograde', dv2, 5.972e24)
chaser = cartesian_to_elements(chaser_pos, chaser_vel, 5.972e24)
arrival_executed = True
print(f"\n *** ARRIVAL BURN at t={sim_time:.1f}s ***")
print(f" Chaser: r={vmag(chaser_pos):.0f} nu={math.degrees(chaser['nu']):.1f}° "
f"a={chaser['a']:.0f} e={chaser['e']:.6f}")
dump_state("AFTER ARRIVAL BURN", chaser, target, chaser_pos, chaser_vel, target_pos, target_vel, sim_time)
# Final comparison
c_sep = vmag(vsub(chaser_pos, target_pos))
c_r = vmag(chaser_pos)
t_r = vmag(target_pos)
c_vel = vmag(chaser_vel)
t_vel = vmag(target_vel)
print(f"\n=== FINAL STATE ===")
print(f"Chaser: r={c_r:.0f} m, nu={chaser['nu']:.6f} rad ({math.degrees(chaser['nu']):.1f}°)")
print(f" pos={chaser_pos}, vel={chaser_vel}")
print(f"Target: r={t_r:.0f} m, nu={target['nu']:.6f} rad ({math.degrees(target['nu']):.1f}°)")
print(f" pos={target_pos}, vel={target_vel}")
print(f"Separation: {c_sep:.0f} m")
print(f"Speed: chaser={c_vel:.2f} target={t_vel:.2f} m/s")
print(f"Radius error: {abs(c_r - t_r):.6f} m")
print(f"Chaser eccentricity: {chaser['e']:.15f}")
print(f"Target eccentricity: {target['e']:.15f}")
break
def update_simulation(chaser, target, sim_time, dt, dv1, dv2, wait_time, arrival_time,
chaser_pos, chaser_vel, target_pos, target_vel,
chaser_executed, arrival_executed):
"""Propagate one timestep for both spacecraft."""
chaser = propagate(chaser, dt, 5.972e24)
target = propagate(target, dt, 5.972e24)
chaser_pos, chaser_vel = orbital_to_cartesian(chaser, 5.972e24)
target_pos, target_vel = orbital_to_cartesian(target, 5.972e24)
return chaser, target, chaser_pos, chaser_vel, target_pos, target_vel
if __name__ == '__main__':
main()

732
scripts/source_documentation_analysis.md

@ -1,732 +0,0 @@
# Source Documentation Analysis
## Executive Summary
Comprehensive examination of all source files in `src/` directory (15 files) to identify under-documented interface structs, enums, and functions. Focus areas: frequently used structures/enums and functions containing complex algorithms.
**Total under-documented items identified:** 24
**By priority:**
- **High (frequently used, complex algorithms):** 10 items
- **Medium (frequently used, important logic):** 10 items
- **Low (UI/rendering, less critical):** 4 items
**By type:**
- **Structs:** 5
- **Enums:** 2
- **Functions:** 17
---
## HIGH PRIORITY - Frequently Used Core Structures
### 1. OrbitalElements (orbital_mechanics.h)
**Status:** Under-documented despite being used throughout the codebase
**Missing documentation:**
- Union field usage explanation (why semi_major_axis vs semi_latus_rectum)
- Parameter meanings and units for all 6 elements
- Which elements are required vs optional in config files
- Relationship to classical Keplerian elements
- Parabolic orbit handling rationale
**Usage frequency:** Extremely high - used in simulation.cpp, maneuver.cpp, renderer.cpp, config_loader.cpp
**Current code:**
```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;
};
```
---
### 2. CelestialBody (simulation.h)
**Status:** Partially documented but missing key details
**Missing documentation:**
- When local vs global coordinates are used
- SOI radius calculation method and significance
- Parent index hierarchy rules
- Color format specification (RGB range)
- Velocity frame of reference
**Usage frequency:** Extremely high - core simulation structure
**Current code:**
```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)
};
```
---
### 3. SimulationState (simulation.h)
**Status:** Missing documentation
**Missing documentation:**
- Memory ownership (who allocates/frees)
- Thread safety (if applicable)
- Maximum time step recommendations
- Config name field usage
- Relationship between max_* and actual counts
**Usage frequency:** Extremely high - passed to almost every function
**Current code:**
```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];
};
```
---
## HIGH PRIORITY - Core Enums
### 4. BurnDirection (maneuver.h)
**Status:** Missing documentation
**Missing documentation:**
- Local frame definition (which axes)
- Physical meaning of each direction
- When to use each direction (mission planning context)
- BURN_CUSTOM usage and requirements
**Usage frequency:** High - used in maneuver.cpp, ui_renderer.cpp, renderer.cpp
**Current code:**
```cpp
enum BurnDirection {
BURN_PROGRADE,
BURN_RETROGRADE,
BURN_NORMAL,
BURN_ANTINORMAL,
BURN_RADIAL_IN,
BURN_RADIAL_OUT,
BURN_CUSTOM
};
```
---
### 5. TriggerType (maneuver.h)
**Status:** Missing documentation
**Missing documentation:**
- Time trigger precision requirements
- True anomaly trigger tolerance (0.01 rad)
- Which trigger is more accurate for different scenarios
- Recommended use cases
**Usage frequency:** Medium - used in maneuver.cpp, config_loader.cpp
**Current code:**
```cpp
enum TriggerType {
TRIGGER_TIME,
TRIGGER_TRUE_ANOMALY
};
```
---
## HIGH PRIORITY - Complex Algorithm Functions
### 6. propagate_orbital_elements() (orbital_mechanics.h/cpp)
**Status:** Documented in technical_reference.md but missing in source
**Missing in source code:**
- Parameter units explanation
- Edge case handling (near-parabolic, near-circular)
- Numerical stability guarantees
- Expected accuracy per time step
- When to use vs RK4_step
**Current comment:** None in header, minimal in implementation
**Function signature:**
```cpp
OrbitalElements propagate_orbital_elements(OrbitalElements elements, double dt, double parent_mass);
```
**Algorithm complexity:** High - handles elliptical, parabolic, and hyperbolic orbits with different propagation methods
---
### 7. cartesian_to_orbital_elements() (orbital_mechanics.h/cpp)
**Status:** Documented in technical_reference.md but missing in source
**Missing in source code:**
- Near-parabolic detection threshold
- Near-circular handling rationale
- Numerical stability measures
- Output guarantees for edge cases
- When reconstruction is needed (after burns, SOI transitions)
**Current comment:** None in header, minimal in implementation
**Function signature:**
```cpp
OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, double parent_mass);
```
**Algorithm complexity:** High - requires vector cross products, eccentricity vector calculation, and special case handling
---
### 8. solve_kepler_elliptical() and solve_kepler_hyperbolic() (orbital_mechanics.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Initial guess strategy explanation
- Convergence criteria and tolerance
- Maximum iterations justification
- Expected iteration count for typical cases
- Failure mode behavior
**Current comment:** Initial guess formula in get_initial_trial_value() but no explanation
**Function signatures:**
```cpp
double solve_kepler_elliptical(double mean_anomaly, double eccentricity);
double solve_kepler_hyperbolic(double mean_anomaly, double eccentricity);
```
**Algorithm complexity:** High - Newton-Raphson iterative solver with convergence checking
---
### 9. find_dominant_body() (simulation.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- SOI transition detection algorithm
- Performance characteristics (O(n) worst case)
- Root body handling logic
- When transitions are checked (every frame)
- Performance impact of many bodies
**Current comment:** None
**Function signature:**
```cpp
int find_dominant_body(SimulationState* sim, int body_index);
```
**Algorithm complexity:** Medium - iterates through all bodies to find closest SOI-dominant body
---
### 10. check_maneuver_trigger() (maneuver.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- True anomaly trigger prediction algorithm
- scheduled_dt calculation method
- Wraparound crossing detection logic
- Time precision requirements
- Performance considerations
**Current comment:** Some inline comments but no function-level documentation
**Function signature:**
```cpp
bool check_maneuver_trigger(Maneuver* maneuver, Spacecraft* craft, SimulationState* sim);
```
**Algorithm complexity:** Medium - requires orbital element to true anomaly conversion and crossing detection
---
## MEDIUM PRIORITY - Frequently Used Functions
### 11. initialize_orbital_objects() (simulation.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Initialization sequence and order
- Error handling behavior
- When to call (after config load, before simulation)
- What happens on failed initialization
**Current comment:** "Initialize orbital objects from orbital elements" - too brief
**Function signature:**
```cpp
void initialize_orbital_objects(SimulationState* sim);
```
---
### 12. update_simulation() (simulation.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Complete update sequence explanation
- Sub-function call order rationale
- Time step integration
- Maneuver execution timing within frame
- When to call (main loop only)
**Current comment:** None
**Function signature:**
```cpp
void update_simulation(SimulationState* sim);
```
---
### 13. execute_pending_maneuvers() (simulation.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Exact position burn execution rationale
- scheduled_dt usage
- Remaining time propagation
- Spacecraft handling flag purpose
- Performance implications
**Current comment:** None
**Function signature:**
```cpp
void execute_pending_maneuvers(SimulationState* sim);
```
---
### 14. orbital_elements_to_cartesian() (orbital_mechanics.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Rotation matrix derivation (z-x-z Euler)
- Coordinate frame transformations
- Output frame of reference
- Numerical stability measures
- Performance characteristics
**Current comment:** None
**Function signature:**
```cpp
void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass, Vec3* out_pos, Vec3* out_vel);
```
**Algorithm complexity:** Medium - requires Euler angle rotations and anomaly conversions
---
## MEDIUM PRIORITY - Configuration/Validation
### 15. load_system_config() (config_loader.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- TOML format requirements
- Error handling behavior
- Memory allocation strategy
- Validation sequence
- What fields are required vs optional
**Current comment:** None
**Function signature:**
```cpp
bool load_system_config(SimulationState* sim, const char* filepath);
```
---
### 16. run_all_config_validations() (config_validator.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Complete validation list
- Validation order and dependencies
- Error reporting format
- Which validations are critical vs warnings
- Performance impact
**Current comment:** None
**Function signature:**
```cpp
bool run_all_config_validations(SimulationState* sim);
```
---
## MEDIUM PRIORITY - Physics Module
### 17. mat3_rotation_orbital() (physics.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Euler angle sequence (z-x-z)
- Angle order and units
- Output matrix structure
- Use cases
- Numerical stability
**Current comment:** None
**Function signature:**
```cpp
Mat3 mat3_rotation_orbital(double omega, double i, double Omega);
```
---
### 18. rk4_step() (physics.h/cpp)
**Status:** Documented in technical_reference.md but missing in source
**Missing in source code:**
- When it's available but not used
- Comparison to analytical propagation
- Performance vs accuracy trade-offs
- Integration with orbital_mechanics
**Current comment:** "RK4 inferior to analytical propagation" comment exists but insufficient
**Function signature:**
```cpp
void rk4_step(Vec3* out_accel, Vec3 position, Vec3 velocity, double dt, double parent_mass);
```
---
## MEDIUM PRIORITY - Maneuver Module
### 19. preview_burn_result() (maneuver.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- UI preview usage
- State vector reconstruction method
- Accuracy guarantees
- Performance characteristics
- When to call (before execution)
**Current comment:** None
**Function signature:**
```cpp
void preview_burn_result(Spacecraft* craft, BurnDirection direction, double delta_v,
Vec3* out_new_local_pos, Vec3* out_new_local_vel, SimulationState* sim);
```
---
### 20. calculate_hohmann_transfer() (maneuver.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- Hohmann transfer assumptions
- Required conditions (coplanar, circular)
- Output field meanings
- Accuracy limitations
- When to use vs other transfers
**Current comment:** None
**Function signature:**
```cpp
void calculate_hohmann_transfer(SimulationState* sim, int current_body_index,
int target_body_index, HohmannTransfer* out_transfer);
```
---
## LOW PRIORITY - UI/Renderer
### 21. RenderState (renderer.h)
**Status:** Missing documentation
**Missing documentation:**
- Texture array usage (6 textures for burn directions)
- Camera follow mode state
- Offset preservation logic
- Last target index change detection
**Current comment:** None
**Current code:**
```cpp
struct RenderState {
Camera3D camera;
double distance_scale; // Scale factor for distances
double size_scale; // Scale factor for body sizes
int selected_body_index; // -1 = no selection
bool camera_target_enabled; // Whether camera follows selected body
Vector3 camera_offset; // Offset from target when following body
int last_target_index; // Tracks body index for change detection
Texture2D maneuver_textures[6]; // One per burn direction
bool texture_loaded;
};
```
---
### 22. UIState (ui_renderer.h)
**Status:** Missing documentation
**Missing documentation:**
- Buffer management strategy
- Cache invalidation rules
- Memory allocation/deallocation
- State persistence between frames
- Thread safety (if applicable)
**Current comment:** None
**Current code:**
```cpp
struct UIState {
int body_list_scroll;
int body_list_active;
int selected_craft_index;
ManeuverDialogState maneuver_dialog;
int maneuver_list_active;
int maneuver_list_scroll;
int selected_maneuver_index;
int cached_body_count;
int cached_craft_count;
int cached_maneuver_count;
int cached_executed_count;
int body_list_buffer_size;
char* body_list_buffer;
int craft_list_buffer_size;
char* craft_list_buffer;
int body_dropdown_buffer_size;
char* body_dropdown_buffer;
int maneuver_list_buffer_size;
char* maneuver_list_buffer;
};
```
---
### 23. ManeuverDialogState (ui_renderer.h)
**Status:** Missing documentation
**Missing documentation:**
- Tab state machine
- Preview calculation triggers
- Error message management
- Delete confirmation flow
- Field validation state
**Current comment:** None
**Current code:**
```cpp
struct ManeuverDialogState {
bool open;
ManeuverDialogTab active_tab;
int craft_index;
char name[64];
int direction_active;
double delta_v;
int trigger_type_active;
double trigger_value;
int edit_maneuver_index;
int target_body_index;
int current_body_index;
int transfer_body_index;
bool show_hohmann_preview;
HohmannTransfer last_calc;
bool show_preview;
bool preview_valid;
OrbitalElements preview_elements;
char error_message[256];
bool show_error;
bool show_delete_confirm;
};
```
---
## TEST UTILITIES
### 24. OrbitTracker (test_utilities.h/cpp)
**Status:** Missing documentation
**Missing documentation:**
- 3D angle calculation method
- Quadrant transition counting
- Orbit completion criteria
- min_time_days significance
- When to use vs simple angle tracking
**Current comment:** None
**Current code:**
```cpp
struct OrbitTracker {
double initial_angle;
double previous_angle;
int quadrant_transitions;
bool orbit_completed;
double time_at_completion;
int body_index;
double min_time_days;
// Orbital elements for 3D angle calculation
double inclination;
double longitude_of_ascending_node;
double argument_of_periapsis;
bool has_orbital_elements;
};
```
---
## Most Critical Gaps
The following 5 items are most critical for code maintainability and developer onboarding:
1. **OrbitalElements** - Core data structure used throughout the codebase
2. **propagate_orbital_elements()** - Main propagation algorithm with complex orbital mechanics
3. **cartesian_to_orbital_elements()** - Element reconstruction after burns and SOI transitions
4. **BurnDirection** - Mission planning interface with physical significance
5. **find_dominant_body()** - SOI transition detection algorithm
---
## Recommendations
### Immediate Actions (High Priority)
1. Add inline documentation to all HIGH PRIORITY items
2. Document parameter units and expected ranges
3. Explain algorithm choices and edge case handling
4. Add usage examples for complex functions
### Short-term Actions (Medium Priority)
1. Add function-level comments for MEDIUM PRIORITY items
2. Document initialization sequences and dependencies
3. Add performance notes where relevant
4. Document error handling behavior
### Long-term Actions (Low Priority)
1. Add comprehensive documentation for UI/Renderer structures
2. Document memory management patterns
3. Add thread safety notes (if applicable)
4. Create usage examples for complex UI state machines
---
## Impact Assessment
**Current state:** Code relies on technical_reference.md for documentation, but this is not visible to developers reading source files directly.
**Expected improvement:** Adding inline documentation will:
- Reduce onboarding time for new developers
- Improve code maintainability
- Reduce bugs from misunderstood algorithms
- Enable better code reviews
- Facilitate future feature additions
**Documentation coverage target:** 100% of public API functions and complex internal functions should have inline documentation explaining:
- Purpose and algorithm
- Parameter units and ranges
- Return value meaning
- Edge cases and error conditions
- Performance characteristics (where relevant)
---
## Methodology
This analysis was performed by:
1. Reading all 15 source files in `src/` directory
2. Identifying structs, enums, and functions
3. Assessing usage frequency across codebase
4. Evaluating algorithm complexity
5. Checking existing documentation (inline comments, technical_reference.md)
6. Prioritizing based on impact and frequency of use
**Files examined:**
- config_loader.h/cpp
- config_validator.h/cpp
- orbital_mechanics.h/cpp
- physics.h/cpp
- maneuver.h/cpp
- simulation.h/cpp
- spacecraft.h
- renderer.h/cpp
- ui_renderer.h/cpp
- main.cpp
- test_utilities.h/cpp
**Total lines of code analyzed:** ~8,500 lines

2
src/config_loader.h

@ -2,7 +2,7 @@
#define CONFIG_LOADER_H #define CONFIG_LOADER_H
#include "simulation.h" #include "simulation.h"
#include "orbital_objects.h" #include "spacecraft.h"
#include "../ext/tomlc17/src/tomlc17.h" #include "../ext/tomlc17/src/tomlc17.h"
#include "config_validator.h" #include "config_validator.h"

2
src/config_validator.h

@ -2,7 +2,7 @@
#define CONFIG_VALIDATOR_H #define CONFIG_VALIDATOR_H
#include "simulation.h" #include "simulation.h"
#include "orbital_objects.h" #include "spacecraft.h"
#define MIN_MASS_RATIO 1000.0 #define MIN_MASS_RATIO 1000.0
#define NESTED_ORBIT_FRACTION 5.0 #define NESTED_ORBIT_FRACTION 5.0

84
src/maneuver.cpp

@ -1,6 +1,6 @@
#include "maneuver.h" #include "maneuver.h"
#include "physics.h" #include "physics.h"
#include "orbital_objects.h" #include "spacecraft.h"
#include "simulation.h" #include "simulation.h"
#include "orbital_mechanics.h" #include "orbital_mechanics.h"
#include <cmath> #include <cmath>
@ -110,51 +110,23 @@ OrbitalElements preview_burn_result(const Spacecraft* craft, BurnDirection direc
return cartesian_to_orbital_elements(pos, new_vel, parent_mass); return cartesian_to_orbital_elements(pos, new_vel, parent_mass);
} }
// Elliptical orbits only: uses analytical mean anomaly delta to compute // Check if target angle lies between current and next angles (accounting for wraparound)
// exact time to target true anomaly, eliminating per-frame propagation. static bool angle_between(double current, double next, double target) {
// TODO: add parabolic (Barker's equation) and hyperbolic branches. double curr_norm = normalize_angle(current);
// TODO: Combine the two trigger code paths into one. Both ultimately need double next_norm = normalize_angle(next);
// to compute "dt from step start until the burn fires" and schedule a sub-step. double target_norm = normalize_angle(target);
// For TRIGGER_TRUE_ANOMALY, precalculate the absolute time at which the
// anomaly is reached (based on current orbital elements) and cache it on the
// Maneuver struct, so both paths reduce to: if (sim->time + sim->dt >= trigger_time)
// { scheduled_dt = min(trigger_time - sim->time, sim->dt); return true; }
//
// Complication: maneuvers for a single craft are stored in a flat array with no
// guaranteed order. When an earlier maneuver fires and changes the orbit, the
// cached trigger time for subsequent maneuvers becomes stale and must be
// recomputed. The current per-frame approach self-corrects by recomputing from
// the craft's current orbital elements each step.
bool check_maneuver_trigger(Maneuver* maneuver, Spacecraft* craft, SimulationState* sim) {
switch (maneuver->trigger_type) {
case TRIGGER_TIME: {
// Fire at the step that contains the trigger time.
// The orbit state is at sim->time (start of current step).
// We propagate forward to trigger_value, burn, then propagate
// the remaining time to reach sim->time + sim->dt.
if (sim->time > maneuver->trigger_value) {
// Trigger is before the start of this step — clamp to 0
// (should have fired in an earlier step; fire immediately)
maneuver->scheduled_dt = 0.0;
return true;
}
if (sim->time + sim->dt <= maneuver->trigger_value) {
return false;
}
double dt_to_burn = maneuver->trigger_value - sim->time;
// Clamp to valid range [0, sim->dt] if (curr_norm <= next_norm) {
if (dt_to_burn < 0.0) { return (target_norm >= curr_norm) && (target_norm <= next_norm);
dt_to_burn = 0.0; } else {
return (target_norm >= curr_norm) || (target_norm <= next_norm);
} }
if (dt_to_burn > sim->dt) {
dt_to_burn = sim->dt;
} }
maneuver->scheduled_dt = dt_to_burn; bool check_maneuver_trigger(Maneuver* maneuver, Spacecraft* craft, SimulationState* sim) {
return true; switch (maneuver->trigger_type) {
} case TRIGGER_TIME:
return sim->time >= maneuver->trigger_value;
case TRIGGER_TRUE_ANOMALY: { case TRIGGER_TRUE_ANOMALY: {
if (craft->parent_index < 0 || craft->parent_index >= sim->body_count) { if (craft->parent_index < 0 || craft->parent_index >= sim->body_count) {
@ -166,11 +138,37 @@ bool check_maneuver_trigger(Maneuver* maneuver, Spacecraft* craft, SimulationSta
double target_nu = normalize_angle(maneuver->trigger_value); double target_nu = normalize_angle(maneuver->trigger_value);
double current_diff = angular_distance(current_nu, target_nu); double current_diff = angular_distance(current_nu, target_nu);
#if 0
printf("INFO: Trigger check for '%s':\n", maneuver->name);
printf("INFO: Position: (%.2e, %.2e, %.2e) m\n",
r.x, r.y, r.z);
printf("INFO: current_nu (from orbit): %.6f rad\n", current_nu);
printf("INFO: target_nu: %.6f rad\n", target_nu);
printf("INFO: angular_distance: %.6f rad\n", current_diff);
#endif
if (current_diff < 0.01) { if (current_diff < 0.01) {
maneuver->scheduled_dt = 0.0; maneuver->scheduled_dt = 0.0;
return true; return true;
} }
OrbitalElements future_elements = propagate_orbital_elements(craft->orbit, sim->dt, parent->mass);
double future_nu = normalize_angle(future_elements.true_anomaly);
bool between = angle_between(current_nu, future_nu, target_nu);
if (!between) {
return false;
}
double future_diff = angular_distance(future_nu, target_nu);
bool wraparound_crossing = (current_nu > 5.0 && future_nu < 1.0) ||
(current_nu < 1.0 && future_nu > 5.0);
if (future_diff > current_diff && !wraparound_crossing) {
return false;
}
double a = craft->orbit.semi_major_axis; double a = craft->orbit.semi_major_axis;
double e = craft->orbit.eccentricity; double e = craft->orbit.eccentricity;
double mu = G * parent->mass; double mu = G * parent->mass;
@ -190,7 +188,7 @@ bool check_maneuver_trigger(Maneuver* maneuver, Spacecraft* craft, SimulationSta
dt_needed += M_period / n; dt_needed += M_period / n;
} }
if (dt_needed <= 0.0 || dt_needed > sim->dt) { if (dt_needed > sim->dt) {
return false; return false;
} }

2
src/maneuver.h

@ -2,7 +2,7 @@
#define MANEUVER_H #define MANEUVER_H
#include "physics.h" #include "physics.h"
#include "orbital_objects.h" #include "spacecraft.h"
struct SimulationState; struct SimulationState;

51
src/orbital_mechanics.cpp

@ -240,44 +240,10 @@ OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, doub
double r_mag = vec3_magnitude(r_vec); double r_mag = vec3_magnitude(r_vec);
double r_dot_e = vec3_dot(r_vec, e_vec); double r_dot_e = vec3_dot(r_vec, e_vec);
// Ascending node vector: n = k × h_vec (k is unit Z vector) // Near-circular: e ≈ 0, true anomaly undefined
Vec3 n_vec = {0.0, 0.0, 1.0};
Vec3 n = vec3_cross(n_vec, h_vec);
double n_mag = vec3_magnitude(n);
// True anomaly: angle from periapsis to position
double true_anomaly; double true_anomaly;
if (e < 1e-10) { if (e < 1e-10) {
// Circular orbit: no periapsis direction. Use argument of latitude true_anomaly = 0.0;
// (position angle in orbital plane) as true anomaly, with omega=0.
// For nearly-coplanar orbits, the ascending node is numerically
// unstable. Use the inclination to decide which reference to use.
double true_anomaly_from_position;
double sin_i = (h > 1e-10) ? n_mag / h : 1.0;
if (sin_i > 1e-6 && n_mag > 1e-10) {
// Well-defined ascending node: compute argument of latitude
double x_AN = n.x / n_mag;
double y_AN = n.y / n_mag;
// y_AN in orbital plane = (h × n) / |h × n|
double h_cross_n_x = h_vec.y * 0.0 - h_vec.z * n.y;
double h_cross_n_y = h_vec.z * n.x - h_vec.x * 0.0;
double h_cross_n_z = h_vec.x * n.y - h_vec.y * n.x;
double hcn_mag = sqrt(h_cross_n_x*h_cross_n_x + h_cross_n_y*h_cross_n_y + h_cross_n_z*h_cross_n_z);
if (hcn_mag > 1e-10) {
h_cross_n_x /= hcn_mag;
h_cross_n_y /= hcn_mag;
h_cross_n_z /= hcn_mag;
}
double r_xAN = r_vec.x * x_AN + r_vec.y * y_AN;
double r_yAN = r_vec.x * h_cross_n_x + r_vec.y * h_cross_n_y + r_vec.z * h_cross_n_z;
true_anomaly_from_position = atan2(r_yAN, r_xAN);
} else {
// Nearly coplanar: ascending node is numerically unstable.
// Use X-axis as reference. For coplanar orbits this gives
// the argument of latitude = atan2(y, x).
true_anomaly_from_position = atan2(r_vec.y, r_vec.x);
}
true_anomaly = normalize_angle(true_anomaly_from_position);
} else { } else {
double cos_nu = r_dot_e / (r_mag * e); double cos_nu = r_dot_e / (r_mag * e);
cos_nu = fmax(-1.0, fmin(1.0, cos_nu)); cos_nu = fmax(-1.0, fmin(1.0, cos_nu));
@ -312,6 +278,11 @@ OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, doub
i = 0.0; i = 0.0;
} }
// Ascending node vector: n = k × h_vec (k is unit Z vector)
Vec3 n_vec = {0.0, 0.0, 1.0};
Vec3 n = vec3_cross(n_vec, h_vec);
double n_mag = vec3_magnitude(n);
// Longitude of ascending node from n vector // Longitude of ascending node from n vector
double Omega; double Omega;
if (n_mag > 1e-10) { if (n_mag > 1e-10) {
@ -323,7 +294,6 @@ OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, doub
Omega = 0.0; Omega = 0.0;
} }
// Argument of periapsis: ω = atan2(n×e·h, e·n) // Argument of periapsis: ω = atan2(n×e·h, e·n)
// For coplanar orbits, use longitude of periapsis (angle of eccentricity vector)
double omega; double omega;
double inclination_threshold = 0.01; double inclination_threshold = 0.01;
@ -334,12 +304,6 @@ OrbitalElements cartesian_to_orbital_elements(Vec3 position, Vec3 velocity, doub
double sin_omega = vec3_dot(n_cross_e, h_vec) / (e * n_mag * h); double sin_omega = vec3_dot(n_cross_e, h_vec) / (e * n_mag * h);
omega = atan2(sin_omega, cos_omega); omega = atan2(sin_omega, cos_omega);
if (omega < 0.0) {
omega += 2.0 * M_PI;
}
} else if (e > 1e-10) {
// Coplanar or near-circular: use longitude of periapsis
omega = atan2(e_vec.y, e_vec.x);
if (omega < 0.0) { if (omega < 0.0) {
omega += 2.0 * M_PI; omega += 2.0 * M_PI;
} }
@ -404,6 +368,7 @@ OrbitalElements propagate_orbital_elements(const OrbitalElements& elements, doub
OrbitalElements result = elements; OrbitalElements result = elements;
result.true_anomaly = 2.0 * atan(sqrt((1.0 + e) / (1.0 - e)) * tan(E_new / 2.0)); result.true_anomaly = 2.0 * atan(sqrt((1.0 + e) / (1.0 - e)) * tan(E_new / 2.0));
return result; return result;
} else { // e >= 1.0 (hyperbolic) } else { // e >= 1.0 (hyperbolic)
double n = sqrt(mu / pow(-a, 3.0)); double n = sqrt(mu / pow(-a, 3.0));

49
src/orbital_objects.h

@ -1,49 +0,0 @@
#ifndef ORBITAL_OBJECTS_H
#define ORBITAL_OBJECTS_H
#include "physics.h"
#include "orbital_mechanics.h"
// Represents a planet, star, moon, or other gravitational body in the
// simulation. Supports hierarchical orbital mechanics with parent-child
// relationships and sphere of influence calculations.
struct CelestialBody {
char name[64];
double mass; // kg
double radius; // meters
int parent_index; // index of gravitational parent (-1 for root body like Sun)
float color[3]; // RGB color for rendering
// Orbital elements from config (Keplerian elements)
OrbitalElements orbit;
// Global frame (from origin)
Vec3 global_position; // meters from origin
Vec3 global_velocity; // m/s
// Local frame (relative to parent)
Vec3 local_position; // meters from parent
Vec3 local_velocity; // m/s relative to parent
double soi_radius; // sphere of influence radius (meters)
};
struct Spacecraft { // Spacecraft or probe in the simulation
char name[64];
double mass;
int parent_index;
// Orbital elements from config
OrbitalElements orbit;
// Global frame (from origin)
Vec3 global_position;
Vec3 global_velocity;
// Local frame (relative to parent)
Vec3 local_position;
Vec3 local_velocity;
};
#endif // ORBITAL_OBJECTS_H

2
src/renderer.h

@ -2,7 +2,7 @@
#define RENDERER_H #define RENDERER_H
#include "simulation.h" #include "simulation.h"
#include "orbital_objects.h" #include "spacecraft.h"
#include "maneuver.h" #include "maneuver.h"
#include "raylib.h" #include "raylib.h"

206
src/rendezvous.cpp

@ -1,206 +0,0 @@
#include "rendezvous.h"
#include <math.h>
#include <float.h>
// Mean motion: n = sqrt(mu / a^3)
static double calc_mean_motion(double radius, double mass) {
double mu = G * mass;
return sqrt(mu / pow(radius, 3));
}
// Hohmann transfer time (half orbit of transfer ellipse)
static double hohmann_transfer_time(double r1, double r2, double mass) {
double mu = G * mass;
double a_transfer = (r1 + r2) / 2.0;
double T_transfer = 2.0 * M_PI * sqrt(pow(a_transfer, 3) / mu);
return T_transfer / 2.0;
}
// Calculate required angular separation at first burn
// For Hohmann transfer: target should be at specific angle when chaser burns
// Returns: required angular separation (chaser - target) in radians
// Negative value means chaser should be behind target, positive means ahead
static double required_separation(double r1, double r2, double mass) {
double transfer_time = hohmann_transfer_time(r1, r2, mass);
double n2 = calc_mean_motion(r2, mass);
double target_angle = n2 * transfer_time;
// Chaser travels π radians in transfer orbit, target travels target_angle
// For rendezvous: chaser_pos + π = target_pos + target_angle
// Therefore: chaser_pos - target_pos = target_angle - π
// Negative value means chaser should be behind target
return target_angle - M_PI;
}
// Normalize angle to [0, 2π)
static double normalize_angle_2pi(double angle) {
while (angle < 0.0) {
angle += 2.0 * M_PI;
}
while (angle >= 2.0 * M_PI) {
angle -= 2.0 * M_PI;
}
return angle;
}
// Normalize angle to [-π, π] for shortest path
static double normalize_angle_pi(double angle) {
angle = normalize_angle_2pi(angle);
while (angle > M_PI) {
angle -= 2.0 * M_PI;
}
while (angle < -M_PI) {
angle += 2.0 * M_PI;
}
return angle;
}
// Calculate wait time before starting Hohmann transfer
// Determines how long to wait before executing the first burn so that
// both chaser and target arrive at the interception point simultaneously.
// Returns: wait time in seconds. Positive = wait, negative = transfer already late
double calculate_wait_time_for_hohmann(
double initial_orbit_radius,
double target_orbit_radius,
double angular_separation,
double central_mass
) {
double required_sep = required_separation(initial_orbit_radius, target_orbit_radius, central_mass);
double n1 = calc_mean_motion(initial_orbit_radius, central_mass);
double n2 = calc_mean_motion(target_orbit_radius, central_mass);
double rel_angular_vel = n1 - n2;
// Normalize current separation to [-pi, pi]
// Positive = chaser ahead of target, negative = chaser behind target
double current_sep = normalize_angle_pi(angular_separation);
// Normalize required separation to [-pi, pi]
required_sep = normalize_angle_pi(required_sep);
// Angle to close: difference between required and current separation
// If current_sep > required_sep, chaser is too far ahead (negative wait time)
// If current_sep < required_sep, chaser is too far behind (positive wait time)
double angle_to_close = required_sep - current_sep;
// Wait time = angle_to_close / relative_angular_velocity
return angle_to_close / rel_angular_vel;
}
// Calculate required angular separation for Hohmann transfer
// Computes the ideal angle between chaser and target at the moment
// of first burn to ensure simultaneous arrival at target orbit.
// Returns: required angular separation in radians (-2π, 2π)
double calculate_required_separation_for_hohmann(
double initial_orbit_radius,
double target_orbit_radius,
double central_mass
) {
double required_sep = required_separation(initial_orbit_radius, target_orbit_radius, central_mass);
return normalize_angle_pi(required_sep);
}
// Verify spacecraft is on correct Hohmann transfer orbit
// Checks if current orbit matches expected Hohmann transfer parameters.
// Returns: true if orbit is on Hohmann transfer, false otherwise
bool verify_hohmann_transfer_orbit(
const OrbitalElements* orbit,
double r1,
double r2,
double tolerance
) {
double expected_a = (r1 + r2) / 2.0;
double actual_a = orbit->semi_major_axis;
double diff = fabs(actual_a - expected_a);
return diff < tolerance;
}
// Check if Hohmann transfer is complete
// Determines if transfer time has elapsed and spacecraft is at target radius.
// Returns: true if transfer is complete, false otherwise
bool hohmann_transfer_complete(
double transfer_start_time,
double current_time,
double transfer_time,
double target_radius,
double current_radius,
double tolerance
) {
// Check if enough time has elapsed
if (current_time < transfer_start_time + transfer_time - tolerance) {
return false;
}
// Check if at target radius
double radius_diff = fabs(current_radius - target_radius);
return radius_diff < tolerance;
}
// Validate parameters for Hohmann transfer
// Checks if the transfer parameters are valid before calculation.
// Returns: true if parameters are valid, false otherwise
bool validate_hohmann_transfer_parameters(
double initial_orbit_radius,
double target_orbit_radius,
double central_mass
) {
// Check for positive radii
if (initial_orbit_radius <= 0.0 || target_orbit_radius <= 0.0) {
return false;
}
// Check for positive mass
if (central_mass <= 0.0) {
return false;
}
// Check for different orbits (no relative motion if equal)
if (initial_orbit_radius == target_orbit_radius) {
return false;
}
return true;
}
// Calculate relative orbit period
// Computes the time for the chaser to complete one full relative orbit
// with respect to the target (time between consecutive phasing opportunities).
// Returns: relative orbit period in seconds
double calculate_relative_orbit_period(
double initial_orbit_radius,
double target_orbit_radius,
double central_mass
) {
double n1 = calc_mean_motion(initial_orbit_radius, central_mass);
double n2 = calc_mean_motion(target_orbit_radius, central_mass);
double rel_angular_vel = fabs(n1 - n2);
return 2.0 * M_PI / rel_angular_vel;
}
// Calculate next valid wait time for Hohmann transfer
// Like calculate_wait_time_for_hohmann(), but always returns a non-negative
// value by advancing to the next phasing opportunity if needed.
// Returns: non-negative wait time in seconds (time to wait before executing transfer)
double calculate_next_hohmann_wait_time(
double initial_orbit_radius,
double target_orbit_radius,
double angular_separation,
double central_mass,
double min_wait_time
) {
double wait_time = calculate_wait_time_for_hohmann(
initial_orbit_radius, target_orbit_radius, angular_separation, central_mass
);
double rel_period = calculate_relative_orbit_period(
initial_orbit_radius, target_orbit_radius, central_mass
);
// Add relative orbit periods until wait_time >= min_wait_time
while (wait_time < min_wait_time) {
wait_time += rel_period;
}
return wait_time;
}

95
src/rendezvous.h

@ -1,95 +0,0 @@
#ifndef RENDEZVOUS_H
#define RENDEZVOUS_H
#include "orbital_mechanics.h"
#include "simulation.h"
#include "maneuver.h"
// Hohmann transfer calculations for coplanar circular orbits
// Provides functions to plan and execute Hohmann transfers between
// circular coplanar orbits, including phasing calculations to
// intercept moving targets.
// Calculate wait time before starting Hohmann transfer
// Determines how long to wait before executing the first burn so that
// both chaser and target arrive at the interception point simultaneously.
// Returns: wait time in seconds. Positive = wait, negative = transfer already late
double calculate_wait_time_for_hohmann(
double initial_orbit_radius, // Current circular orbit radius (meters)
double target_orbit_radius, // Target circular orbit radius (meters)
double angular_separation, // Current angle from chaser to target (radians, [0, 2π))
double central_mass // Mass of central body (kg)
);
// Calculate required angular separation for Hohmann transfer
// Computes the ideal angle between chaser and target at the moment
// of first burn to ensure simultaneous arrival at target orbit.
// Returns: required angular separation in radians (-2π, 2π)
double calculate_required_separation_for_hohmann(
double initial_orbit_radius, // Current circular orbit radius (meters)
double target_orbit_radius, // Target circular orbit radius (meters)
double central_mass // Mass of central body (kg)
);
// Create first burn maneuver for Hohmann transfer
// Adds a prograde burn maneuver to enter the Hohmann transfer orbit.
// Returns: maneuver index on success, -1 on failure
int create_hohmann_departure_maneuver(
SimulationState* sim, // Simulation state to add maneuver to
int chaser_index, // Index of chaser spacecraft
double delta_v, // Delta-v magnitude from calculate_hohmann_transfer (m/s)
double trigger_true_anomaly // True anomaly for trigger (-1 for immediate)
);
// Create second burn maneuver for Hohmann transfer
// Adds a circularization burn maneuver to match target orbit.
// Returns: maneuver index on success, -1 on failure
int create_hohmann_arrival_maneuver(
SimulationState* sim, // Simulation state to add maneuver to
int chaser_index, // Index of chaser spacecraft
double delta_v, // Delta-v magnitude from calculate_hohmann_transfer (m/s)
double trigger_true_anomaly // True anomaly at target orbit for trigger
);
// Verify spacecraft is on correct Hohmann transfer orbit
// Checks if current orbit matches expected Hohmann transfer parameters.
// Returns: true if orbit is on Hohmann transfer, false otherwise
bool verify_hohmann_transfer_orbit(
const OrbitalElements* orbit, // Current orbit after first burn
double r1, // Initial orbit radius (meters)
double r2, // Target orbit radius (meters)
double tolerance // Acceptable deviation in semi-major axis (meters)
);
// Validate parameters for Hohmann transfer
// Checks if the transfer parameters are valid before calculation.
// Returns: true if parameters are valid, false otherwise
bool validate_hohmann_transfer_parameters(
double initial_orbit_radius, // Current circular orbit radius (meters)
double target_orbit_radius, // Target circular orbit radius (meters)
double central_mass // Mass of central body (kg)
);
// Calculate relative orbit period
// Computes the time for the chaser to complete one full relative orbit
// with respect to the target (time between consecutive phasing opportunities).
// Returns: relative orbit period in seconds
double calculate_relative_orbit_period(
double initial_orbit_radius, // Current circular orbit radius (meters)
double target_orbit_radius, // Target circular orbit radius (meters)
double central_mass // Mass of central body (kg)
);
// Calculate next valid wait time for Hohmann transfer
// Like calculate_wait_time_for_hohmann(), but always returns a non-negative
// value by advancing to the next phasing opportunity if needed.
// Returns: non-negative wait time in seconds (time to wait before executing transfer)
double calculate_next_hohmann_wait_time(
double initial_orbit_radius, // Current circular orbit radius (meters)
double target_orbit_radius, // Target circular orbit radius (meters)
double angular_separation, // Current angle from chaser to target (radians)
double central_mass, // Mass of central body (kg)
double min_wait_time // Minimum acceptable wait time (seconds)
);
#endif // RENDEZVOUS_HOHMANN_H

69
src/simulation.cpp

@ -1,5 +1,5 @@
#include "simulation.h" #include "simulation.h"
#include "orbital_objects.h" #include "spacecraft.h"
#include "maneuver.h" #include "maneuver.h"
#include "orbital_mechanics.h" #include "orbital_mechanics.h"
#include <cassert> #include <cassert>
@ -8,6 +8,15 @@
#include <cstdio> #include <cstdio>
#include <cmath> #include <cmath>
// FIXME: I'm not a fan of this
static bool spacecraft_handled_this_frame[256];
static void reset_spacecraft_tracking(int max_craft) {
for (int i = 0; i < max_craft && i < 256; i++) {
spacecraft_handled_this_frame[i] = false;
}
}
// Create a new simulation // Create a new simulation
SimulationState* create_simulation(int max_bodies, int max_craft, int max_maneuvers, double time_step) { SimulationState* create_simulation(int max_bodies, int max_craft, int max_maneuvers, double time_step) {
SimulationState* sim = (SimulationState*)malloc(sizeof(SimulationState)); SimulationState* sim = (SimulationState*)malloc(sizeof(SimulationState));
@ -163,8 +172,10 @@ void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_ax
} }
void update_simulation(SimulationState* sim) { void update_simulation(SimulationState* sim) {
reset_spacecraft_tracking(sim->max_craft);
update_bodies_physics(sim); update_bodies_physics(sim);
compute_global_coordinates(sim); compute_global_coordinates(sim);
execute_pending_maneuvers(sim);
update_spacecraft_physics(sim); update_spacecraft_physics(sim);
compute_spacecraft_globals(sim); compute_spacecraft_globals(sim);
sim->time += sim->dt; sim->time += sim->dt;
@ -287,6 +298,10 @@ void update_spacecraft_physics(SimulationState* sim) {
continue; continue;
} }
if (spacecraft_handled_this_frame[i]) {
continue;
}
CelestialBody* parent = &sim->bodies[craft->parent_index]; CelestialBody* parent = &sim->bodies[craft->parent_index];
Vec3 expected_pos, expected_vel; Vec3 expected_pos, expected_vel;
@ -297,42 +312,46 @@ void update_spacecraft_physics(SimulationState* sim) {
craft->orbit = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, parent->mass); craft->orbit = cartesian_to_orbital_elements(craft->local_position, craft->local_velocity, parent->mass);
} }
// Check if any pending maneuver fires this frame. If so, propagate craft->orbit = propagate_orbital_elements(craft->orbit, sim->dt, parent->mass);
// to the burn time, execute the burn, then propagate the remaining orbital_elements_to_cartesian(craft->orbit, parent->mass, &craft->local_position, &craft->local_velocity);
// timestep. This enables sub-step interpolation for true-anomaly }
// triggers and sets up future interpolated time triggers. }
bool maneuver_fired = false;
double burn_dt = 0.0; void execute_pending_maneuvers(SimulationState* sim) {
Maneuver* fired_maneuver = NULL; for (int i = 0; i < sim->maneuver_count; i++) {
for (int j = 0; j < sim->maneuver_count; j++) { Maneuver* maneuver = &sim->maneuvers[i];
Maneuver* maneuver = &sim->maneuvers[j];
if (maneuver->executed) { if (maneuver->executed) {
continue; continue;
} }
if (maneuver->craft_index != i) {
if (maneuver->craft_index < 0 || maneuver->craft_index >= sim->craft_count) {
continue;
}
Spacecraft* craft = &sim->spacecraft[maneuver->craft_index];
if (craft->parent_index < 0 || craft->parent_index >= sim->body_count) {
continue; continue;
} }
if (check_maneuver_trigger(maneuver, craft, sim)) { if (check_maneuver_trigger(maneuver, craft, sim)) {
burn_dt = maneuver->scheduled_dt; CelestialBody* parent = &sim->bodies[craft->parent_index];
fired_maneuver = maneuver; double dt_to_burn = maneuver->scheduled_dt;
maneuver_fired = true;
break;
}
}
if (maneuver_fired) { if (dt_to_burn > 0.0) {
craft->orbit = propagate_orbital_elements(craft->orbit, burn_dt, parent->mass); craft->orbit = propagate_orbital_elements(craft->orbit, dt_to_burn, parent->mass);
orbital_elements_to_cartesian(craft->orbit, parent->mass, &craft->local_position, &craft->local_velocity); orbital_elements_to_cartesian(craft->orbit, parent->mass, &craft->local_position, &craft->local_velocity);
}
double burn_time = sim->time + burn_dt; execute_maneuver(maneuver, craft, sim, sim->time + dt_to_burn);
execute_maneuver(fired_maneuver, craft, sim, burn_time);
double remaining_dt = sim->dt - burn_dt; double remaining_dt = sim->dt - dt_to_burn;
craft->orbit = propagate_orbital_elements(craft->orbit, remaining_dt, parent->mass); craft->orbit = propagate_orbital_elements(craft->orbit, remaining_dt, parent->mass);
orbital_elements_to_cartesian(craft->orbit, parent->mass, &craft->local_position, &craft->local_velocity); orbital_elements_to_cartesian(craft->orbit, parent->mass, &craft->local_position, &craft->local_velocity);
} else {
craft->orbit = propagate_orbital_elements(craft->orbit, sim->dt, parent->mass); spacecraft_handled_this_frame[maneuver->craft_index] = true;
orbital_elements_to_cartesian(craft->orbit, parent->mass, &craft->local_position, &craft->local_velocity); maneuver->scheduled_dt = 0.0;
} }
} }
} }

30
src/simulation.h

@ -1,13 +1,34 @@
#ifndef BODIES_H #ifndef BODIES_H
#define BODIES_H #define BODIES_H
#include "orbital_objects.h" #include "physics.h"
#include "maneuver.h" #include "orbital_mechanics.h"
// Forward declarations struct Spacecraft;
struct SimulationState;
struct Maneuver; struct Maneuver;
// Celestial body structure
struct CelestialBody {
char name[64];
double mass; // kg
double radius; // meters
int parent_index; // index of gravitational parent (-1 for root body like Sun)
float color[3]; // RGB color for rendering
// Orbital elements from config (Keplerian elements)
OrbitalElements orbit;
// Global frame (from origin)
Vec3 global_position; // meters from origin
Vec3 global_velocity; // m/s
// Local frame (relative to parent)
Vec3 local_position; // meters from parent
Vec3 local_velocity; // m/s relative to parent
double soi_radius; // sphere of influence radius (meters)
};
// Simulation state // Simulation state
struct SimulationState { struct SimulationState {
CelestialBody* bodies; CelestialBody* bodies;
@ -47,6 +68,7 @@ void compute_global_coordinates(SimulationState* sim);
// Simulation update helper functions // Simulation update helper functions
void update_bodies_physics(SimulationState* sim); void update_bodies_physics(SimulationState* sim);
void update_spacecraft_physics(SimulationState* sim); void update_spacecraft_physics(SimulationState* sim);
void execute_pending_maneuvers(SimulationState* sim);
void compute_spacecraft_globals(SimulationState* sim); void compute_spacecraft_globals(SimulationState* sim);
// Initialize orbital objects from orbital elements // Initialize orbital objects from orbital elements

5
src/spacecraft.cpp

@ -0,0 +1,5 @@
#include "spacecraft.h"
#include "simulation.h"
#include <cstdio>
#include <cstring>
#include <cmath>

24
src/spacecraft.h

@ -0,0 +1,24 @@
#ifndef SPACECRAFT_H
#define SPACECRAFT_H
#include "physics.h"
#include "orbital_mechanics.h"
struct Spacecraft {
char name[64];
double mass;
int parent_index;
// Orbital elements from config
OrbitalElements orbit;
// Global frame (from origin)
Vec3 global_position;
Vec3 global_velocity;
// Local frame (relative to parent)
Vec3 local_position;
Vec3 local_velocity;
};
#endif

52
src/test_utilities.cpp

@ -1,7 +1,6 @@
#include "test_utilities.h" #include "test_utilities.h"
#include <cstdlib> #include <cstdlib>
#include <cmath> #include <cmath>
#include <cstdio>
double calculate_kinetic_energy(CelestialBody* body) { double calculate_kinetic_energy(CelestialBody* body) {
double v_squared = body->global_velocity.x * body->global_velocity.x + double v_squared = body->global_velocity.x * body->global_velocity.x +
@ -165,54 +164,3 @@ bool compare_vec3(Vec3 a, Vec3 b, double tolerance) {
fabs(a.y - b.y) <= tolerance && fabs(a.y - b.y) <= tolerance &&
fabs(a.z - b.z) <= tolerance; fabs(a.z - b.z) <= tolerance;
} }
int dump_simulation_state(SimulationState* sim, const char* label,
char* buffer, int buffer_size) {
int offset = 0;
offset += snprintf(buffer + offset, buffer_size - offset,
"\n=== %s (t=%.0f s) ===\n", label, sim->time);
offset += snprintf(buffer + offset, buffer_size - offset,
"Bodies (%d):\n", sim->body_count);
for (int i = 0; i < sim->body_count; i++) {
offset += snprintf(buffer + offset, buffer_size - offset,
" [%d] %s: mass=%.2e kg\n",
i, sim->bodies[i].name, sim->bodies[i].mass);
}
offset += snprintf(buffer + offset, buffer_size - offset,
"Spacecraft (%d):\n", sim->craft_count);
for (int i = 0; i < sim->craft_count; i++) {
Spacecraft* s = &sim->spacecraft[i];
double r = sqrt(s->local_position.x*s->local_position.x +
s->local_position.y*s->local_position.y +
s->local_position.z*s->local_position.z);
double v = sqrt(s->local_velocity.x*s->local_velocity.x +
s->local_velocity.y*s->local_velocity.y +
s->local_velocity.z*s->local_velocity.z);
offset += snprintf(buffer + offset, buffer_size - offset,
" [%d] %s: r=%.1f v=%.1f nu=%.5f a=%.1f e=%.6f, omega=%.6f\n",
i, s->name, r, v,
s->orbit.true_anomaly,
s->orbit.semi_major_axis,
s->orbit.eccentricity,
s->orbit.argument_of_periapsis);
offset += snprintf(buffer + offset, buffer_size - offset,
" pos=(%.1f, %.1f, %.1f) vel=(%.1f, %.1f, %.1f)\n",
s->local_position.x, s->local_position.y, s->local_position.z,
s->local_velocity.x, s->local_velocity.y, s->local_velocity.z);
}
offset += snprintf(buffer + offset, buffer_size - offset,
"Maneuvers (%d):\n", sim->maneuver_count);
for (int i = 0; i < sim->maneuver_count; i++) {
Maneuver* m = &sim->maneuvers[i];
offset += snprintf(buffer + offset, buffer_size - offset,
" [%d] %s: craft=%d dir=%d dv=%.4f trigger=%d val=%.2f exec=%d\n",
i, m->name, m->craft_index, m->direction, m->delta_v,
m->trigger_type, m->trigger_value, m->executed);
}
return offset;
}

6
src/test_utilities.h

@ -46,10 +46,4 @@ void destroy_orbit_tracker(OrbitTracker* tracker);
bool compare_double(double a, double b, double tolerance); bool compare_double(double a, double b, double tolerance);
bool compare_vec3(Vec3 a, Vec3 b, double tolerance); bool compare_vec3(Vec3 a, Vec3 b, double tolerance);
// Write simulation state to a caller-allocated buffer.
// Returns number of characters written (excluding null terminator).
// Caller must ensure buffer is large enough.
int dump_simulation_state(SimulationState* sim, const char* label,
char* buffer, int buffer_size);
#endif #endif

16
src/ui_renderer.cpp

@ -1,5 +1,5 @@
#include "ui_renderer.h" #include "ui_renderer.h"
#include "orbital_objects.h" #include "spacecraft.h"
#include "maneuver.h" #include "maneuver.h"
#include <cstdio> #include <cstdio>
#include <cstring> #include <cstring>
@ -668,7 +668,7 @@ void render_maneuver_create_tab(SimulationState* sim, UIState* ui_state, int x,
update_orbital_preview(sim, ui_state); update_orbital_preview(sim, ui_state);
if (ui_state->maneuver_dialog.show_preview) { if (ui_state->maneuver_dialog.show_preview) {
Rectangle preview_bounds = {(float)x + 8, (float)current_y + 50, (float)width - 16, 110}; Rectangle preview_bounds = {(float)x + 8, (float)current_y + 10, (float)width - 16, 110};
render_orbital_preview(sim, ui_state, preview_bounds); render_orbital_preview(sim, ui_state, preview_bounds);
} }
} }
@ -751,7 +751,7 @@ void render_orbital_preview(SimulationState* sim, UIState* ui_state, Rectangle b
y += 30; y += 30;
Rectangle status_bounds = {(float)x, (float)y, (float)width, 25}; Rectangle status_bounds = {(float)x, (float)y, (float)width, 25};
const char* status_text = ui_state->maneuver_dialog.preview_valid ? "[OK] All parameters valid" : "[X] Invalid parameters"; const char* status_text = ui_state->maneuver_dialog.preview_valid ? " All parameters valid" : " Invalid parameters";
GuiLabel(status_bounds, status_text); GuiLabel(status_bounds, status_text);
} }
@ -824,7 +824,7 @@ void render_maneuver_hohmann_tab(SimulationState* sim, UIState* ui_state, int x,
if (sim->craft_count > 0) { if (sim->craft_count > 0) {
static bool edit_mode_spacecraft = false; static bool edit_mode_spacecraft = false;
if (GuiDropdownBox({(float)x + 120, (float)current_y, (float)(width - 130), 25}, if (GuiDropdownBox({(float)x + 120, (float)(current_y - 35), (float)(width - 130), 25},
ui_state->craft_list_buffer, &ui_state->maneuver_dialog.craft_index, edit_mode_spacecraft)) { ui_state->craft_list_buffer, &ui_state->maneuver_dialog.craft_index, edit_mode_spacecraft)) {
edit_mode_spacecraft = !edit_mode_spacecraft; edit_mode_spacecraft = !edit_mode_spacecraft;
} }
@ -837,7 +837,7 @@ void render_maneuver_hohmann_tab(SimulationState* sim, UIState* ui_state, int x,
if (sim->body_count > 0) { if (sim->body_count > 0) {
static bool edit_mode_transfer_parent = false; static bool edit_mode_transfer_parent = false;
if (GuiDropdownBox({(float)x + 120, (float)current_y, (float)(width - 130), 25}, if (GuiDropdownBox({(float)x + 120, (float)(current_y - 35), (float)(width - 130), 25},
ui_state->body_dropdown_buffer, &ui_state->maneuver_dialog.transfer_body_index, edit_mode_transfer_parent)) { ui_state->body_dropdown_buffer, &ui_state->maneuver_dialog.transfer_body_index, edit_mode_transfer_parent)) {
edit_mode_transfer_parent = !edit_mode_transfer_parent; edit_mode_transfer_parent = !edit_mode_transfer_parent;
} }
@ -850,7 +850,7 @@ void render_maneuver_hohmann_tab(SimulationState* sim, UIState* ui_state, int x,
if (sim->body_count > 0) { if (sim->body_count > 0) {
static bool edit_mode_current_body = false; static bool edit_mode_current_body = false;
if (GuiDropdownBox({(float)x + 120, (float)current_y, (float)(width - 130), 25}, if (GuiDropdownBox({(float)x + 120, (float)(current_y - 35), (float)(width - 130), 25},
ui_state->body_dropdown_buffer, &ui_state->maneuver_dialog.current_body_index, edit_mode_current_body)) { ui_state->body_dropdown_buffer, &ui_state->maneuver_dialog.current_body_index, edit_mode_current_body)) {
edit_mode_current_body = !edit_mode_current_body; edit_mode_current_body = !edit_mode_current_body;
} }
@ -863,7 +863,7 @@ void render_maneuver_hohmann_tab(SimulationState* sim, UIState* ui_state, int x,
if (sim->body_count > 0) { if (sim->body_count > 0) {
static bool edit_mode_target_body = false; static bool edit_mode_target_body = false;
if (GuiDropdownBox({(float)x + 120, (float)current_y, (float)(width - 130), 25}, if (GuiDropdownBox({(float)x + 120, (float)(current_y - 35), (float)(width - 130), 25},
ui_state->body_dropdown_buffer, &ui_state->maneuver_dialog.target_body_index, edit_mode_target_body)) { ui_state->body_dropdown_buffer, &ui_state->maneuver_dialog.target_body_index, edit_mode_target_body)) {
edit_mode_target_body = !edit_mode_target_body; edit_mode_target_body = !edit_mode_target_body;
} }
@ -1431,7 +1431,7 @@ void render_maneuver_edit_tab(SimulationState* sim, UIState* ui_state, int x, in
update_orbital_preview(sim, ui_state); update_orbital_preview(sim, ui_state);
if (ui_state->maneuver_dialog.show_preview) { if (ui_state->maneuver_dialog.show_preview) {
Rectangle preview_bounds = {(float)x + 8, (float)current_y + 50, (float)width - 16, 110}; Rectangle preview_bounds = {(float)x + 8, (float)current_y + 10, (float)width - 16, 110};
render_orbital_preview(sim, ui_state, preview_bounds); render_orbital_preview(sim, ui_state, preview_bounds);
} }

2
tests/test_cartesian_to_elements_advanced.cpp

@ -2,7 +2,7 @@
#include <catch2/matchers/catch_matchers_floating_point.hpp> #include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <cmath> #include <cmath>
#include "../src/orbital_mechanics.h" #include "../src/orbital_mechanics.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/test_utilities.h" #include "../src/test_utilities.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include "../src/simulation.h" #include "../src/simulation.h"

2
tests/test_hybrid_burns.cpp

@ -3,7 +3,7 @@
#include "../src/physics.h" #include "../src/physics.h"
#include "../src/orbital_mechanics.h" #include "../src/orbital_mechanics.h"
#include "../src/simulation.h" #include "../src/simulation.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/maneuver.h" #include "../src/maneuver.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include "../src/test_utilities.h" #include "../src/test_utilities.h"

2
tests/test_hybrid_energy_conservation.cpp

@ -3,7 +3,7 @@
#include "../src/physics.h" #include "../src/physics.h"
#include "../src/orbital_mechanics.h" #include "../src/orbital_mechanics.h"
#include "../src/simulation.h" #include "../src/simulation.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/maneuver.h" #include "../src/maneuver.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include <cmath> #include <cmath>

9
tests/test_maneuver_planning.cpp

@ -1,16 +1,10 @@
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "../src/physics.h" #include "../src/physics.h"
#include "../src/simulation.h" #include "../src/simulation.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/maneuver.h" #include "../src/maneuver.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include "../src/rendezvous.h"
#include "../src/test_utilities.h"
#include <cmath> #include <cmath>
#include <cstring>
using Catch::Matchers::WithinAbs;
TEST_CASE("Maneuver loading from config", "[maneuver][config]") { TEST_CASE("Maneuver loading from config", "[maneuver][config]") {
const double TIME_STEP = 60.0; const double TIME_STEP = 60.0;
@ -149,4 +143,3 @@ TEST_CASE("Duplicate maneuver names fail config load", "[maneuver][config][error
destroy_simulation(sim); destroy_simulation(sim);
} }

2
tests/test_maneuvers.cpp

@ -1,7 +1,7 @@
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include "../src/physics.h" #include "../src/physics.h"
#include "../src/simulation.h" #include "../src/simulation.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/maneuver.h" #include "../src/maneuver.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include "../src/test_utilities.h" #include "../src/test_utilities.h"

13
tests/test_omega_debug.cpp

@ -2,7 +2,7 @@
#include <catch2/matchers/catch_matchers_floating_point.hpp> #include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "../src/physics.h" #include "../src/physics.h"
#include "../src/simulation.h" #include "../src/simulation.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/maneuver.h" #include "../src/maneuver.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include "../src/orbital_mechanics.h" #include "../src/orbital_mechanics.h"
@ -55,11 +55,10 @@ TEST_CASE("Omega calculation after prograde burn", "[omega][debug]") {
Vec3 e_vec_new = vec3_scale(vec3_sub(vec3_scale(vel, r), vec3_scale(pos, vec3_magnitude(vel) / mu)), 1.0 / mu); Vec3 e_vec_new = vec3_scale(vec3_sub(vec3_scale(vel, r), vec3_scale(pos, vec3_magnitude(vel) / mu)), 1.0 / mu);
INFO(" new e_vec = (" << e_vec_new.x << ", " << e_vec_new.y << ", " << e_vec_new.z << ")"); INFO(" new e_vec = (" << e_vec_new.x << ", " << e_vec_new.y << ", " << e_vec_new.z << ")");
// For zero-inclination orbit, omega is computed from the eccentricity vector // For zero-inclination orbit, omega should stay 0 (or 2π which is equivalent)
// (longitude of periapsis since ascending node is undefined) // If omega becomes π, that's a bug
// The key constraint is that omega should be in [0, 2π) bool omega_correct = (new_elements.argument_of_periapsis < M_PI / 2.0) ||
bool omega_in_range = (new_elements.argument_of_periapsis >= 0.0) && (new_elements.argument_of_periapsis > 1.5 * M_PI);
(new_elements.argument_of_periapsis < 2.0 * M_PI);
REQUIRE(omega_in_range); REQUIRE(omega_correct);
} }

2
tests/test_orbit_rendering.cpp

@ -1,7 +1,7 @@
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include "../src/physics.h" #include "../src/physics.h"
#include "../src/simulation.h" #include "../src/simulation.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include <cmath> #include <cmath>

2
tests/test_periapsis_burn.cpp

@ -2,7 +2,7 @@
#include <catch2/matchers/catch_matchers_floating_point.hpp> #include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "../src/physics.h" #include "../src/physics.h"
#include "../src/simulation.h" #include "../src/simulation.h"
#include "../src/orbital_objects.h" #include "../src/spacecraft.h"
#include "../src/maneuver.h" #include "../src/maneuver.h"
#include "../src/config_loader.h" #include "../src/config_loader.h"
#include "../src/orbital_mechanics.h" #include "../src/orbital_mechanics.h"

613
tests/test_rendezvous.cpp

@ -1,613 +0,0 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "../src/physics.h"
#include "../src/orbital_mechanics.h"
#include "../src/simulation.h"
#include "../src/orbital_objects.h"
#include "../src/rendezvous.h"
#include "../src/config_loader.h"
#include "../src/test_utilities.h"
#include <cmath>
#include <cstring>
using Catch::Matchers::WithinAbs;
// ============================================================================
// Helper Functions
// ============================================================================
static int find_spacecraft_by_name(SimulationState* sim, const char* name) {
for (int i = 0; i < sim->craft_count; i++) {
if (strcmp(sim->spacecraft[i].name, name) == 0) {
return i;
}
}
return -1;
}
// ============================================================================
// Test-only output helper
// ============================================================================
struct TestOutput {
char buf[32768];
int offset = 0;
void dump_state(SimulationState* sim, const char* label) {
int n = dump_simulation_state(sim, label, buf + offset, sizeof(buf) - offset);
if (n > 0) offset += n;
int target_idx = -1, chaser_idx = -1;
for (int i = 0; i < sim->craft_count; i++) {
if (strcmp(sim->spacecraft[i].name, "Target_Satellite") == 0)
target_idx = i;
if (strcmp(sim->spacecraft[i].name, "Chaser_Lower") == 0)
chaser_idx = i;
}
if (target_idx >= 0 && chaser_idx >= 0) {
Vec3 target_pos = sim->spacecraft[target_idx].local_position;
Vec3 chaser_pos = sim->spacecraft[chaser_idx].local_position;
double target_angle = atan2(target_pos.y, target_pos.x);
double chaser_angle = atan2(chaser_pos.y, chaser_pos.x);
double angular_sep = chaser_angle - target_angle;
while (angular_sep > M_PI) angular_sep -= 2.0 * M_PI;
while (angular_sep < -M_PI) angular_sep += 2.0 * M_PI;
Vec3 diff = vec3_sub(chaser_pos, target_pos);
double sep_mag = vec3_magnitude(diff);
n = snprintf(buf + offset, sizeof(buf) - offset,
" Angular separation (Chaser-Target): %.6f rad (%.4f deg)\n"
" Separation magnitude: %.2f m\n",
angular_sep, angular_sep * 180.0 / M_PI, sep_mag);
if (n > 0) offset += n;
}
}
};
TEST_CASE("Config loading for Hohmann transfer", "[rendezvous_hohmann][config]") {
const double TIME_STEP = 30.0;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
REQUIRE(sim->body_count == 1);
REQUIRE(std::string(sim->bodies[0].name) == "Earth");
REQUIRE(sim->craft_count == 3);
REQUIRE(std::string(sim->spacecraft[0].name) == "Target_Satellite");
REQUIRE(std::string(sim->spacecraft[1].name) == "Chaser_Lower");
REQUIRE(std::string(sim->spacecraft[2].name) == "Chaser_Higher");
REQUIRE(sim->spacecraft[0].parent_index == 0);
REQUIRE(sim->spacecraft[1].parent_index == 0);
REQUIRE(sim->spacecraft[2].parent_index == 0);
REQUIRE_THAT(sim->spacecraft[0].orbit.semi_major_axis, WithinAbs(6.771e6, 1.0));
REQUIRE_THAT(sim->spacecraft[1].orbit.semi_major_axis, WithinAbs(6.671e6, 1.0));
REQUIRE_THAT(sim->spacecraft[2].orbit.semi_major_axis, WithinAbs(6.871e6, 1.0));
REQUIRE_THAT(sim->spacecraft[0].orbit.true_anomaly, WithinAbs(0.0, 0.001));
REQUIRE_THAT(sim->spacecraft[1].orbit.true_anomaly, WithinAbs(4.71238898038469, 0.001));
REQUIRE_THAT(sim->spacecraft[2].orbit.true_anomaly, WithinAbs(1.5707963267948966, 0.001));
destroy_simulation(sim);
}
TEST_CASE("Calculate wait time for Hohmann transfer (lower to higher)", "[rendezvous_hohmann][phasing]") {
const double TIME_STEP = 30.0;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
int target_idx = find_spacecraft_by_name(sim, "Target_Satellite");
int chaser_lower_idx = find_spacecraft_by_name(sim, "Chaser_Lower");
REQUIRE(target_idx >= 0);
REQUIRE(chaser_lower_idx >= 0);
Spacecraft* target = &sim->spacecraft[target_idx];
Spacecraft* chaser = &sim->spacecraft[chaser_lower_idx];
CelestialBody* earth = &sim->bodies[0];
initialize_orbital_objects(sim);
double r1 = vec3_magnitude(chaser->local_position);
double r2 = vec3_magnitude(target->local_position);
SECTION("Zero angular separation - immediate transfer not possible") {
double angular_separation = 0.0;
double wait_time = calculate_wait_time_for_hohmann(r1, r2, angular_separation, earth->mass);
INFO("Wait time: " << wait_time << " s");
// At 0 separation, chaser is faster and needs to be behind target
// Wait time is negative (should have burned already)
REQUIRE_THAT(wait_time, WithinAbs(-1358.16, 1.0));
}
SECTION("Small angular separation") {
double angular_separation = 0.5;
double wait_time = calculate_wait_time_for_hohmann(r1, r2, angular_separation, earth->mass);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Wait time: " << wait_time << " s");
REQUIRE_THAT(wait_time, WithinAbs(-20909.0, 1.0));
}
SECTION("Large angular separation (near 2pi)") {
double angular_separation = 6.0;
double wait_time = calculate_wait_time_for_hohmann(r1, r2, angular_separation, earth->mass);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Wait time: " << wait_time << " s");
REQUIRE_THAT(wait_time, WithinAbs(9714.9, 1.0));
}
destroy_simulation(sim);
}
TEST_CASE("Calculate wait time for Hohmann transfer (higher to lower)", "[rendezvous_hohmann][phasing]") {
const double TIME_STEP = 30.0;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
int target_idx = find_spacecraft_by_name(sim, "Target_Satellite");
int chaser_higher_idx = find_spacecraft_by_name(sim, "Chaser_Higher");
REQUIRE(target_idx >= 0);
REQUIRE(chaser_higher_idx >= 0);
Spacecraft* target = &sim->spacecraft[target_idx];
Spacecraft* chaser = &sim->spacecraft[chaser_higher_idx];
CelestialBody* earth = &sim->bodies[0];
initialize_orbital_objects(sim);
double r1 = vec3_magnitude(chaser->local_position);
double r2 = vec3_magnitude(target->local_position);
SECTION("Zero angular separation - target must catch up") {
double angular_separation = 0.0;
double wait_time = calculate_wait_time_for_hohmann(r1, r2, angular_separation, earth->mass);
INFO("Wait time: " << wait_time << " s");
// Higher orbit case: chaser is slower, needs to be ahead of target
REQUIRE_THAT(wait_time, WithinAbs(-1414.46, 1.0));
}
SECTION("Small angular separation") {
double angular_separation = 0.3;
double wait_time = calculate_wait_time_for_hohmann(r1, r2, angular_separation, earth->mass);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Wait time: " << wait_time << " s");
REQUIRE_THAT(wait_time, WithinAbs(10757.3, 1.0));
}
destroy_simulation(sim);
}
TEST_CASE("Calculate required separation for Hohmann transfer", "[rendezvous_hohmann][phasing]") {
const double TIME_STEP = 30.0;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
int target_idx = find_spacecraft_by_name(sim, "Target_Satellite");
int chaser_lower_idx = find_spacecraft_by_name(sim, "Chaser_Lower");
int chaser_higher_idx = find_spacecraft_by_name(sim, "Chaser_Higher");
REQUIRE(target_idx >= 0);
REQUIRE(chaser_lower_idx >= 0);
REQUIRE(chaser_higher_idx >= 0);
Spacecraft* target = &sim->spacecraft[target_idx];
Spacecraft* chaser_lower = &sim->spacecraft[chaser_lower_idx];
Spacecraft* chaser_higher = &sim->spacecraft[chaser_higher_idx];
CelestialBody* earth = &sim->bodies[0];
initialize_orbital_objects(sim);
double r_lower = vec3_magnitude(chaser_lower->local_position);
double r_target = vec3_magnitude(target->local_position);
double r_higher = vec3_magnitude(chaser_higher->local_position);
SECTION("Lower to higher transfer") {
double required_separation = calculate_required_separation_for_hohmann(r_lower, r_target, earth->mass);
INFO("Required separation: " << required_separation << " rad");
INFO("Required separation (deg): " << required_separation * 180.0 / M_PI << " deg");
REQUIRE_THAT(required_separation, WithinAbs(-0.034734, 0.001)); // Chaser faster, needs to be behind
}
SECTION("Higher to lower transfer") {
double required_separation = calculate_required_separation_for_hohmann(r_higher, r_target, earth->mass);
INFO("Required separation: " << required_separation << " rad");
INFO("Required separation (deg): " << required_separation * 180.0 / M_PI << " deg");
REQUIRE_THAT(required_separation, WithinAbs(0.0348625, 0.001)); // Chaser slower, needs to be ahead
}
SECTION("Equal radii - no transfer needed") {
double required_separation = calculate_required_separation_for_hohmann(r_target, r_target, earth->mass);
INFO("Required separation: " << required_separation << " rad");
REQUIRE_THAT(required_separation, WithinAbs(0.0, 0.001));
}
destroy_simulation(sim);
}
// ============================================================================
// New Test Cases for Validation and Next Wait Time
// ============================================================================
TEST_CASE("Validate Hohmann transfer parameters", "[rendezvous_hohmann][validation]") {
const double TIME_STEP = 30.0;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
initialize_orbital_objects(sim);
Spacecraft* chaser_lower = &sim->spacecraft[1];
Spacecraft* chaser_higher = &sim->spacecraft[2];
Spacecraft* target = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[0];
double r_lower = vec3_magnitude(chaser_lower->local_position);
double r_higher = vec3_magnitude(chaser_higher->local_position);
double r_target = vec3_magnitude(target->local_position);
SECTION("Valid lower to higher transfer") {
bool valid = validate_hohmann_transfer_parameters(r_lower, r_target, earth->mass);
INFO("Valid lower to higher: " << (valid ? "true" : "false"));
REQUIRE(valid == true);
}
SECTION("Valid higher to lower transfer") {
bool valid = validate_hohmann_transfer_parameters(r_higher, r_target, earth->mass);
INFO("Valid higher to lower: " << (valid ? "true" : "false"));
REQUIRE(valid == true);
}
SECTION("Equal radii - invalid (no relative motion)") {
bool valid = validate_hohmann_transfer_parameters(r_target, r_target, earth->mass);
INFO("Equal radii valid: " << (valid ? "true" : "false"));
REQUIRE(valid == false);
}
SECTION("Negative initial radius - invalid") {
bool valid = validate_hohmann_transfer_parameters(-1000.0, r_target, earth->mass);
INFO("Negative radius valid: " << (valid ? "true" : "false"));
REQUIRE(valid == false);
}
SECTION("Negative target radius - invalid") {
bool valid = validate_hohmann_transfer_parameters(r_lower, -1000.0, earth->mass);
INFO("Negative target radius valid: " << (valid ? "true" : "false"));
REQUIRE(valid == false);
}
SECTION("Zero central mass - invalid") {
bool valid = validate_hohmann_transfer_parameters(r_lower, r_target, 0.0);
INFO("Zero mass valid: " << (valid ? "true" : "false"));
REQUIRE(valid == false);
}
SECTION("Negative central mass - invalid") {
bool valid = validate_hohmann_transfer_parameters(r_lower, r_target, -1.0);
INFO("Negative mass valid: " << (valid ? "true" : "false"));
REQUIRE(valid == false);
}
destroy_simulation(sim);
}
TEST_CASE("Calculate relative orbit period", "[rendezvous_hohmann][phasing]") {
const double TIME_STEP = 30.0;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
initialize_orbital_objects(sim);
Spacecraft* chaser_lower = &sim->spacecraft[1];
Spacecraft* chaser_higher = &sim->spacecraft[2];
Spacecraft* target = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[0];
double r_lower = vec3_magnitude(chaser_lower->local_position);
double r_higher = vec3_magnitude(chaser_higher->local_position);
double r_target = vec3_magnitude(target->local_position);
SECTION("Lower to higher transfer") {
double rel_period = calculate_relative_orbit_period(r_lower, r_target, earth->mass);
INFO("Relative orbit period: " << rel_period << " s");
INFO("Relative orbit period: " << rel_period / 3600.0 << " h");
REQUIRE(rel_period > 0.0);
REQUIRE_THAT(rel_period, WithinAbs(245683.24, 1.0));
}
SECTION("Higher to lower transfer") {
double rel_period = calculate_relative_orbit_period(r_higher, r_target, earth->mass);
INFO("Relative orbit period: " << rel_period << " s");
INFO("Relative orbit period: " << rel_period / 3600.0 << " h");
REQUIRE(rel_period > 0.0);
REQUIRE_THAT(rel_period, WithinAbs(254924.71, 1.0));
}
SECTION("Very small radius difference - long period") {
double r1 = 6770000.0;
double r2 = 6771000.0;
double rel_period = calculate_relative_orbit_period(r1, r2, earth->mass);
INFO("Relative orbit period: " << rel_period << " s");
INFO("Relative orbit period: " << rel_period / 3600.0 / 24.0 << " days");
REQUIRE(rel_period > 0.0);
REQUIRE_THAT(rel_period, WithinAbs(25025208.27, 1.0));
}
SECTION("Large radius difference - short period") {
double r1 = 6571000.0;
double r2 = 7071000.0;
double rel_period = calculate_relative_orbit_period(r1, r2, earth->mass);
INFO("Relative orbit period: " << rel_period << " s");
INFO("Relative orbit period: " << rel_period / 3600.0 << " h");
REQUIRE(rel_period > 0.0);
REQUIRE_THAT(rel_period, WithinAbs(50889.08, 1.0));
}
destroy_simulation(sim);
}
TEST_CASE("Calculate next valid wait time for Hohmann transfer", "[rendezvous_hohmann][phasing]") {
const double TIME_STEP = 30.0;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
initialize_orbital_objects(sim);
Spacecraft* chaser_lower = &sim->spacecraft[1];
Spacecraft* chaser_higher = &sim->spacecraft[2];
Spacecraft* target = &sim->spacecraft[0];
CelestialBody* earth = &sim->bodies[0];
double r_lower = vec3_magnitude(chaser_lower->local_position);
double r_higher = vec3_magnitude(chaser_higher->local_position);
double r_target = vec3_magnitude(target->local_position);
SECTION("Lower to higher - positive wait time") {
double angular_separation = 0.0;
double next_wait = calculate_next_hohmann_wait_time(r_lower, r_target, angular_separation, earth->mass, TIME_STEP);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Next wait time: " << next_wait << " s");
INFO("Next wait time: " << next_wait / 3600.0 << " h");
REQUIRE(next_wait >= TIME_STEP);
REQUIRE_THAT(next_wait, WithinAbs(244325.08, 1.0));
}
SECTION("Higher to lower - negative wait time becomes next cycle") {
double angular_separation = 0.0;
double next_wait = calculate_next_hohmann_wait_time(r_higher, r_target, angular_separation, earth->mass, TIME_STEP);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Next wait time: " << next_wait << " s");
INFO("Next wait time: " << next_wait / 3600.0 << " h");
REQUIRE(next_wait >= TIME_STEP);
REQUIRE_THAT(next_wait, WithinAbs(253510.25, 1.0));
}
SECTION("Higher to lower - positive wait time (no adjustment needed)") {
double angular_separation = 0.3;
double next_wait = calculate_next_hohmann_wait_time(r_higher, r_target, angular_separation, earth->mass, TIME_STEP);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Next wait time: " << next_wait << " s");
REQUIRE(next_wait >= TIME_STEP);
REQUIRE_THAT(next_wait, WithinAbs(10757.30, 1.0));
}
SECTION("At required separation - forced to next cycle by min_wait") {
double req_sep = calculate_required_separation_for_hohmann(r_lower, r_target, earth->mass);
double next_wait = calculate_next_hohmann_wait_time(r_lower, r_target, req_sep, earth->mass, TIME_STEP);
INFO("Required separation: " << req_sep << " rad");
INFO("Next wait time: " << next_wait << " s");
INFO("Next wait time: " << next_wait / 3600.0 << " h");
REQUIRE(next_wait >= TIME_STEP);
REQUIRE_THAT(next_wait, WithinAbs(245683.24, 1.0));
}
SECTION("Force skip to next cycle with large min_wait") {
double angular_separation = 0.5;
double min_wait = 30000.0;
double next_wait = calculate_next_hohmann_wait_time(r_lower, r_target, angular_separation, earth->mass, min_wait);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Min wait: " << min_wait << " s");
INFO("Next wait time: " << next_wait << " s");
REQUIRE(next_wait >= min_wait);
REQUIRE_THAT(next_wait, WithinAbs(224774.23, 1.0));
}
SECTION("Zero min_wait allows immediate transfer") {
double angular_separation = 0.3;
double next_wait = calculate_next_hohmann_wait_time(r_higher, r_target, angular_separation, earth->mass, 0.0);
INFO("Angular separation: " << angular_separation << " rad");
INFO("Next wait time: " << next_wait << " s");
REQUIRE(next_wait >= 0.0);
REQUIRE_THAT(next_wait, WithinAbs(10757.30, 1.0));
}
destroy_simulation(sim);
}
SCENARIO("Hohmann transfer rendezvous with validation", "[rendezvous_hohmann][integration]") {
const double TIME_STEP = 0.1;
SimulationState* sim = create_simulation(3, 5, 10, TIME_STEP);
REQUIRE(load_system_config(sim, "tests/test_rendezvous.toml"));
int target_idx = find_spacecraft_by_name(sim, "Target_Satellite");
int chaser_lower_idx = find_spacecraft_by_name(sim, "Chaser_Lower");
REQUIRE(target_idx >= 0);
REQUIRE(chaser_lower_idx >= 0);
Spacecraft* target = &sim->spacecraft[target_idx];
Spacecraft* chaser = &sim->spacecraft[chaser_lower_idx];
CelestialBody* earth = &sim->bodies[0];
initialize_orbital_objects(sim);
SECTION("Execute validated Hohmann transfer and verify rendezvous") {
double r1 = vec3_magnitude(chaser->local_position);
double r2 = vec3_magnitude(target->local_position);
INFO("\n=== HOHMANN TRANSFER RENDZVOUS TEST ===");
INFO("Initial state:");
INFO(" Chaser: r=" << r1 << " m, nu=" << chaser->orbit.true_anomaly << " rad");
INFO(" Target: r=" << r2 << " m, nu=" << target->orbit.true_anomaly << " rad");
// Calculate angular separation
double angular_separation = chaser->orbit.true_anomaly - target->orbit.true_anomaly;
while (angular_separation > M_PI) angular_separation -= 2.0 * M_PI;
while (angular_separation < -M_PI) angular_separation += 2.0 * M_PI;
INFO(" Angular separation: " << angular_separation << " rad");
// Validate parameters before proceeding
REQUIRE(validate_hohmann_transfer_parameters(r1, r2, earth->mass));
INFO(" Validation passed");
// Calculate Hohmann transfer parameters
HohmannTransfer hohmann = calculate_hohmann_transfer(r1, r2, earth->mass);
double wait_time = calculate_next_hohmann_wait_time(r1, r2, angular_separation, earth->mass, TIME_STEP);
double arrival_time = wait_time + hohmann.transfer_time;
INFO("\nTransfer parameters:");
INFO(" Transfer time: " << hohmann.transfer_time << " s");
INFO(" DV1 (departure): " << hohmann.dv1 << " m/s");
INFO(" DV2 (arrival): " << hohmann.dv2 << " m/s");
INFO(" Wait time: " << wait_time << " s");
INFO(" Arrival time: " << arrival_time << " s");
// Create departure maneuver
Maneuver departure = create_maneuver(
"Departure_Burn",
chaser_lower_idx,
BURN_PROGRADE,
hohmann.dv1,
TRIGGER_TIME,
wait_time
);
int dep_idx = add_maneuver_to_simulation(sim, &departure);
REQUIRE(dep_idx >= 0);
// Create arrival maneuver
Maneuver arrival = create_maneuver(
"Circularization_Burn",
chaser_lower_idx,
BURN_PROGRADE,
hohmann.dv2,
TRIGGER_TIME,
arrival_time
);
int arr_idx = add_maneuver_to_simulation(sim, &arrival);
REQUIRE(arr_idx >= 0);
// Compute expected step count from analytical times
const double expected_duration = wait_time + hohmann.transfer_time;
const int expected_steps = static_cast<int>(expected_duration / sim->dt);
// Safety limit: 1 year — prevents infinite loops if calculation is wrong
const double SAFETY_LIMIT_SECONDS = 3600.0 * 24.0 * 365.0;
TestOutput out;
bool transfer_complete = false;
for (int i = 0; i < expected_steps + 1000; i++) {
update_simulation(sim);
if (sim->time > SAFETY_LIMIT_SECONDS) {
INFO("\n=== SAFETY LIMIT REACHED ===");
INFO("sim->time: " << sim->time << " s (limit: " << SAFETY_LIMIT_SECONDS << " s)");
break;
}
if (i == 0) out.dump_state(sim, "T=0 (initial)");
if (i == static_cast<int>(wait_time / sim->dt)) out.dump_state(sim, "JUST BEFORE DEPARTURE");
if (i == static_cast<int>(wait_time / sim->dt) + 1) out.dump_state(sim, "AFTER DEPARTURE BURN");
if (i == static_cast<int>(arrival_time / sim->dt) - 1) out.dump_state(sim, "JUST BEFORE ARRIVAL BURN");
if (sim->maneuvers[arr_idx].executed && !transfer_complete) {
out.dump_state(sim, "AFTER ARRIVAL BURN");
transfer_complete = true;
break;
}
}
INFO(out.buf);
// Verify rendezvous quality
double final_radius = vec3_magnitude(chaser->local_position);
double radius_error = fabs(final_radius - r2);
Vec3 chaser_vel = chaser->local_velocity;
Vec3 target_vel = target->local_velocity;
double chaser_speed = vec3_magnitude(chaser_vel);
double target_speed = vec3_magnitude(target_vel);
Vec3 separation = vec3_sub(chaser->local_position, target->local_position);
double separation_distance = vec3_magnitude(separation);
double relative_velocity = vec3_magnitude(vec3_sub(chaser_vel, target_vel));
INFO("\n=== FINAL STATE ===");
INFO("Final time: " << sim->time << " s");
INFO("Departure executed: " << sim->maneuvers[dep_idx].executed);
INFO("Arrival executed: " << sim->maneuvers[arr_idx].executed);
INFO("\nVerification:");
INFO(" Radius error: " << radius_error << " m");
INFO(" Chaser eccentricity: " << chaser->orbit.eccentricity);
INFO(" Chaser speed: " << chaser_speed << " m/s");
INFO(" Target speed: " << target_speed << " m/s");
INFO(" Separation: " << separation_distance << " m");
INFO(" Relative velocity: " << relative_velocity << " m/s");
// Verify maneuvers executed
REQUIRE(sim->maneuvers[dep_idx].executed);
REQUIRE(sim->maneuvers[arr_idx].executed);
REQUIRE_THAT(radius_error, WithinAbs(10.0, 10.0));
REQUIRE_THAT(chaser->orbit.eccentricity, WithinAbs(0.0, 0.001));
REQUIRE_THAT(chaser_speed, WithinAbs(target_speed, 1.0));
REQUIRE_THAT(separation_distance, WithinAbs(0.0, 100.0));
}
destroy_simulation(sim);
}

70
tests/test_rendezvous.toml

@ -1,70 +0,0 @@
# Test Configuration: Hohmann Transfer Between Coplanar Circular Orbits
# Three spacecraft in circular, coplanar LEO orbits around Earth
# Tests Hohmann transfer phasing calculations in both directions
#
# Configuration:
# - Target_Satellite: middle orbit (reference)
# - Chaser_Lower: lower orbit, will transfer UP to target
# - Chaser_Higher: higher orbit, will transfer DOWN to target
[[bodies]]
name = "Earth"
mass = 5.972e24
radius = 6.371e6
parent_index = -1
color = { r = 0.0, g = 0.5, b = 1.0 }
orbit = {
semi_major_axis = 0.0,
eccentricity = 0.0,
true_anomaly = 0.0
}
# ========== TARGET SPACECRAFT ==========
# Circular LEO orbit at 400 km altitude
# Reference orbit for Hohmann transfers
[[spacecraft]]
name = "Target_Satellite"
mass = 500.0
parent_index = 0
orbit = {
semi_major_axis = 6.771e6,
eccentricity = 0.0,
true_anomaly = 0.0,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
# ========== CHASER LOWER ==========
# Circular LEO orbit at 300 km altitude (lower than target)
# Will perform Hohmann transfer UP to target orbit
# Starts 90 degrees behind target to test phasing
[[spacecraft]]
name = "Chaser_Lower"
mass = 500.0
parent_index = 0
orbit = {
semi_major_axis = 6.671e6,
eccentricity = 0.0,
true_anomaly = 4.71238898038469,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
# ========== CHASER HIGHER ==========
# Circular LEO orbit at 500 km altitude (higher than target)
# Will perform Hohmann transfer DOWN to target orbit
# Starts 90 degrees ahead of target to test phasing
[[spacecraft]]
name = "Chaser_Higher"
mass = 500.0
parent_index = 0
orbit = {
semi_major_axis = 6.871e6,
eccentricity = 0.0,
true_anomaly = 1.5707963267948966,
inclination = 0.0,
longitude_of_ascending_node = 0.0,
argument_of_periapsis = 0.0
}
Loading…
Cancel
Save