From dade860e89b722d38222cfe4e857aa87e2a6cea7 Mon Sep 17 00:00:00 2001 From: cinnaboot Date: Sat, 24 Jan 2026 14:55:26 -0500 Subject: [PATCH] Phase 4: Update initialization to use orbital mechanics module --- docs/unified_orbital_elements_plan.md | 57 +++++++++---------- src/config_loader.cpp | 2 +- src/maneuver.cpp | 2 +- src/simulation.cpp | 82 +++++++++++---------------- src/simulation.h | 9 +-- src/test_utilities.cpp | 14 ++--- 6 files changed, 76 insertions(+), 90 deletions(-) diff --git a/docs/unified_orbital_elements_plan.md b/docs/unified_orbital_elements_plan.md index abea434..1d964f0 100644 --- a/docs/unified_orbital_elements_plan.md +++ b/docs/unified_orbital_elements_plan.md @@ -245,32 +245,32 @@ struct Spacecraft { - For elliptical orbits (e < 1): `semi_major_axis > 0` - For hyperbolic orbits (e > 1): `semi_major_axis` negative or handle separately -### Phase 4: Update Initialization - -10. **Replace `initialize_bodies()` with `initialize_orbital_objects()` (simulation.cpp)** - - For each body: - - If `parent_index >= 0`: - - Call `orbital_elements_to_cartesian(body->orbit, parent->mass, &local_pos, &local_vel)` - - Set `body->local_position = local_pos`, `body->local_velocity = local_vel` - - Compute global: `body->global_position = vec3_add(parent->global_position, local_pos)` - - Compute global: `body->global_velocity = vec3_add(parent->global_velocity, local_vel)` - - Calculate SOI: `body->soi_radius = calculate_soi_radius(body, parent)` - - Else (root body): - - Set all positions/velocities to zero - - Set `soi_radius = 1e15` - - - For each spacecraft: - - Call `orbital_elements_to_cartesian(craft->orbit, parent->mass, &local_pos, &local_vel)` - - Set `craft->local_position = local_pos`, `craft->local_velocity = local_vel` - - Compute global: `craft->global_position = vec3_add(parent->global_position, local_pos)` - - Compute global: `craft->global_velocity = vec3_add(parent->global_velocity, local_vel)` - -11. **Remove old `calc_orbital_velocity()` function** (simulation.cpp) - - No longer needed, replaced by orbital mechanics module - -12. **Update all calls to `initialize_bodies()` → `initialize_orbital_objects()`** - - `config_loader.cpp`: Update call after config load - - Tests: Update if they call initialization directly +### Phase 4: Update Initialization ✅ COMPLETE + +10. ✅ **Replace `initialize_bodies()` with `initialize_orbital_objects()` (simulation.cpp)** + - For each body: + - If `parent_index >= 0`: + - Call `orbital_elements_to_cartesian(body->orbit, parent->mass, &local_pos, &local_vel)` + - Set `body->local_position = local_pos`, `body->local_velocity = local_vel` + - Compute global: `body->global_position = vec3_add(parent->global_position, local_pos)` + - Compute global: `body->global_velocity = vec3_add(parent->global_velocity, local_vel)` + - Calculate SOI: `body->soi_radius = calculate_soi_radius(body, parent)` + - Else (root body): + - Set all positions/velocities to zero + - Set `soi_radius = 1e15` + + - For each spacecraft: + - Call `orbital_elements_to_cartesian(craft->orbit, parent->mass, &local_pos, &local_vel)` + - Set `craft->local_position = local_pos`, `craft->local_velocity = local_vel` + - Compute global: `craft->global_position = vec3_add(parent->global_position, local_pos)` + - Compute global: `craft->global_velocity = vec3_add(parent->global_velocity, local_vel)` + +11. ✅ **Remove old `calc_orbital_velocity()` function** (simulation.cpp) + - No longer needed, replaced by orbital mechanics module + +12. ✅ **Update all calls to `initialize_bodies()` → `initialize_orbital_objects()`** + - `config_loader.cpp`: Update call after config load + - Tests: Update if they call initialization directly ### Phase 5: Update Test Configs @@ -463,12 +463,11 @@ This is a **major breaking change** that affects all configs and code: - Help identify old vs new config formats (if we ever need backward compat) ## Success Criteria - 1. [ ] All test configs updated to new format 2. [ ] Config loader successfully parses `orbit` tables 3. [ ] `orbital_elements_to_cartesian()` correctly converts all orbit types -4. [ ] `initialize_orbital_objects()` sets up bodies and spacecraft correctly -5. [ ] All references to `position`/`velocity` renamed to `global_*` +4. [x] `initialize_orbital_objects()` sets up bodies and spacecraft correctly +5. [x] All references to `position`/`velocity` renamed to `global_*` 6. [ ] Renderer displays orbits correctly 7. [ ] Full test suite passes (`make test`) 8. [ ] Solar system simulation runs correctly diff --git a/src/config_loader.cpp b/src/config_loader.cpp index 3ca1508..1bb0c5c 100644 --- a/src/config_loader.cpp +++ b/src/config_loader.cpp @@ -259,7 +259,7 @@ bool load_system_config(SimulationState* sim, const char* filepath) { toml_free(result); - initialize_bodies(sim); + initialize_orbital_objects(sim); printf("Loaded %d bodies from %s\n", body_count, filepath); return true; diff --git a/src/maneuver.cpp b/src/maneuver.cpp index 5fa1eae..9484d26 100644 --- a/src/maneuver.cpp +++ b/src/maneuver.cpp @@ -63,7 +63,7 @@ void apply_impulsive_burn(Spacecraft* craft, BurnDirection direction, double del void apply_custom_burn(Spacecraft* craft, Vec3 delta_v_local) { craft->local_velocity = vec3_add(craft->local_velocity, delta_v_local); - craft->velocity = vec3_add(craft->velocity, delta_v_local); + craft->global_velocity = vec3_add(craft->global_velocity, delta_v_local); } double calculate_orbital_velocity(Spacecraft* craft, double parent_mass) { diff --git a/src/simulation.cpp b/src/simulation.cpp index aa0b621..5bf29fa 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,6 +1,7 @@ #include "simulation.h" #include "spacecraft.h" #include "maneuver.h" +#include "orbital_mechanics.h" #include #include #include @@ -167,41 +168,6 @@ void update_simulation(SimulationState* sim) { sim->time += sim->dt; } -// Calculate orbital velocity using vis-viva equation -// Returns velocity vector for body relative to parent -static Vec3 calc_orbital_velocity(CelestialBody* body, CelestialBody* parent) { - Vec3 r = vec3_sub(body->global_position, parent->global_position); - double distance = vec3_magnitude(r); - double e = body->orbit.eccentricity; - double a = body->orbit.semi_major_axis; - double v_squared; - - if (fabs(e) < 0.0001) { - v_squared = G * parent->mass / a; - } else if (fabs(e -1.0) < 0.0001) { - v_squared = 2.0 * G * parent->mass / distance; - } else { - v_squared = G * parent->mass * (2.0 / distance - 1.0 / a); - } - - assert(v_squared >= 0); - double speed = (double) sqrt(v_squared); - - // FIXME: this whole section is just wrong and will break with inclined orbits - // it's also a failure of our testing that we're not catching it - Vec3 z_axis = {0.0, 0.0, 1.0}; - Vec3 vel_dir = vec3_cross(z_axis, r); - - if (vec3_magnitude(vel_dir) < 0.01) { - Vec3 x_axis = {1.0, 0.0, 0.0}; - vel_dir = vec3_cross(z_axis, r); - } - - vel_dir = vec3_normalize(vel_dir); - Vec3 velocity = vec3_scale(vel_dir, speed); - return vec3_add(velocity, parent->global_velocity); -} - // Calculate SOI radius for a single body // r_soi = a * (m/M)^(2/5) where a = semi-major axis, m = body mass, M = parent mass // Returns SOI radius in meters @@ -211,26 +177,46 @@ double calculate_soi_radius(CelestialBody* body, CelestialBody* parent) { return body->orbit.semi_major_axis * pow(mass_ratio, 0.4); // 2/5 = 0.4 } -// Combined initialization - sets velocities, SOI radii, and local coordinates in single loop -void initialize_bodies(SimulationState* sim) { +// Initialize orbital objects from orbital elements +// Converts orbital elements to local position/velocity and computes global coordinates +void initialize_orbital_objects(SimulationState* sim) { for (int i = 0; i < sim->body_count; i++) { CelestialBody* body = &sim->bodies[i]; CelestialBody* parent = NULL; - // Set parent pointer if not root body if (body->parent_index >= 0 && body->parent_index < sim->body_count) { parent = &sim->bodies[body->parent_index]; - body->global_velocity = calc_orbital_velocity(body, parent); - body->local_position = vec3_sub(body->global_position, parent->global_position); - body->local_velocity = vec3_sub(body->global_velocity, parent->global_velocity); + + Vec3 local_pos, local_vel; + orbital_elements_to_cartesian(body->orbit, parent->mass, &local_pos, &local_vel); + + body->local_position = local_pos; + body->local_velocity = local_vel; + body->global_position = vec3_add(parent->global_position, local_pos); + body->global_velocity = vec3_add(parent->global_velocity, local_vel); body->soi_radius = calculate_soi_radius(body, parent); - } else { // root body + } else { + body->local_position = {0.0, 0.0, 0.0}; + body->local_velocity = {0.0, 0.0, 0.0}; + body->global_position = {0.0, 0.0, 0.0}; body->global_velocity = {0.0, 0.0, 0.0}; - body->local_position = body->global_position; - body->local_velocity = body->global_velocity; + body->soi_radius = 1e15; + } + } + + 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]; + + Vec3 local_pos, local_vel; + orbital_elements_to_cartesian(craft->orbit, parent->mass, &local_pos, &local_vel); - // Root body (like Sun) has infinite SOI, use a large value - body->soi_radius = 1e15; // 1000 AU in meters + craft->local_position = local_pos; + craft->local_velocity = local_vel; + craft->global_position = vec3_add(parent->global_position, local_pos); + craft->global_velocity = vec3_add(parent->global_velocity, local_vel); } } } @@ -342,13 +328,13 @@ void compute_global_coordinates(SimulationState* sim) { } } -OrbitalMetrics calculate_orbital_elements(CelestialBody* body, CelestialBody* primary, +OrbitalAnalysis calculate_orbital_elements(CelestialBody* body, CelestialBody* primary, CelestialBody* optional_ref_body, double current_time) { const double AU = 1.496e11; const double SECONDS_PER_DAY = 86400.0; const double M_sun = primary->mass; - OrbitalMetrics elem; + OrbitalAnalysis elem; elem.time_days = current_time / SECONDS_PER_DAY; Vec3 r_vec = vec3_sub(body->global_position, primary->global_position); diff --git a/src/simulation.h b/src/simulation.h index d82223d..c8714b6 100644 --- a/src/simulation.h +++ b/src/simulation.h @@ -71,11 +71,12 @@ void update_spacecraft_physics(SimulationState* sim); void execute_pending_maneuvers(SimulationState* sim); void compute_spacecraft_globals(SimulationState* sim); -// Combined initialization - sets velocities, SOI radii, and local coordinates in single loop -void initialize_bodies(SimulationState* sim); +// Initialize orbital objects from orbital elements +// Converts orbital elements to local position/velocity and computes global coordinates +void initialize_orbital_objects(SimulationState* sim); // Orbital elements calculation (for output/analysis) -struct OrbitalMetrics { +struct OrbitalAnalysis { double time_days; double semi_major_axis_au; double eccentricity; @@ -85,7 +86,7 @@ struct OrbitalMetrics { double velocity_magnitude; }; -OrbitalMetrics calculate_orbital_elements(CelestialBody* body, CelestialBody* primary, +OrbitalAnalysis calculate_orbital_elements(CelestialBody* body, CelestialBody* primary, CelestialBody* optional_ref_body, double current_time); #endif diff --git a/src/test_utilities.cpp b/src/test_utilities.cpp index e53738e..430df48 100644 --- a/src/test_utilities.cpp +++ b/src/test_utilities.cpp @@ -3,14 +3,14 @@ #include double calculate_kinetic_energy(CelestialBody* body) { - double v_squared = body->velocity.x * body->velocity.x + - body->velocity.y * body->velocity.y + - body->velocity.z * body->velocity.z; + double v_squared = body->global_velocity.x * body->global_velocity.x + + body->global_velocity.y * body->global_velocity.y + + body->global_velocity.z * body->global_velocity.z; return 0.5 * body->mass * v_squared; } double calculate_potential_energy_pair(CelestialBody* body1, CelestialBody* body2) { - double distance = vec3_distance(body1->position, body2->position); + double distance = vec3_distance(body1->global_position, body2->global_position); if (distance < 1.0) distance = 1.0; return -G * body1->mass * body2->mass / distance; } @@ -33,9 +33,9 @@ double calculate_system_total_energy(SimulationState* sim) { OrbitalMetrics calculate_orbital_metrics(CelestialBody* body, CelestialBody* parent) { OrbitalMetrics metrics; - Vec3 relative_pos = vec3_sub(body->position, parent->position); + Vec3 relative_pos = vec3_sub(body->global_position, parent->global_position); metrics.orbital_radius = vec3_magnitude(relative_pos); - metrics.velocity_magnitude = vec3_magnitude(body->velocity); + metrics.velocity_magnitude = vec3_magnitude(body->global_velocity); metrics.angular_position = atan2(relative_pos.y, relative_pos.x); metrics.kinetic_energy = calculate_kinetic_energy(body); @@ -80,7 +80,7 @@ void reset_orbit_tracker(OrbitTracker* tracker) { void update_orbit_tracker(OrbitTracker* tracker, CelestialBody* body, CelestialBody* parent, double current_time) { if (tracker->orbit_completed) return; - Vec3 relative_pos = vec3_sub(body->position, parent->position); + Vec3 relative_pos = vec3_sub(body->global_position, parent->global_position); double current_angle = atan2(relative_pos.y, relative_pos.x); if (tracker->quadrant_transitions == 0) {