diff --git a/docs/refactor_update_simulation.md b/docs/refactor_update_simulation.md new file mode 100644 index 0000000..232ccda --- /dev/null +++ b/docs/refactor_update_simulation.md @@ -0,0 +1,62 @@ +# Refactoring Plan: Break Up update_simulation + +## Current State +`update_simulation()` function is long with multiple responsibilities: +1. Body RK4 updates and SOI transitions (~60 lines) +2. Body global coordinate computation (~5 lines) +3. Spacecraft RK4 updates (~15 lines) +4. Maneuver execution checking (~20 lines) +5. Spacecraft global coordinate computation (~15 lines) + +## Proposed Helper Functions + +### 1. `update_bodies_physics(SimulationState* sim)` +- Handle body RK4 integration +- Handle SOI transitions +- ~60 lines of current body update logic + +### 2. `compute_body_globals(SimulationState* sim)` +- Compute global coordinates for all bodies +- Reuse existing function or inline it +- ~5 lines + +### 3. `update_spacecraft_physics(SimulationState* sim)` +- RK4 integration for spacecraft +- ~15 lines + +### 4. `execute_pending_maneuvers(SimulationState* sim)` +- Check maneuver triggers +- Execute triggered maneuvers +- ~20 lines + +### 5. `compute_spacecraft_globals(SimulationState* sim)` +- Compute global coordinates for spacecraft +- ~15 lines + +### 6. Simplified `update_simulation(SimulationState* sim)` +```cpp +void update_simulation(SimulationState* sim) { + update_bodies_physics(sim); + compute_body_globals(sim); + update_spacecraft_physics(sim); + execute_pending_maneuvers(sim); + compute_spacecraft_globals(sim); + sim->time += sim->dt; +} +``` + +## Implementation Order +1. Move body update logic to `update_bodies_physics()` +2. Move spacecraft physics logic to `update_spacecraft_physics()` +3. Create `execute_pending_maneuvers()` from existing maneuver loop +4. Move spacecraft globals logic to `compute_spacecraft_globals()` +5. Simplify `update_simulation()` to call helpers +6. Update `simulation.h` with helper declarations +7. Test to ensure behavior is identical + +## Benefits +- Single responsibility for each function +- Easier to test individual components +- Cleaner `update_simulation()` that reads like a sequence +- Easier to debug issues in specific areas +- Functions can be reused or tested independently diff --git a/src/simulation.cpp b/src/simulation.cpp index f5f999e..020f034 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -159,6 +159,82 @@ void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_ax } void update_simulation(SimulationState* sim) { + update_bodies_physics(sim); + compute_global_coordinates(sim); + update_spacecraft_physics(sim); + execute_pending_maneuvers(sim); + compute_spacecraft_globals(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->position, parent->position); + double distance = vec3_magnitude(r); + double e = body->eccentricity; + double a = body->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); + + 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->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 +double calculate_soi_radius(CelestialBody* body, CelestialBody* parent) { + assert(body != nullptr && parent != nullptr); + double mass_ratio = body->mass / parent->mass; + return body->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) { + 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->velocity = calc_orbital_velocity(body, parent); + body->local_position = vec3_sub(body->position, parent->position); + body->local_velocity = vec3_sub(body->velocity, parent->velocity); + body->soi_radius = calculate_soi_radius(body, parent); + } else { // root body + body->velocity = {0.0, 0.0, 0.0}; + body->local_position = body->position; + body->local_velocity = body->velocity; + + // Root body (like Sun) has infinite SOI, use a large value + body->soi_radius = 1e15; // 1000 AU in meters + } + } +} + +// Simulation update helper functions +void update_bodies_physics(SimulationState* sim) { for (int i = 0; i < sim->body_count; i++) { CelestialBody* body = &sim->bodies[i]; @@ -169,27 +245,22 @@ void update_simulation(SimulationState* sim) { int new_parent = find_dominant_body(sim, i); if (new_parent != body->parent_index) { - // Convert current local coordinates to global coordinates using old parent if (body->parent_index >= 0 && body->parent_index < sim->body_count) { CelestialBody* old_parent = &sim->bodies[body->parent_index]; body->position = vec3_add(body->local_position, old_parent->position); body->velocity = vec3_add(body->local_velocity, old_parent->velocity); } else { - // old_parent is root (Sun): local = global body->position = body->local_position; body->velocity = body->local_velocity; } - // Update parent index body->parent_index = new_parent; - // Convert global coordinates to local coordinates using new parent if (body->parent_index >= 0 && body->parent_index < sim->body_count) { CelestialBody* new_parent_body = &sim->bodies[body->parent_index]; body->local_position = vec3_sub(body->position, new_parent_body->position); body->local_velocity = vec3_sub(body->velocity, new_parent_body->velocity); } else { - // new_parent is root (Sun): global = local body->local_position = body->position; body->local_velocity = body->velocity; } @@ -197,14 +268,14 @@ void update_simulation(SimulationState* sim) { if (body->parent_index >= 0 && body->parent_index < sim->body_count) { CelestialBody* parent = &sim->bodies[body->parent_index]; + rk4_step(&body->local_position, &body->local_velocity, sim->dt, body->mass, parent->mass); } } +} - compute_global_coordinates(sim); - - // Update spacecraft +void update_spacecraft_physics(SimulationState* sim) { for (int i = 0; i < sim->craft_count; i++) { Spacecraft* craft = &sim->spacecraft[i]; @@ -217,30 +288,32 @@ void update_simulation(SimulationState* sim) { rk4_step(&craft->local_position, &craft->local_velocity, sim->dt, craft->mass, parent->mass); } - - // Execute pending maneuvers (before computing globals) +} + +void execute_pending_maneuvers(SimulationState* sim) { for (int i = 0; i < sim->maneuver_count; i++) { Maneuver* maneuver = &sim->maneuvers[i]; - + if (maneuver->executed) { continue; } - + if (maneuver->craft_index < 0 || maneuver->craft_index >= sim->craft_count) { continue; } - + Spacecraft* craft = &sim->spacecraft[maneuver->craft_index]; - + if (check_maneuver_trigger(maneuver, craft, sim)) { execute_maneuver(maneuver, craft, sim->time); } } - - // Compute spacecraft global coordinates (after maneuvers) +} + +void compute_spacecraft_globals(SimulationState* sim) { 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->position = vec3_add(parent->position, craft->local_position); @@ -250,75 +323,6 @@ void update_simulation(SimulationState* sim) { craft->velocity = craft->local_velocity; } } - - 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->position, parent->position); - double distance = vec3_magnitude(r); - double e = body->eccentricity; - double a = body->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); - - 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->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 -double calculate_soi_radius(CelestialBody* body, CelestialBody* parent) { - assert(body != nullptr && parent != nullptr); - double mass_ratio = body->mass / parent->mass; - return body->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) { - 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->velocity = calc_orbital_velocity(body, parent); - body->local_position = vec3_sub(body->position, parent->position); - body->local_velocity = vec3_sub(body->velocity, parent->velocity); - body->soi_radius = calculate_soi_radius(body, parent); - } else { // root body - body->velocity = {0.0, 0.0, 0.0}; - body->local_position = body->position; - body->local_velocity = body->velocity; - - // Root body (like Sun) has infinite SOI, use a large value - body->soi_radius = 1e15; // 1000 AU in meters - } - } } void compute_global_coordinates(SimulationState* sim) { diff --git a/src/simulation.h b/src/simulation.h index c5e9df6..59bf24d 100644 --- a/src/simulation.h +++ b/src/simulation.h @@ -57,6 +57,12 @@ double calculate_soi_radius(CelestialBody* body, CelestialBody* parent); void compute_global_coordinates(SimulationState* sim); +// Simulation update helper functions +void update_bodies_physics(SimulationState* sim); +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);