From c5946acdc47a496dc3da67260ec0e15fa1d72af7 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sun, 25 Jan 2026 09:04:31 -0500 Subject: [PATCH] WIP: Phase 4-6 - Initialization, validation, and test updates - Added orbital_mechanics module with orbital_elements_to_cartesian() - Updated initialize_orbital_objects() to use orbital mechanics - Added validate_initial_positions() for post-initialization checking - Fixed test files to use global_position/global_velocity - Updated config loader to support spacecraft altitude parameter - Fixed orbital_mechanics.cpp velocity calculation bug (removed duplicate scaling) - Updated Makefile to include orbital_mechanics.o in test build - Renamed simulation.h OrbitalMetrics to OrbitalAnalysis to avoid conflict - Added docs/parabolic_union_implementation.md for parabolic orbit support plan Note: Test configs still need manual fix for orbit table TOML syntax --- Makefile | 1 + docs/parabolic_union_implementation.md | 152 +++++++++++++++++++++++ src/config_loader.cpp | 41 +++--- src/orbital_mechanics.cpp | 9 -- src/simulation.cpp | 45 ++++++- src/simulation.h | 3 + tests/test_hyperbolic_orbit.cpp | 38 +++--- tests/test_invalid_parent_assignment.cpp | 8 +- tests/test_maneuvers.cpp | 12 +- tests/test_moon_orbits.cpp | 20 +-- tests/test_parabolic_orbit.cpp | 22 ++-- tests/test_root_body_transitions.cpp | 12 +- tests/test_soi_transition.cpp | 6 +- 13 files changed, 276 insertions(+), 93 deletions(-) create mode 100644 docs/parabolic_union_implementation.md diff --git a/Makefile b/Makefile index 3b9d6fa..1083503 100644 --- a/Makefile +++ b/Makefile @@ -75,6 +75,7 @@ test-build: $(BUILD_DIR) $(C_OBJECTS) $(CPP_OBJECTS) $(TEST_OBJECTS) $(CXX) $(C_OBJECTS) $(TEST_OBJECTS) \ build/test_utilities.o \ build/physics.o \ + build/orbital_mechanics.o \ build/simulation.o \ build/config_loader.o \ build/maneuver.o \ diff --git a/docs/parabolic_union_implementation.md b/docs/parabolic_union_implementation.md new file mode 100644 index 0000000..8d94027 --- /dev/null +++ b/docs/parabolic_union_implementation.md @@ -0,0 +1,152 @@ +# Parabolic Orbit Union Implementation Plan + +## Overview +Add support for parabolic orbits (e≈1.0) using semi-latus rectum parameter `p` instead of the current hacky `semi_major_axis = 1.0e30` infinity approximation. + +## Problem +Current implementation uses `semi_major_axis = 1.0e30` to approximate infinity for parabolic orbits, causing: +1. Numerical precision issues with extremely large distances (~6.68e18 AU) +2. Velocities approaching zero (1.6e-08 km/s instead of ~42 km/s escape velocity) +3. Test failures due to floating-point equality (final_distance ≈ initial_distance) + +## Solution +Use a union in `OrbitalElements` struct to support both `semi_major_axis` (for elliptical/hyperbolic) and `semi_latus_rectum` (for parabolic). + +## Mathematical Background +For parabolic orbits (e=1.0), the semi-major axis is theoretically infinity. Using semi-latus rectum `p` is mathematically correct: + +Position: `r = p / (1 + cos(ν))` +Velocity: `v = √(2μ / r)` + +Where: +- `p` = semi-latus rectum +- `ν` = true anomaly +- `μ` = GM (gravitational parameter) + +For parabolic orbits: `p = 2q` where `q` is perihelion distance + +## Implementation Steps + +### Phase 1: Update OrbitalElements Struct +**File: `src/orbital_mechanics.h`** + +```cpp +struct OrbitalElements { + union { + double semi_major_axis; // for elliptical (e<1) and hyperbolic (e>1) + double semi_latus_rectum; // for parabolic (e≈1) + }; + double eccentricity; + double true_anomaly; + double inclination; + double longitude_of_ascending_node; + double argument_of_periapsis; +}; +``` + +### Phase 2: Update Config Loader +**File: `src/config_loader.cpp`** + +Add to `parse_toml_body()` and `parse_toml_spacecraft()`: + +1. Parse both `semi_major_axis` and `semi_latus_rectum` from orbit table +2. Initialize union field based on which is specified +3. Validate exactly one is present per eccentricity range + +**Validation Logic:** +```cpp +bool has_semi_major = (semi_major.type == TOML_FP64); +bool has_semi_latus = (semi_latus.type == TOML_FP64); + +if (fabs(elements.eccentricity - 1.0) < 0.005) { + // Parabolic orbit - requires semi_latus_rectum + if (!has_semi_latus) { + printf("Error: Parabolic orbit requires 'semi_latus_rectum'\n"); + return false; + } + if (has_semi_major) { + printf("Error: Parabolic orbit cannot have 'semi_major_axis'\n"); + return false; + } + elements.semi_latus_rectum = semi_latus.u.fp64; +} else { + // Elliptical or hyperbolic - requires semi_major_axis + if (!has_semi_major) { + printf("Error: Elliptical/hyperbolic orbit requires 'semi_major_axis'\n"); + return false; + } + if (has_semi_latus) { + printf("Error: Elliptical/hyperbolic orbit cannot have 'semi_latus_rectum'\n"); + return false; + } + elements.semi_major_axis = semi_major.u.fp64; +} +``` + +### Phase 3: Update orbital_mechanics.cpp +**File: `src/orbital_mechanics.cpp`** + +Update parabolic case (line 21-23): + +```cpp +} else if (fabs(e - 1.0) < 0.005) { + double p = elements.semi_latus_rectum; + r = p / (1.0 + cos(nu)); + v_mag = sqrt(2.0 * mu / r); +} +``` + +Remove the `2.0 * a` approximation that requires `a=1.0e30`. + +### Phase 4: Update Test Configs +**File: `tests/configs/parabolic_comet.toml`** + +Replace `semi_major_axis = 1.0e30` with `semi_latus_rectum = 1.496e11` (p = 1 AU): + +```toml +[[bodies]] +name = "ParabolicComet" +mass = 1.0e14 +radius = 5.0e3 +parent_index = 0 +color = { r = 0.7, g = 0.8, b = 0.9 } +orbit = { + semi_latus_rectum = 1.496e11, + eccentricity = 1.0, + true_anomaly = 0.0 +} +``` + +### Phase 5: Update Documentation +**File: `docs/technical_reference.md`** + +1. Update `OrbitalElements` struct documentation to show union +2. Add note about `semi_latus_rectum` being required for parabolic orbits (e≈1.0) +3. Document `semi_latus_rectum` in config format section + +**File: `docs/unified_orbital_elements_plan.md`** + +Mark union implementation as complete in Phase 7 status. + +## Validation Steps + +1. Build: `make clean && make` +2. Run parabolic test: `./orbit_test '[parabolic]'` +3. Verify velocity is correct: should be ~42,127 m/s escape velocity at 1 AU +4. Verify energy is ~0 (parabolic orbits have total energy = 0) + +## Decisions Made + +### Default Behavior +No backward compatibility for `semi_major_axis` on parabolic orbits - require explicit `semi_latus_rectum` for all parabolic configs. This is cleaner than trying to auto-convert `p = 2*a`. + +### Spacecraft Altitude Parameter +Spacecraft `altitude` parameter is not supported for parabolic orbits in this implementation. If user specifies `altitude` with `eccentricity ≈ 1.0`, the config loader will require `semi_latus_rectum` instead and reject `altitude` or `semi_major_axis`. Added to future todos for Phase 8+. + +### Parabolic Detection Tolerance +Using `|e - 1.0| < 0.005` as threshold for detecting parabolic orbits. This matches tolerance used elsewhere in the codebase. + +## Future Enhancements (TODO) + +- Spacecraft `altitude` parameter for parabolic orbits: parse `altitude` and convert to `semi_latus_rectum = parent_radius + altitude` when eccentricity is parabolic +- Consider adding explicit `perihelion` parameter to config file, then derive `semi_latus_rectum = 2 * perihelion` for parabolic orbits diff --git a/src/config_loader.cpp b/src/config_loader.cpp index 1bb0c5c..2c7b49c 100644 --- a/src/config_loader.cpp +++ b/src/config_loader.cpp @@ -191,41 +191,24 @@ bool load_system_config(SimulationState* sim, const char* filepath) { } } - // Validate parent-child distances - for (int i = 0; i < body_count; i++) { - if (sim->bodies[i].parent_index >= 0) { - CelestialBody* body = &sim->bodies[i]; - CelestialBody* parent = &sim->bodies[body->parent_index]; - - double distance = vec3_distance(body->global_position, parent->global_position); - double min_distance = parent->radius + body->radius; - - if (distance < min_distance) { - printf("Error: Body '%s' (index %d) too close to parent '%s' (index %d)\n", - body->name, i, parent->name, body->parent_index); - printf(" Distance: %.2e m\n", distance); - printf(" Minimum required: %.2e m (parent radius + body radius)\n", min_distance); - toml_free(result); - return false; - } - } - } - // Validate orbital elements for (int i = 0; i < body_count; i++) { CelestialBody* body = &sim->bodies[i]; - // Validate semi_major_axis - if (fabs(body->orbit.semi_major_axis) < 1e-10) { + // Skip validation for root bodies (parent_index=-1) + if (body->parent_index < 0) { + continue; + } + + if (body->orbit.semi_major_axis == 0.0) { printf("Error: Body '%s' has invalid semi_major_axis: %.2e (must not be zero)\n", body->name, body->orbit.semi_major_axis); toml_free(result); return false; } - // Validate eccentricity if (body->orbit.eccentricity < 0.0) { - printf("Error: Body '%s' has invalid eccentricity: %.3f (must be >= 0)\n", + printf("Error: Body '%s' has invalid eccentricity: %.2e (must be >= 0)\n", body->name, body->orbit.eccentricity); toml_free(result); return false; @@ -260,6 +243,11 @@ bool load_system_config(SimulationState* sim, const char* filepath) { toml_free(result); initialize_orbital_objects(sim); + + if (!validate_initial_positions(sim)) { + printf("Error: Initial position validation failed\n"); + return false; + } printf("Loaded %d bodies from %s\n", body_count, filepath); return true; @@ -307,6 +295,11 @@ static bool load_spacecraft_from_toml(SimulationState* sim, toml_result_t result for (int i = 0; i < sim->craft_count; i++) { Spacecraft* craft = &sim->spacecraft[i]; + if (craft->parent_index >= 0 && craft->parent_index < sim->body_count) { + CelestialBody* parent = &sim->bodies[craft->parent_index]; + craft->orbit.semi_major_axis += parent->radius; + } + // Validate semi_major_axis if (fabs(craft->orbit.semi_major_axis) < 1e-10) { printf("Error: Spacecraft '%s' has invalid semi_major_axis: %.2e (must not be zero)\n", diff --git a/src/orbital_mechanics.cpp b/src/orbital_mechanics.cpp index 6ba52bc..f44fac8 100644 --- a/src/orbital_mechanics.cpp +++ b/src/orbital_mechanics.cpp @@ -41,25 +41,16 @@ void orbital_elements_to_cartesian(OrbitalElements elements, double parent_mass, vy_orbital = v_mag; } else if (e < 1.0) { double p = a * (1.0 - e * e); - double h = sqrt(mu * p); vx_orbital = -sqrt(mu / p) * sin_nu; vy_orbital = sqrt(mu / p) * (e + cos_nu); - vx_orbital *= h; - vy_orbital *= h; } else if (fabs(e - 1.0) < 1e-10) { double p = 2.0 * a; - double h = sqrt(mu * p); vx_orbital = -sqrt(mu / p) * sin_nu; vy_orbital = sqrt(mu / p) * (1.0 + cos_nu); - vx_orbital *= h; - vy_orbital *= h; } else { double p = a * (1.0 - e * e); - double h = sqrt(mu * p); vx_orbital = -sqrt(mu / p) * sin_nu; vy_orbital = sqrt(mu / p) * (e + cos_nu); - vx_orbital *= h; - vy_orbital *= h; } Vec3 velocity = {vx_orbital, vy_orbital, 0.0}; diff --git a/src/simulation.cpp b/src/simulation.cpp index 5bf29fa..4a94012 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -177,6 +177,47 @@ double calculate_soi_radius(CelestialBody* body, CelestialBody* parent) { return body->orbit.semi_major_axis * pow(mass_ratio, 0.4); // 2/5 = 0.4 } +// Validate initial positions to ensure bodies aren't inside their parent +bool validate_initial_positions(SimulationState* sim) { + for (int i = 0; i < sim->body_count; i++) { + CelestialBody* body = &sim->bodies[i]; + + if (body->parent_index >= 0 && body->parent_index < sim->body_count) { + CelestialBody* parent = &sim->bodies[body->parent_index]; + double distance = vec3_magnitude(vec3_sub(body->global_position, parent->global_position)); + double min_distance = parent->radius + body->radius; + + if (distance < min_distance) { + printf("Error: Body '%s' (index %d) too close to parent '%s' (index %d)\n", + body->name, i, parent->name, body->parent_index); + printf(" Distance: %.2e m\n", distance); + printf(" Minimum required: %.2e m (parent radius + body radius)\n", min_distance); + return false; + } + } + } + + for (int i = 0; i < sim->craft_count; i++) { + Spacecraft* craft = &sim->spacecraft[i]; + + if (craft->parent_index >= 0 && craft->parent_index < sim->body_count) { + CelestialBody* parent = &sim->bodies[craft->parent_index]; + double distance = vec3_magnitude(vec3_sub(craft->global_position, parent->global_position)); + double min_distance = parent->radius; + + if (distance < min_distance) { + printf("Error: Spacecraft '%s' too close to parent '%s'\n", + craft->name, parent->name); + printf(" Distance: %.2e m\n", distance); + printf(" Minimum required: %.2e m (parent radius)\n", min_distance); + return false; + } + } + } + + return true; +} + // Initialize orbital objects from orbital elements // Converts orbital elements to local position/velocity and computes global coordinates void initialize_orbital_objects(SimulationState* sim) { @@ -211,7 +252,9 @@ void initialize_orbital_objects(SimulationState* sim) { CelestialBody* parent = &sim->bodies[craft->parent_index]; Vec3 local_pos, local_vel; - orbital_elements_to_cartesian(craft->orbit, parent->mass, &local_pos, &local_vel); + OrbitalElements elements = craft->orbit; + + orbital_elements_to_cartesian(elements, parent->mass, &local_pos, &local_vel); craft->local_position = local_pos; craft->local_velocity = local_vel; diff --git a/src/simulation.h b/src/simulation.h index c8714b6..a73d27f 100644 --- a/src/simulation.h +++ b/src/simulation.h @@ -75,6 +75,9 @@ void compute_spacecraft_globals(SimulationState* sim); // Converts orbital elements to local position/velocity and computes global coordinates void initialize_orbital_objects(SimulationState* sim); +// Validate initial positions to ensure bodies aren't inside their parent +bool validate_initial_positions(SimulationState* sim); + // Orbital elements calculation (for output/analysis) struct OrbitalAnalysis { double time_days; diff --git a/tests/test_hyperbolic_orbit.cpp b/tests/test_hyperbolic_orbit.cpp index 83454a0..6a87a5d 100644 --- a/tests/test_hyperbolic_orbit.cpp +++ b/tests/test_hyperbolic_orbit.cpp @@ -19,9 +19,9 @@ TEST_CASE("Hyperbolic orbit - energy and escape trajectory", "[hyperbolic][energ const int COMET_INDEX = 1; const int SUN_INDEX = 0; - Vec3 initial_position = sim->bodies[COMET_INDEX].position; + Vec3 initial_position = sim->bodies[COMET_INDEX].global_position; double initial_distance = vec3_magnitude(initial_position); - double initial_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double initial_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity); double initial_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); double initial_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], @@ -49,8 +49,8 @@ TEST_CASE("Hyperbolic orbit - energy and escape trajectory", "[hyperbolic][energ int step_count = 0; while (sim->time < max_time) { if (step_count % 1000 == 0) { - double current_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); - double current_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double current_distance = vec3_magnitude(sim->bodies[COMET_INDEX].global_position); + double current_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity); double current_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); double current_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], &sim->bodies[SUN_INDEX]); @@ -65,8 +65,8 @@ TEST_CASE("Hyperbolic orbit - energy and escape trajectory", "[hyperbolic][energ step_count++; } - double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); - double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].global_position); + double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity); double final_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); double final_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], @@ -120,8 +120,8 @@ TEST_CASE("Hyperbolic orbit initial conditions", "[hyperbolic][initial]") { CelestialBody* comet = &sim->bodies[COMET_INDEX]; CelestialBody* sun = &sim->bodies[SUN_INDEX]; - double distance = vec3_magnitude(vec3_sub(comet->position, sun->position)); - double velocity = vec3_magnitude(comet->velocity); + double distance = vec3_magnitude(vec3_sub(comet->global_position, sun->global_position)); + double velocity = vec3_magnitude(comet->global_velocity); double escape_velocity = sqrt(2.0 * G * sun->mass / distance); double circular_velocity = sqrt(G * sun->mass / distance); @@ -136,14 +136,14 @@ TEST_CASE("Hyperbolic orbit initial conditions", "[hyperbolic][initial]") { REQUIRE(velocity > escape_velocity); - INFO("Eccentricity: " << comet->eccentricity); + INFO("Eccentricity: " << comet->orbit.eccentricity); - REQUIRE(comet->eccentricity > 1.0); - REQUIRE(fabs(comet->eccentricity - 1.5) < 0.01); + REQUIRE(comet->orbit.eccentricity > 1.0); + REQUIRE(fabs(comet->orbit.eccentricity - 1.5) < 0.01); - INFO("Semi-major axis: " << comet->semi_major_axis / 1.496e11 << " AU"); + INFO("Semi-major axis: " << comet->orbit.semi_major_axis / 1.496e11 << " AU"); - REQUIRE(comet->semi_major_axis < 0.0); + REQUIRE(comet->orbit.semi_major_axis < 0.0); double initial_kinetic = calculate_kinetic_energy(comet); double initial_potential = calculate_potential_energy_pair(comet, sun); @@ -153,9 +153,9 @@ TEST_CASE("Hyperbolic orbit initial conditions", "[hyperbolic][initial]") { REQUIRE(total_energy > 0.0); - Vec3 r_vec = vec3_sub(comet->position, sun->position); + Vec3 r_vec = vec3_sub(comet->global_position, sun->global_position); double r = vec3_magnitude(r_vec); - double a = comet->semi_major_axis; + double a = comet->orbit.semi_major_axis; double expected_v_squared = G * sun->mass * (2.0 / r - 1.0 / a); double expected_velocity = sqrt(expected_v_squared); @@ -182,7 +182,7 @@ TEST_CASE("Hyperbolic orbit asymptotic velocity", "[hyperbolic][asymptotic]") { const int SUN_INDEX = 0; CelestialBody* sun = &sim->bodies[SUN_INDEX]; - double a = sim->bodies[COMET_INDEX].semi_major_axis; + double a = sim->bodies[COMET_INDEX].orbit.semi_major_axis; double expected_v_infinity = sqrt(2.0 * G * sun->mass / fabs(a)); @@ -196,7 +196,7 @@ TEST_CASE("Hyperbolic orbit asymptotic velocity", "[hyperbolic][asymptotic]") { step_count++; if (step_count % 10000 == 0) { - double distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); + double distance = vec3_magnitude(sim->bodies[COMET_INDEX].global_position); if (distance > 20.0 * AU) { INFO("Stopping simulation at distance: " << distance / AU << " AU"); break; @@ -204,8 +204,8 @@ TEST_CASE("Hyperbolic orbit asymptotic velocity", "[hyperbolic][asymptotic]") { } } - double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); - double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].global_position); + double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity); INFO("Final distance: " << final_distance / AU << " AU"); INFO("Final velocity: " << final_velocity / 1000.0 << " km/s"); diff --git a/tests/test_invalid_parent_assignment.cpp b/tests/test_invalid_parent_assignment.cpp index 6584439..79e4a98 100644 --- a/tests/test_invalid_parent_assignment.cpp +++ b/tests/test_invalid_parent_assignment.cpp @@ -95,8 +95,8 @@ TEST_CASE("Invalid parent: detect placeholder config values", const int EARTH_IDX = 1; const int SPACECRAFT_IDX = 3; - Vec3 craft_pos = sim->bodies[SPACECRAFT_IDX].position; - Vec3 earth_pos = sim->bodies[EARTH_IDX].position; + Vec3 craft_pos = sim->bodies[SPACECRAFT_IDX].global_position; + Vec3 earth_pos = sim->bodies[EARTH_IDX].global_position; double distance = vec3_distance(craft_pos, earth_pos); double min_distance = sim->bodies[EARTH_IDX].radius + sim->bodies[SPACECRAFT_IDX].radius; @@ -126,8 +126,8 @@ TEST_CASE("Mutual SOI: similar mass planets within SOI boundary", double planet_a_soi = sim->bodies[PLANET_A_IDX].soi_radius; double planet_b_soi = sim->bodies[PLANET_B_IDX].soi_radius; - double separation = vec3_distance(sim->bodies[PLANET_A_IDX].position, - sim->bodies[PLANET_B_IDX].position); + double separation = vec3_distance(sim->bodies[PLANET_A_IDX].global_position, + sim->bodies[PLANET_B_IDX].global_position); INFO("PlanetA SOI: " << planet_a_soi / 1e9 << " million km"); INFO("PlanetB SOI: " << planet_b_soi / 1e9 << " million km"); diff --git a/tests/test_maneuvers.cpp b/tests/test_maneuvers.cpp index 044fd2d..ab70b3d 100644 --- a/tests/test_maneuvers.cpp +++ b/tests/test_maneuvers.cpp @@ -31,7 +31,7 @@ TEST_CASE("Prograde burn increases orbital energy", "[spacecraft][burn][prograde Spacecraft* craft = &sim->spacecraft[0]; CelestialBody* earth = &sim->bodies[1]; - double initial_distance = vec3_distance(craft->position, earth->position); + double initial_distance = vec3_distance(craft->global_position, earth->global_position); double initial_velocity = vec3_magnitude(craft->local_velocity); apply_impulsive_burn(craft, BURN_PROGRADE, 100.0); @@ -45,7 +45,7 @@ TEST_CASE("Prograde burn increases orbital energy", "[spacecraft][burn][prograde sim_time += TIME_STEP; } - double final_distance = vec3_distance(craft->position, earth->position); + double final_distance = vec3_distance(craft->global_position, earth->global_position); REQUIRE(final_distance > initial_distance); destroy_simulation(sim); @@ -61,7 +61,7 @@ TEST_CASE("Retrograde burn decreases orbital energy", "[spacecraft][burn][retrog Spacecraft* craft = &sim->spacecraft[0]; CelestialBody* earth = &sim->bodies[1]; - double initial_distance = vec3_distance(craft->position, earth->position); + double initial_distance = vec3_distance(craft->global_position, earth->global_position); double initial_velocity = vec3_magnitude(craft->local_velocity); apply_impulsive_burn(craft, BURN_RETROGRADE, 100.0); @@ -75,7 +75,7 @@ TEST_CASE("Retrograde burn decreases orbital energy", "[spacecraft][burn][retrog sim_time += TIME_STEP; } - double final_distance = vec3_distance(craft->position, earth->position); + double final_distance = vec3_distance(craft->global_position, earth->global_position); REQUIRE(final_distance < initial_distance); destroy_simulation(sim); @@ -139,7 +139,7 @@ TEST_CASE("Spacecraft propagation maintains stability", "[spacecraft][propagatio Spacecraft* craft = &sim->spacecraft[0]; CelestialBody* earth = &sim->bodies[1]; - double initial_distance = vec3_distance(craft->position, earth->position); + double initial_distance = vec3_distance(craft->global_position, earth->global_position); double total_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY; double sim_time = 0.0; @@ -148,7 +148,7 @@ TEST_CASE("Spacecraft propagation maintains stability", "[spacecraft][propagatio sim_time += TIME_STEP; } - double final_distance = vec3_distance(craft->position, earth->position); + double final_distance = vec3_distance(craft->global_position, earth->global_position); double distance_drift_percent = fabs((final_distance - initial_distance) / initial_distance) * 100.0; INFO("Initial distance: " << initial_distance << " m"); diff --git a/tests/test_moon_orbits.cpp b/tests/test_moon_orbits.cpp index d25febc..2da4f0d 100644 --- a/tests/test_moon_orbits.cpp +++ b/tests/test_moon_orbits.cpp @@ -25,8 +25,8 @@ TEST_CASE("Moon orbital stability around Earth", "[moon][earth]") { int initial_parent = sim->bodies[MOON_INDEX].parent_index; Vec3 initial_pos_relative_to_earth = vec3_sub( - sim->bodies[MOON_INDEX].position, - sim->bodies[EARTH_INDEX].position + sim->bodies[MOON_INDEX].global_position, + sim->bodies[EARTH_INDEX].global_position ); double initial_distance = vec3_magnitude(initial_pos_relative_to_earth); @@ -46,8 +46,8 @@ TEST_CASE("Moon orbital stability around Earth", "[moon][earth]") { } Vec3 current_pos_relative_to_earth = vec3_sub( - sim->bodies[MOON_INDEX].position, - sim->bodies[EARTH_INDEX].position + sim->bodies[MOON_INDEX].global_position, + sim->bodies[EARTH_INDEX].global_position ); double current_distance = vec3_magnitude(current_pos_relative_to_earth); @@ -71,8 +71,8 @@ TEST_CASE("Moon orbital stability around Earth", "[moon][earth]") { REQUIRE(period_error_days < 3.0); Vec3 final_pos_relative_to_earth = vec3_sub( - sim->bodies[MOON_INDEX].position, - sim->bodies[EARTH_INDEX].position + sim->bodies[MOON_INDEX].global_position, + sim->bodies[EARTH_INDEX].global_position ); double final_distance = vec3_magnitude(final_pos_relative_to_earth); @@ -172,8 +172,8 @@ TEST_CASE("Titan orbital stability around Saturn", "[moon][saturn]") { OrbitTracker* tracker = create_orbit_tracker_with_min_time(TITAN_INDEX, 10.0); Vec3 initial_pos_relative_to_saturn = vec3_sub( - sim->bodies[TITAN_INDEX].position, - sim->bodies[SATURN_INDEX].position + sim->bodies[TITAN_INDEX].global_position, + sim->bodies[SATURN_INDEX].global_position ); double initial_distance = vec3_magnitude(initial_pos_relative_to_saturn); @@ -185,8 +185,8 @@ TEST_CASE("Titan orbital stability around Saturn", "[moon][saturn]") { REQUIRE(sim->bodies[TITAN_INDEX].parent_index == SATURN_INDEX); Vec3 current_pos_relative_to_saturn = vec3_sub( - sim->bodies[TITAN_INDEX].position, - sim->bodies[SATURN_INDEX].position + sim->bodies[TITAN_INDEX].global_position, + sim->bodies[SATURN_INDEX].global_position ); double current_distance = vec3_magnitude(current_pos_relative_to_saturn); diff --git a/tests/test_parabolic_orbit.cpp b/tests/test_parabolic_orbit.cpp index aaf93d8..407ec80 100644 --- a/tests/test_parabolic_orbit.cpp +++ b/tests/test_parabolic_orbit.cpp @@ -19,9 +19,9 @@ TEST_CASE("Parabolic orbit - energy and escape trajectory", "[parabolic][energy] const int COMET_INDEX = 1; const int SUN_INDEX = 0; - Vec3 initial_position = sim->bodies[COMET_INDEX].position; + Vec3 initial_position = sim->bodies[COMET_INDEX].global_position; double initial_distance = vec3_magnitude(initial_position); - double initial_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double initial_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity); double initial_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); double initial_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], @@ -29,7 +29,7 @@ TEST_CASE("Parabolic orbit - energy and escape trajectory", "[parabolic][energy] double initial_total_energy = initial_kinetic + initial_potential; INFO("Initial distance: " << initial_distance / AU << " AU"); - INFO("Initial velocity: " << vec3_magnitude(sim->bodies[COMET_INDEX].velocity) / 1000.0 << " km/s"); + INFO("Initial velocity: " << vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity) / 1000.0 << " km/s"); INFO("Initial kinetic energy: " << initial_kinetic); INFO("Initial potential energy: " << initial_potential); INFO("Initial total energy: " << initial_total_energy); @@ -44,8 +44,8 @@ TEST_CASE("Parabolic orbit - energy and escape trajectory", "[parabolic][energy] int step_count = 0; while (sim->time < max_time) { if (step_count % 1000 == 0) { - double current_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); - double current_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double current_distance = vec3_magnitude(sim->bodies[COMET_INDEX].global_position); + double current_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity); double current_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); double current_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], &sim->bodies[SUN_INDEX]); @@ -60,8 +60,8 @@ TEST_CASE("Parabolic orbit - energy and escape trajectory", "[parabolic][energy] step_count++; } - double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].position); - double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].velocity); + double final_distance = vec3_magnitude(sim->bodies[COMET_INDEX].global_position); + double final_velocity = vec3_magnitude(sim->bodies[COMET_INDEX].global_velocity); double final_kinetic = calculate_kinetic_energy(&sim->bodies[COMET_INDEX]); double final_potential = calculate_potential_energy_pair(&sim->bodies[COMET_INDEX], @@ -113,8 +113,8 @@ TEST_CASE("Parabolic orbit initial conditions", "[parabolic][initial]") { CelestialBody* comet = &sim->bodies[COMET_INDEX]; CelestialBody* sun = &sim->bodies[SUN_INDEX]; - double distance = vec3_magnitude(vec3_sub(comet->position, sun->position)); - double velocity = vec3_magnitude(comet->velocity); + double distance = vec3_magnitude(vec3_sub(comet->global_position, sun->global_position)); + double velocity = vec3_magnitude(comet->global_velocity); double escape_velocity = sqrt(2.0 * G * sun->mass / distance); double circular_velocity = sqrt(G * sun->mass / distance); @@ -129,9 +129,9 @@ TEST_CASE("Parabolic orbit initial conditions", "[parabolic][initial]") { REQUIRE(velocity_error < 0.001); - INFO("Eccentricity: " << comet->eccentricity); + INFO("Eccentricity: " << comet->orbit.eccentricity); - REQUIRE(fabs(comet->eccentricity - 1.0) < 0.0001); + REQUIRE(fabs(comet->orbit.eccentricity - 1.0) < 0.0001); destroy_simulation(sim); } diff --git a/tests/test_root_body_transitions.cpp b/tests/test_root_body_transitions.cpp index 3082edd..b85d0de 100644 --- a/tests/test_root_body_transitions.cpp +++ b/tests/test_root_body_transitions.cpp @@ -60,9 +60,9 @@ TEST_CASE("Root body transition - Earth to Sun", "[root][transition]") { strcpy(event.new_name, "Sun"); } - event.x_position = sim->bodies[PROBE_INDEX].position.x; - event.y_position = sim->bodies[PROBE_INDEX].position.y; - event.z_position = sim->bodies[PROBE_INDEX].position.z; + event.x_position = sim->bodies[PROBE_INDEX].global_position.x; + event.y_position = sim->bodies[PROBE_INDEX].global_position.y; + event.z_position = sim->bodies[PROBE_INDEX].global_position.z; transitions.push_back(event); previous_parent = current_parent; @@ -134,9 +134,9 @@ TEST_CASE("Root body round-trip - Earth -> Sun -> Mars -> Sun", "[root][round-tr strcpy(event.new_name, "Sun"); } - event.x_position = sim->bodies[PROBE_INDEX].position.x; - event.y_position = sim->bodies[PROBE_INDEX].position.y; - event.z_position = sim->bodies[PROBE_INDEX].position.z; + event.x_position = sim->bodies[PROBE_INDEX].global_position.x; + event.y_position = sim->bodies[PROBE_INDEX].global_position.y; + event.z_position = sim->bodies[PROBE_INDEX].global_position.z; transitions.push_back(event); previous_parent = current_parent; diff --git a/tests/test_soi_transition.cpp b/tests/test_soi_transition.cpp index cd7e52e..8e10b02 100644 --- a/tests/test_soi_transition.cpp +++ b/tests/test_soi_transition.cpp @@ -49,9 +49,9 @@ TEST_CASE("SOI transition - Sun to Mars", "[soi][transition]") { change.time_days = sim->time / SECONDS_PER_DAY; change.old_parent = previous_parent; change.new_parent = current_parent; - change.distance_to_mars_au = vec3_distance(sim->bodies[SMALL_BODY_INDEX].position, - sim->bodies[MARS_INDEX].position) / AU; - change.distance_to_sun_au = vec3_magnitude(sim->bodies[SMALL_BODY_INDEX].position) / AU; + change.distance_to_mars_au = vec3_distance(sim->bodies[SMALL_BODY_INDEX].global_position, + sim->bodies[MARS_INDEX].global_position) / AU; + change.distance_to_sun_au = vec3_magnitude(sim->bodies[SMALL_BODY_INDEX].global_position) / AU; parent_changes.push_back(change); previous_parent = current_parent; }