#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;
    sim->config_name[0] = '\0';
    return sim;
}

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

// Add a body to the simulation at runtime
int add_body_to_simulation(SimulationState* sim, CelestialBody* body) {
    if (sim->body_count >= sim->max_bodies) {
        printf("Error: Cannot add body - simulation full (%d/%d)\n",
               sim->body_count, sim->max_bodies);
        return -1;
    }

    int new_idx = sim->body_count;
    sim->bodies[new_idx] = *body;
    sim->body_count++;

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

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

    sim->bodies[new_idx].position = body->position;
    sim->bodies[new_idx].velocity = body->velocity;

    return new_idx;
}

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 parent_idx = body->parent_index;

    // If parent is not root (not Sun): only check if still within parent's SOI
    if (parent_idx != 0) {
        if (parent_idx < 0 || parent_idx >= sim->body_count) {
            return -1;
        }

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

        // Stay with parent if within SOI, otherwise go to Sun
        if (distance < parent->soi_radius) {
            return parent_idx;
        } else {
            return 0;
        }
    }

    // Parent is root (Sun): check all bodies for SOI containment
    int new_parent = 0;
    double min_distance = INFINITY;

    for (int i = 0; i < sim->body_count; i++) {
        if (i == body_index) continue;

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

        // If within SOI and closer than current, switch to this body
        if (distance < potential->soi_radius && distance < min_distance) {
            min_distance = distance;
            new_parent = i;
        }
    }

    return new_parent;
}
// 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;
}

// 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) {
    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;
        } 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;
}