claudes_game/src/simulation.cpp

#include "simulation.h"
#include <cassert>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>

// Create a new simulation
SimulationState* create_simulation(int max_bodies, double time_step) {
    SimulationState* sim = (SimulationState*)malloc(sizeof(SimulationState));
    sim->bodies = (CelestialBody*)malloc(sizeof(CelestialBody) * max_bodies);
    sim->body_count = 0;
    sim->max_bodies = max_bodies;
    sim->time = 0.0;
    sim->dt = time_step;
    return sim;
}

// Destroy simulation and free memory
void destroy_simulation(SimulationState* sim) {
    if (sim) {
        if (sim->bodies) {
            free(sim->bodies);
        }
        free(sim);
    }
}

// Find which body is gravitationally dominant for the given body
int find_dominant_body(SimulationState* sim, int body_index) {
    if (body_index < 0 || body_index >= sim->body_count) {
        return -1;
    }

    CelestialBody* body = &sim->bodies[body_index];
    int dominant = body->parent_index;

    // Check if we're outside current parent's SOI
    bool outside_current_soi = false;
    if (dominant >= 0 && dominant < sim->body_count) {
        CelestialBody* current_parent = &sim->bodies[dominant];
        double dist_to_current = vec3_distance(body->position, current_parent->position);
        outside_current_soi = dist_to_current > current_parent->soi_radius;
    }

    // Check all other bodies to see if we're within their SOI
    for (int i = 0; i < sim->body_count; i++) {
        if (i == body_index) continue;

        CelestialBody* potential_parent = &sim->bodies[i];
        double distance = vec3_distance(body->position, potential_parent->position);

        // If we're outside current parent's SOI, consider all bodies
        // Otherwise, only consider bodies within their SOI
        bool can_switch = outside_current_soi || distance < potential_parent->soi_radius;

        if (can_switch && i != dominant) {
            if (dominant == -1) {
                dominant = i;
            } else {
                CelestialBody* current_parent = &sim->bodies[dominant];
                double dist_to_current = vec3_distance(body->position, current_parent->position);

                if (outside_current_soi) {
                    // Outside current SOI: switch to closest body
                    if (distance < dist_to_current) {
                        dominant = i;
                    }
                } else {
                    // Inside current SOI: apply hysteresis to prevent oscillations
                    if (distance < dist_to_current * 0.5) {
                        dominant = i;
                    }
                }
            }
        }
    }

    return dominant;
}

// Update sphere of influence radius using Hill sphere approximation
// r_soi = a * (m/M)^(2/5) where a = semi-major axis, m = body mass, M = parent mass
void update_soi(CelestialBody* body, CelestialBody* parent, double semi_major_axis) {
    if (parent == NULL || parent->mass <= 0.0) {
        // Root body (like Sun) has infinite SOI, use a large value
        body->soi_radius = 1e15;  // 1000 AU in meters
        return;
    }

    double mass_ratio = body->mass / parent->mass;
    body->soi_radius = semi_major_axis * pow(mass_ratio, 0.4);  // 2/5 = 0.4
}

void update_simulation(SimulationState* sim) {
    for (int i = 0; i < sim->body_count; i++) {
        CelestialBody* body = &sim->bodies[i];

        if (body->parent_index == -1) {
            continue;
        }

        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;
            }
        }

        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);

    sim->time += sim->dt;
}

static void compute_orbital_velocity_from_vis_viva(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 - 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(r, z_axis);

    if (vec3_magnitude(vel_dir) < 0.01) {
        Vec3 x_axis = {1.0, 0.0, 0.0};
        vel_dir = vec3_cross(r, x_axis);
    }

    vel_dir = vec3_normalize(vel_dir);
    body->velocity = vec3_scale(vel_dir, speed);
    body->velocity = vec3_add(body->velocity, parent->velocity);
}

// FIXME: remove this function since we're already looping in calculate_initial_velocities
static void set_child_bodies_velocity(SimulationState* sim) {
    for (int i = 0; i < sim->body_count; i++) {
        CelestialBody* body = &sim->bodies[i];

        if (body->parent_index == -1) {
            continue;
        }

        if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
            CelestialBody* parent = &sim->bodies[body->parent_index];
            compute_orbital_velocity_from_vis_viva(body, parent);
        }
    }
}

void calculate_initial_velocities(SimulationState* sim) {
    for (int i = 0; i < sim->body_count; i++) {
        if (sim->bodies[i].parent_index == -1) {
            sim->bodies[i].velocity = {0.0, 0.0, 0.0};
        }
    }

    set_child_bodies_velocity(sim);
}

// Calculate SOI radii for all bodies
void calculate_soi_radii(SimulationState* sim) {
    for (int i = 0; i < sim->body_count; i++) {
        CelestialBody* body = &sim->bodies[i];

        if (body->parent_index == -1) {
            body->soi_radius = 1e15;
        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
            CelestialBody* parent = &sim->bodies[body->parent_index];
            update_soi(body, parent, body->semi_major_axis);
        }
    }
}

// FIXME: actually all the above functions should be called inside this loop
//  or the loop should be in config_loader
void initialize_local_coordinates(SimulationState* sim) {
    for (int i = 0; i < sim->body_count; i++) {
        CelestialBody* body = &sim->bodies[i];

        if (body->parent_index == -1) {
            body->local_position = body->position;
            body->local_velocity = body->velocity;
        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
            CelestialBody* parent = &sim->bodies[body->parent_index];
            body->local_position = vec3_sub(body->position, parent->position);
            body->local_velocity = vec3_sub(body->velocity, parent->velocity);
        }
    }
}

void compute_global_coordinates(SimulationState* sim) {
    for (int i = 0; i < sim->body_count; i++) {
        CelestialBody* body = &sim->bodies[i];

        if (body->parent_index == -1) {
            body->position = body->local_position;
            body->velocity = body->local_velocity;
        // FIXME: parent_index > body_count is an error in config file, we
        //  should either verify that in config_loader, or have an assert here
        } else if (body->parent_index >= 0 && body->parent_index < sim->body_count) {
            CelestialBody* parent = &sim->bodies[body->parent_index];
            body->position = vec3_add(body->local_position, parent->position);
            body->velocity = vec3_add(body->local_velocity, parent->velocity);
        }
    }
}

OrbitalElements 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;

    OrbitalElements elem;
    elem.time_days = current_time / SECONDS_PER_DAY;

    Vec3 r_vec = vec3_sub(body->position, primary->position);
    double r = vec3_magnitude(r_vec);
    double v = vec3_magnitude(body->velocity);

    elem.distance_to_sun_au = r / AU;
    elem.velocity_magnitude = v;

    if (optional_ref_body) {
        double dist_ref = vec3_distance(body->position, optional_ref_body->position);
        elem.distance_to_ref_body_au = dist_ref / AU;
    } else {
        elem.distance_to_ref_body_au = -1.0;
    }

    elem.specific_energy = (v * v) / 2.0 - (G * M_sun) / r;

    if (elem.specific_energy < 0) {
        elem.semi_major_axis_au = -(G * M_sun) / (2.0 * elem.specific_energy) / AU;

        double v_squared = v * v;
        double r_dot_v = r_vec.x * body->velocity.x + r_vec.y * body->velocity.y + r_vec.z * body->velocity.z;

        Vec3 e_vec;
        e_vec.x = (v_squared - G * M_sun / r) * r_vec.x - r_dot_v * body->velocity.x;
        e_vec.y = (v_squared - G * M_sun / r) * r_vec.y - r_dot_v * body->velocity.y;
        e_vec.z = (v_squared - G * M_sun / r) * r_vec.z - r_dot_v * body->velocity.z;

        double e_mag = vec3_magnitude(e_vec) / (G * M_sun);
        elem.eccentricity = e_mag;
    } else {
        elem.semi_major_axis_au = 0.0;
        elem.eccentricity = 1.0;
    }

    return elem;
}