@ -36,27 +36,20 @@ git stash pop # Apply debug changes
### Step 1.1: Add spacecraft to `tests/configs/earth_mars_simple.toml` ### Step 1.1: Add spacecraft to `tests/configs/earth_mars_simple.toml`
Append to config file: Add Spacecraft body to config with placeholder position/velocity (set at runtime by `initialize_spacecraft_leo()` ).
```toml
[[bodies]] **Implementation:** See `tests/configs/earth_mars_simple.toml` for full config
name = "Spacecraft"
mass = 1.0
radius = 1000.0
# Position and velocity will be initialized at runtime for LEO orbit
position = { x = 0.0, y = 0.0, z = 0.0 }
velocity = { x = 0.0, y = 0.0, z = 0.0 }
parent_index = 1 # Earth
color = { r = 1.0, g = 0.0, b = 0.5 }
eccentricity = 0.0
# Semi-major axis will be: Earth radius + 200km
semi_major_axis = 6.571e6 # Placeholder, will be set during initialization
```
**Note**: Position/velocity are placeholders; will be calculated by `initialize_spacecraft_leo()` at runtime. **Key parameters:**
- mass = 1.0 kg (test particle)
- radius = 1000.0 m
- parent_index = 1 (Earth)
- color = magenta (r=1.0, g=0.0, b=0.5)
- position/velocity: Placeholders (0,0,0)
**TODO**: Future config file format should support: **TODO**: Future config format should support:
- Earth-relative position (e.g., `{ altitude_km = 200.0 }` ) - Earth-relative position: `{ altitude_km = 200.0 }`
- Earth-relative velocity (e.g., `{ orbit_type = "circular" }` ) - Earth-relative orbit: `{ orbit_type = "circular" }`
- More intuitive spacecraft mission parameters - More intuitive spacecraft mission parameters
--- ---
@ -65,334 +58,86 @@ semi_major_axis = 6.571e6 # Placeholder, will be set during initialization
### Step 2.1: Add function declarations to `src/mission_planning.h` ### Step 2.1: Add function declarations to `src/mission_planning.h`
```cpp **Implementation:** See `src/mission_planning.h`
// Initialize spacecraft in circular LEO around parent body
void initialize_spacecraft_leo(CelestialBody* spacecraft, CelestialBody* parent,
double altitude_m);
// Apply patched conics impulse burn for Hohmann transfer **Functions:**
void apply_transfer_burn(SimulationState* sim, int spacecraft_idx, - `initialize_spacecraft_leo()` - Initialize spacecraft in circular LEO around parent body
int departure_idx, TransferParameters* params); - `apply_transfer_burn()` - Apply patched conics impulse burn for Hohmann transfer
- `calculate_phase_angle()` - Calculate current phase angle between two bodies (in degrees)
// Helper: Calculate current phase angle between two bodies (in degrees)
double calculate_phase_angle(SimulationState* sim, int departure_idx, int arrival_idx);
```
### Step 2.2: Implement `initialize_spacecraft_leo()` in `src/mission_planning.cpp` ### Step 2.2: Implement `initialize_spacecraft_leo()` in `src/mission_planning.cpp`
**Algorithm**: **Implementation:** `src/mission_planning.cpp:20-56`
```cpp
void initialize_spacecraft_leo(CelestialBody* spacecraft, CelestialBody* parent,
double altitude_m) {
// Calculate orbital radius (distance from Earth center)
double orbital_radius = parent->radius + altitude_m;
// Position spacecraft radially outward from Earth-Sun line
// Get vector from Sun to Earth
Vec3 sun_to_earth = vec3_sub(parent->position,
(Vec3){0.0, 0.0, 0.0}); // Sun at origin
Vec3 direction = vec3_normalize(sun_to_earth);
// Position: Earth position + offset radially outward **Algorithm:**
Vec3 offset = vec3_scale(direction, orbital_radius); - Calculate orbital radius = parent radius + altitude
spacecraft->position = vec3_add(parent->position, offset); - Position spacecraft radially outward from Sun (any angular position acceptable)
- Calculate circular LEO velocity: v = sqrt(G * M_parent / r)
- Set prograde orientation (tangential to Earth-Sun line)
- Set both local and global coordinates correctly
// Initialize local coordinates (relative to parent) **Key Points:**
spacecraft->local_position = offset; - LEO orbit is circular at 200km altitude (~7,788 m/s)
spacecraft->local_velocity = (Vec3){0.0, 0.0, 0.0}; // Will be set below - Spacecraft velocity = Earth velocity + LEO velocity
- Local velocity = LEO velocity only (relative to Earth)
// Calculate circular LEO velocity magnitude
double v_leo = sqrt(G * parent->mass / orbital_radius);
// Direction: tangential to Earth-Sun line (prograde)
// If sun_to_earth = (x, y, 0), then tangent = (-y, x, 0)
Vec3 leo_tangent = (Vec3){-direction.y, direction.x, 0.0};
Vec3 leo_velocity = vec3_scale(leo_tangent, v_leo);
// Spacecraft velocity = Earth velocity + LEO velocity
spacecraft->velocity = vec3_add(parent->velocity, leo_velocity);
// Local velocity relative to Earth = LEO velocity only
spacecraft->local_velocity = leo_velocity;
// Update semi-major axis for reference
spacecraft->semi_major_axis = orbital_radius;
// SOI will be calculated by config loader
}
```
**Key Points**:
- Spacecraft positioned radially outward from Sun (any position is acceptable)
- LEO orbit is circular at 200km altitude
- Prograde orientation (same direction as Earth's orbital velocity)
- Both local and global coordinates set correctly
### Step 2.3: Implement `calculate_phase_angle()` in `src/mission_planning.cpp` ### Step 2.3: Implement `calculate_phase_angle()` in `src/mission_planning.cpp`
**Algorithm**: **Implementation:** `src/mission_planning.cpp:58-78`
```cpp
double calculate_phase_angle(SimulationState* sim, int departure_idx, int arrival_idx) { **Algorithm:**
CelestialBody* departure = &sim->bodies[departure_idx]; - Calculate angular positions of departure and arrival bodies relative to Sun
CelestialBody* arrival = &sim->bodies[arrival_idx]; - Compute phase difference: θ_arrival - θ_departure
CelestialBody* sun = &sim->bodies[0]; // Assume Sun at index 0 - Normalize to [0°, 360°) range
- Return phase angle in degrees
// Calculate angular positions relative to Sun
double theta_depart = calculate_angular_position(departure, sun);
double theta_arrival = calculate_angular_position(arrival, sun);
// Calculate phase difference
double phase_rad = theta_arrival - theta_depart;
// Normalize to [0, 2π)
while (phase_rad < 0.0 ) {
phase_rad += 2.0 * M_PI;
}
while (phase_rad >= 2.0 * M_PI) {
phase_rad -= 2.0 * M_PI;
}
// Convert to degrees
return phase_rad * 180.0 / M_PI;
}
```
### Step 2.4: Implement `apply_transfer_burn()` in `src/mission_planning.cpp` ### Step 2.4: Implement `apply_transfer_burn()` in `src/mission_planning.cpp`
**Algorithm (Patched Conics Approach)**: **Implementation:** `src/mission_planning.cpp:80-116`
```cpp
void apply_transfer_burn(SimulationState* sim, int spacecraft_idx, **Algorithm (Patched Conics Approach):**
int departure_idx, TransferParameters* params) { - Calculate required heliocentric transfer velocity magnitude from params
CelestialBody* spacecraft = &sim->bodies[spacecraft_idx]; - Determine prograde direction (tangential to departure-Sun line)
CelestialBody* departure = &sim->bodies[departure_idx]; - Compute delta-v: Δv = v_transfer - v_current (vector subtraction)
CelestialBody* sun = &sim->bodies[0]; // Assume Sun at index 0 - Apply impulse to spacecraft velocity
- Update local velocity relative to departure body
// Calculate required heliocentric transfer velocity - Print burn information for debugging
// v_transfer = params->departure_velocity
// Direction: prograde (tangential to Earth-Sun line)
Vec3 sun_to_earth = vec3_sub(departure->position, sun->position);
Vec3 sun_to_earth_norm = vec3_normalize(sun_to_earth);
// Tangent direction (prograde): (-y, x, 0)
Vec3 transfer_dir = (Vec3){-sun_to_earth_norm.y, sun_to_earth_norm.x, 0.0};
Vec3 v_transfer_helio = vec3_scale(transfer_dir, params->departure_velocity);
// Current heliocentric velocity
Vec3 current_helio = spacecraft->velocity;
// Calculate total Δv to apply
Vec3 delta_v = vec3_sub(v_transfer_helio, current_helio);
// Apply impulse burn
spacecraft->velocity = vec3_add(spacecraft->velocity, delta_v);
// Update local velocity
spacecraft->local_velocity = vec3_sub(spacecraft->velocity, departure->velocity);
// Print burn information
printf("Transfer burn applied:\n");
printf(" Current heliocentric velocity: (%.2f, %.2f, %.2f) m/s\n",
current_helio.x, current_helio.y, current_helio.z);
printf(" Target heliocentric velocity: (%.2f, %.2f, %.2f) m/s\n",
v_transfer_helio.x, v_transfer_helio.y, v_transfer_helio.z);
printf(" Delta-v: (%.2f, %.2f, %.2f) m/s\n",
delta_v.x, delta_v.y, delta_v.z);
printf(" Delta-v magnitude: %.2f m/s (%.3f km/s)\n",
vec3_magnitude(delta_v), vec3_magnitude(delta_v) / 1000.0);
}
```
**Note**: This is a simplified single-impulse approach. A true patched conics calculation would: **Note:** Simplified single-impulse approximation. True patched conics would:
1. Calculate Δv to reach SOI boundary (escape trajectory) 1. Calculate Δv to reach SOI boundary
2. Calculate velocity at SOI boundary 2. Calculate velocity at SOI boundary
3. Add transfer Δv at SOI boundary 3. Add transfer Δv at SOI boundary
4. Combine into equivalent single impulse 4. Combine into equivalent single impulse
For initial implementation, we'll use single impulse as approximation.
--- ---
## Phase 3: Comprehensive Test Case ## Phase 3: Comprehensive Test Case
### Step 3.1: Create new test in `tests/test_hohmann_transfer.cpp` ### Step 3.1: Create new test in `tests/test_hohmann_transfer.cpp`
```cpp **Implementation:** `tests/test_hohmann_transfer.cpp` - See file for full test
TEST_CASE("Earth → Mars Hohmann Transfer with LEO Spacecraft", "[mission][hohmann][config][integration]") {
const double TIME_STEP = 60.0; **Test:** "Earth → Mars Hohmann Transfer with LEO Spacecraft"
const double SECONDS_PER_DAY = 86400.0;
const double LEO_ALTITUDE_M = 200000.0; // 200 km **Test Structure:**
1. Load config with 4 bodies (Sun, Earth, Mars, Spacecraft)
// 1. Load config with LEO spacecraft 2. Initialize spacecraft in 200km LEO around Earth
SimulationState* sim = create_simulation(4, TIME_STEP); 3. Verify LEO orbit stability (parent, position, velocity, energy)
REQUIRE(load_system_config(sim, "tests/configs/earth_mars_simple.toml")); 4. Calculate Hohmann transfer parameters
5. Wait for Earth-Mars launch window (within 1° tolerance)
const int SUN_IDX = 0; 6. Verify phase angle accuracy
const int EARTH_IDX = 1; 7. Apply impulse burn for transfer
const int MARS_IDX = 2; 8. Verify post-burn energy >= 0 (escape trajectory)
const int CRAFT_IDX = 3; 9. Simulate transfer for 110% of expected duration
10. Track SOI transitions (Earth→Sun→Mars)
// Verify spacecraft loaded 11. Verify final parent and energy conservation (< 5 % drift )
REQUIRE(sim->body_count == 4); 12. If Mars SOI entry, verify distance (< 2 × SOI )
REQUIRE(strcmp(sim->bodies[CRAFT_IDX].name, "Spacecraft") == 0);
**Key Assertions:**
// 2. Initialize spacecraft LEO orbit - Config loading: 4 bodies loaded, spacecraft present
initialize_spacecraft_leo(& sim->bodies[CRAFT_IDX], & sim->bodies[EARTH_IDX], - LEO stability: parent=Earth, position < 1km error , velocity < 10m / s error , energy < 0
LEO_ALTITUDE_M); - Launch window: opens in ~94 days, phase error < 1 °
- Transfer: post-burn energy >= 0, Earth→Sun SOI transition, energy conservation
INFO("Spacecraft initialized at %.2f km altitude", LEO_ALTITUDE_M / 1000.0);
INFO("Spacecraft parent: %d (Earth)", sim->bodies[CRAFT_IDX].parent_index);
// 3. Verify initial LEO orbit is stable
REQUIRE(sim->bodies[CRAFT_IDX].parent_index == EARTH_IDX);
double dist_to_earth = vec3_distance(sim->bodies[CRAFT_IDX].position,
sim->bodies[EARTH_IDX].position);
double expected_radius = sim->bodies[EARTH_IDX].radius + LEO_ALTITUDE_M;
REQUIRE(fabs(dist_to_earth - expected_radius) < 1000.0 ) ; / / Within 1 km
// Verify LEO velocity magnitude
double leo_velocity_mag = sqrt(G * sim->bodies[EARTH_IDX].mass / dist_to_earth);
double v_leo_relative = vec3_magnitude(sim->bodies[CRAFT_IDX].local_velocity);
INFO("Expected LEO velocity: %.2f m/s", leo_velocity_mag);
INFO("Actual LEO velocity: %.2f m/s", v_leo_relative);
REQUIRE(fabs(v_leo_relative - leo_velocity_mag) < 10.0 ) ; / / Within 10 m / s
// Verify negative total energy (bound to Earth)
OrbitalMetrics leo_metrics = calculate_orbital_metrics(& sim->bodies[CRAFT_IDX],
&sim->bodies[EARTH_IDX]);
INFO("LEO total energy: %.2e J", leo_metrics.total_energy);
REQUIRE(leo_metrics.total_energy < 0.0 ) ;
// 4. Calculate Hohmann transfer parameters
double r_earth = vec3_distance(sim->bodies[EARTH_IDX].position,
sim->bodies[SUN_IDX].position);
double r_mars = vec3_distance(sim->bodies[MARS_IDX].position,
sim->bodies[SUN_IDX].position);
TransferParameters params = calculate_hohmann_transfer(r_earth, r_mars,
sim->bodies[SUN_IDX].mass);
INFO("Transfer time: %.2f days", params.transfer_time / SECONDS_PER_DAY);
INFO("Required phase angle: %.3f degrees", params.phase_angle_deg);
INFO("Delta-v injection: %.3f km/s", params.delta_v_injection / 1000.0);
// 5. Wait for Earth-Mars launch window
double wait_start_time = sim->time;
wait_for_launch_window(sim, EARTH_IDX, MARS_IDX, params.phase_angle_deg, 1.0);
double wait_duration = sim->time - wait_start_time;
INFO("Launch window opened after %.2f days", wait_duration / SECONDS_PER_DAY);
// 6. Verify launch window accuracy (within 1°)
double current_phase = calculate_phase_angle(sim, EARTH_IDX, MARS_IDX);
double phase_error = fabs(current_phase - params.phase_angle_deg);
if (phase_error > 180.0) phase_error = fabs(phase_error - 360.0);
INFO("Current phase angle: %.3f degrees", current_phase);
INFO("Required phase angle: %.3f degrees", params.phase_angle_deg);
INFO("Phase angle error: %.3f degrees", phase_error);
REQUIRE(phase_error < 1.0 ) ;
// 7. Apply impulse burn for transfer
double pre_burn_time = sim->time;
OrbitalMetrics pre_burn_metrics = calculate_orbital_metrics(& sim->bodies[CRAFT_IDX],
&sim->bodies[SUN_IDX]);
apply_transfer_burn(sim, CRAFT_IDX, EARTH_IDX, ¶ms);
OrbitalMetrics post_burn_metrics = calculate_orbital_metrics(& sim->bodies[CRAFT_IDX],
&sim->bodies[SUN_IDX]);
INFO("Pre-burn heliocentric energy: %.2e J", pre_burn_metrics.total_energy);
INFO("Post-burn heliocentric energy: %.2e J", post_burn_metrics.total_energy);
INFO("Energy added: %.2e J",
post_burn_metrics.total_energy - pre_burn_metrics.total_energy);
// Verify spacecraft is now in escape trajectory (positive or zero energy)
REQUIRE(post_burn_metrics.total_energy >= 0.0);
// 8. Track SOI transitions during transfer
int earth_soi_exit_step = 0;
int sun_soi_enter_step = 0;
int mars_soi_enter_step = 0;
double transfer_duration = params.transfer_time * 1.1;
int max_steps = (int)(transfer_duration / sim->dt);
INFO("Simulating for %.2f days (%d steps)",
transfer_duration / SECONDS_PER_DAY, max_steps);
for (int step = 0; step < max_steps ; step + + ) {
update_simulation(sim);
// Track Earth SOI exit
if (earth_soi_exit_step == 0 & &
sim->bodies[CRAFT_IDX].parent_index != EARTH_IDX) {
earth_soi_exit_step = step;
INFO("Earth SOI exit at step %d (t = %.2f days)",
step, sim->time / SECONDS_PER_DAY);
}
// Track Sun SOI entry (after leaving Earth)
if (earth_soi_exit_step > 0 & & sun_soi_enter_step == 0 & &
sim->bodies[CRAFT_IDX].parent_index == SUN_IDX) {
sun_soi_enter_step = step;
INFO("Sun SOI entry at step %d (t = %.2f days)",
step, sim->time / SECONDS_PER_DAY);
}
// Track Mars SOI entry
if (mars_soi_enter_step == 0 & &
sim->bodies[CRAFT_IDX].parent_index == MARS_IDX) {
mars_soi_enter_step = step;
INFO("Mars SOI entry at step %d (t = %.2f days)",
step, sim->time / SECONDS_PER_DAY);
}
}
// 9. Verify Earth → Sun transition occurred
INFO("Earth SOI exit step: %d", earth_soi_exit_step);
INFO("Sun SOI entry step: %d", sun_soi_enter_step);
REQUIRE(earth_soi_exit_step > 0);
REQUIRE(sun_soi_enter_step > 0);
// Final parent should be Sun or Mars
int final_parent = sim->bodies[CRAFT_IDX].parent_index;
REQUIRE(final_parent == SUN_IDX || final_parent == MARS_IDX);
INFO("Final parent: %d (%s)", final_parent,
final_parent == SUN_IDX ? "Sun" : "Mars");
// 10. Verify spacecraft followed transfer orbit (energy conservation)
OrbitalMetrics final_metrics = calculate_orbital_metrics(& sim->bodies[CRAFT_IDX],
&sim->bodies[SUN_IDX]);
double energy_drift = fabs(final_metrics.total_energy - post_burn_metrics.total_energy);
if (post_burn_metrics.total_energy != 0.0) {
energy_drift /= fabs(post_burn_metrics.total_energy);
}
INFO("Final orbital radius: %.2f AU",
final_metrics.orbital_radius / 1.496e11);
INFO("Final energy: %.2e J", final_metrics.total_energy);
INFO("Expected energy: %.2e J", post_burn_metrics.total_energy);
INFO("Energy drift: %.2f%%", energy_drift * 100.0);
REQUIRE(energy_drift < 0.05 ) ; / / < 5 % energy conservation
// 11. If Mars SOI entry occurred, verify distance
if (mars_soi_enter_step > 0) {
double dist_to_mars = vec3_distance(sim->bodies[CRAFT_IDX].position,
sim->bodies[MARS_IDX].position);
INFO("Distance to Mars: %.2f km", dist_to_mars / 1000.0);
INFO("Mars SOI radius: %.2f km", sim->bodies[MARS_IDX].soi_radius / 1000.0);
REQUIRE(dist_to_mars < 2.0 * sim- > bodies[MARS_IDX].soi_radius);
} else {
INFO("Spacecraft did not enter Mars SOI within simulation time");
INFO("This may be due to phase angle or timing inaccuracies");
}
destroy_simulation(sim);
}
```
--- ---
@ -578,50 +323,20 @@ Spacecraft may not enter Mars SOI due to:
## Test Configuration Reference ## Test Configuration Reference
### earth_mars_simple.toml ### earth_mars_simple.toml
```toml **Implementation:** `tests/configs/earth_mars_simple.toml`
[[bodies]]
name = "Sun" **Bodies:**
mass = 1.989e30 - Sun (index 0): Root body, 1.989e30 kg
radius = 6.96e8 - Earth (index 1): 5.972e24 kg, 1.496e11 m from Sun
position = { x = 0.0, y = 0.0, z = 0.0 } - Mars (index 2): 6.39e23 kg, 2.279e11 m from Sun
parent_index = -1 - Spacecraft (index 3): 1.0 kg, parent=Earth (position/velocity set at runtime)
color = { r = 1.0, g = 1.0, b = 0.0 }
eccentricity = 0.0 **Spacecraft parameters:**
semi_major_axis = 0.0 - mass = 1.0 kg
- radius = 1000.0 m
[[bodies]] - parent_index = 1 (Earth)
name = "Earth" - color = magenta (r=1.0, g=0.0, b=0.5)
mass = 5.972e24 - position/velocity: Placeholders (0,0,0) - set by `initialize_spacecraft_leo()`
radius = 6.371e6
position = { x = 1.496e11, y = 0.0, z = 0.0 }
parent_index = 0
color = { r = 0.0, g = 0.5, b = 1.0 }
eccentricity = 0.0
semi_major_axis = 1.496e11
[[bodies]]
name = "Mars"
mass = 6.39e23
radius = 3.3895e6
position = { x = 2.279e11, y = 0.0, z = 0.0 }
parent_index = 0
color = { r = 0.8, g = 0.3, b = 0.1 }
eccentricity = 0.0
semi_major_axis = 2.279e11
[[bodies]]
name = "Spacecraft"
mass = 1.0
radius = 1000.0
# Position and velocity will be initialized at runtime for LEO orbit
position = { x = 0.0, y = 0.0, z = 0.0 }
velocity = { x = 0.0, y = 0.0, z = 0.0 }
parent_index = 1 # Earth
color = { r = 1.0, g = 0.0, b = 0.5 }
eccentricity = 0.0
# Semi-major axis will be: Earth radius + 200km
semi_major_axis = 6.571e6 # Placeholder, will be set during initialization
```
--- ---
