#include <catch2/catch_test_macros.hpp>
#include "../src/physics.h"
#include "../src/simulation.h"
#include "../src/config_loader.h"
#include "../src/test_utilities.h"
#include <cmath>
#include <vector>

TEST_CASE("Moon orbital stability around Earth", "[moon][earth]") {
    const double TIME_STEP = 60.0;
    const double EXPECTED_PERIOD_DAYS = 27.3;
    const double SECONDS_PER_DAY = 86400.0;
    const double MAX_SIMULATION_DAYS = 35.0;
    const double MOON_DISTANCE_FROM_EARTH = 384400000.0;

    SimulationState* sim = create_simulation(20, 0, 0, TIME_STEP);

    REQUIRE(load_system_config(sim, "tests/test_moon_orbits.toml"));

    const int EARTH_INDEX = 2;
    const int MOON_INDEX = 8;

    double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY;
    OrbitTracker* tracker = create_orbit_tracker_with_min_time(MOON_INDEX, 5.0);

    int initial_parent = sim->bodies[MOON_INDEX].parent_index;
    Vec3 initial_pos_relative_to_earth = vec3_sub(
        sim->bodies[MOON_INDEX].global_position,
        sim->bodies[EARTH_INDEX].global_position
    );
    double initial_distance = vec3_magnitude(initial_pos_relative_to_earth);

    INFO("Moon initial distance from Earth: " << initial_distance << " m");
    INFO("Expected distance: " << MOON_DISTANCE_FROM_EARTH << " m");

    while (sim->time < max_time && !tracker->orbit_completed) {
        update_simulation(sim);

        int current_parent = sim->bodies[MOON_INDEX].parent_index;

        REQUIRE(current_parent == initial_parent);

        if (current_parent != EARTH_INDEX) {
            INFO("Moon parent changed from " << initial_parent << " to " << current_parent);
            REQUIRE(current_parent == EARTH_INDEX);
        }

        Vec3 current_pos_relative_to_earth = vec3_sub(
            sim->bodies[MOON_INDEX].global_position,
            sim->bodies[EARTH_INDEX].global_position
        );
        double current_distance = vec3_magnitude(current_pos_relative_to_earth);

        double distance_drift = fabs(current_distance - initial_distance);
        double drift_percentage = (distance_drift / initial_distance) * 100.0;

        REQUIRE(drift_percentage < 20.0);

        update_orbit_tracker(tracker, &sim->bodies[MOON_INDEX], &sim->bodies[EARTH_INDEX], sim->time);
    }

    REQUIRE(tracker->orbit_completed);

    double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY;
    double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS);

    INFO("Expected Moon period: " << EXPECTED_PERIOD_DAYS << " days");
    INFO("Measured Moon period: " << measured_period_days << " days");
    INFO("Period error: " << period_error_days << " days");

    REQUIRE(period_error_days < 3.0);

    Vec3 final_pos_relative_to_earth = vec3_sub(
        sim->bodies[MOON_INDEX].global_position,
        sim->bodies[EARTH_INDEX].global_position
    );
    double final_distance = vec3_magnitude(final_pos_relative_to_earth);

    double final_drift_percentage = (fabs(final_distance - initial_distance) / initial_distance) * 100.0;

    INFO("Final distance from Earth: " << final_distance << " m");
    INFO("Initial distance from Earth: " << initial_distance << " m");
    INFO("Final drift: " << final_drift_percentage << "%");

    REQUIRE(final_drift_percentage < 10.0);

    destroy_orbit_tracker(tracker);
    destroy_simulation(sim);
}

TEST_CASE("Galilean moons orbital stability around Jupiter", "[moon][jupiter]") {
    const double TIME_STEP = 60.0;
    const double SECONDS_PER_DAY = 86400.0;
    const double MAX_SIMULATION_DAYS = 20.0;

    const double IO_PERIOD_DAYS = 1.77;
    const double EUROPA_PERIOD_DAYS = 3.55;
    const double GANYMEDE_PERIOD_DAYS = 7.15;
    const double CALLISTO_PERIOD_DAYS = 16.69;

    SimulationState* sim = create_simulation(20, 0, 0, TIME_STEP);

    REQUIRE(load_system_config(sim, "tests/test_moon_orbits.toml"));

    const int JUPITER_INDEX = 4;
    const int IO_INDEX = 9;
    const int EUROPA_INDEX = 10;
    const int GANYMEDE_INDEX = 11;
    const int CALLISTO_INDEX = 12;

    OrbitTracker* io_tracker = create_orbit_tracker_with_min_time(IO_INDEX, 1.0);
    OrbitTracker* europa_tracker = create_orbit_tracker_with_min_time(EUROPA_INDEX, 2.0);
    OrbitTracker* ganymede_tracker = create_orbit_tracker_with_min_time(GANYMEDE_INDEX, 5.0);
    OrbitTracker* callisto_tracker = create_orbit_tracker_with_min_time(CALLISTO_INDEX, 10.0);

    double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY;

    while (sim->time < max_time) {
        update_simulation(sim);

        REQUIRE(sim->bodies[IO_INDEX].parent_index == JUPITER_INDEX);
        REQUIRE(sim->bodies[EUROPA_INDEX].parent_index == JUPITER_INDEX);
        REQUIRE(sim->bodies[GANYMEDE_INDEX].parent_index == JUPITER_INDEX);
        REQUIRE(sim->bodies[CALLISTO_INDEX].parent_index == JUPITER_INDEX);

        update_orbit_tracker(io_tracker, &sim->bodies[IO_INDEX], &sim->bodies[JUPITER_INDEX], sim->time);
        update_orbit_tracker(europa_tracker, &sim->bodies[EUROPA_INDEX], &sim->bodies[JUPITER_INDEX], sim->time);
        update_orbit_tracker(ganymede_tracker, &sim->bodies[GANYMEDE_INDEX], &sim->bodies[JUPITER_INDEX], sim->time);
        update_orbit_tracker(callisto_tracker, &sim->bodies[CALLISTO_INDEX], &sim->bodies[JUPITER_INDEX], sim->time);
    }

    REQUIRE(io_tracker->orbit_completed);
    REQUIRE(europa_tracker->orbit_completed);
    REQUIRE(ganymede_tracker->orbit_completed);
    REQUIRE(callisto_tracker->orbit_completed);

    double io_period_days = io_tracker->time_at_completion / SECONDS_PER_DAY;
    double europa_period_days = europa_tracker->time_at_completion / SECONDS_PER_DAY;
    double ganymede_period_days = ganymede_tracker->time_at_completion / SECONDS_PER_DAY;
    double callisto_period_days = callisto_tracker->time_at_completion / SECONDS_PER_DAY;

    INFO("Io period: " << io_period_days << " days (expected: " << IO_PERIOD_DAYS << ")");
    INFO("Europa period: " << europa_period_days << " days (expected: " << EUROPA_PERIOD_DAYS << ")");
    INFO("Ganymede period: " << ganymede_period_days << " days (expected: " << GANYMEDE_PERIOD_DAYS << ")");
    INFO("Callisto period: " << callisto_period_days << " days (expected: " << CALLISTO_PERIOD_DAYS << ")");

    REQUIRE(fabs(io_period_days - IO_PERIOD_DAYS) < 0.5);
    REQUIRE(fabs(europa_period_days - EUROPA_PERIOD_DAYS) < 1.0);
    REQUIRE(fabs(ganymede_period_days - GANYMEDE_PERIOD_DAYS) < 2.0);
    REQUIRE(fabs(callisto_period_days - CALLISTO_PERIOD_DAYS) < 4.0);

    destroy_orbit_tracker(io_tracker);
    destroy_orbit_tracker(europa_tracker);
    destroy_orbit_tracker(ganymede_tracker);
    destroy_orbit_tracker(callisto_tracker);
    destroy_simulation(sim);
}

TEST_CASE("Titan orbital stability around Saturn", "[moon][saturn]") {
    const double TIME_STEP = 60.0;
    const double EXPECTED_PERIOD_DAYS = 15.95;
    const double SECONDS_PER_DAY = 86400.0;
    const double MAX_SIMULATION_DAYS = 25.0;

    SimulationState* sim = create_simulation(20, 0, 0, TIME_STEP);

    REQUIRE(load_system_config(sim, "tests/test_moon_orbits.toml"));

    const int SATURN_INDEX = 5;
    const int TITAN_INDEX = 13;

    OrbitTracker* tracker = create_orbit_tracker_with_min_time(TITAN_INDEX, 10.0);

    Vec3 initial_pos_relative_to_saturn = vec3_sub(
        sim->bodies[TITAN_INDEX].global_position,
        sim->bodies[SATURN_INDEX].global_position
    );
    double initial_distance = vec3_magnitude(initial_pos_relative_to_saturn);

    double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY;

    while (sim->time < max_time && !tracker->orbit_completed) {
        update_simulation(sim);

        REQUIRE(sim->bodies[TITAN_INDEX].parent_index == SATURN_INDEX);

        Vec3 current_pos_relative_to_saturn = vec3_sub(
            sim->bodies[TITAN_INDEX].global_position,
            sim->bodies[SATURN_INDEX].global_position
        );
        double current_distance = vec3_magnitude(current_pos_relative_to_saturn);

        double drift_percentage = (fabs(current_distance - initial_distance) / initial_distance) * 100.0;
        REQUIRE(drift_percentage < 10.0);

        update_orbit_tracker(tracker, &sim->bodies[TITAN_INDEX], &sim->bodies[SATURN_INDEX], sim->time);
    }

    REQUIRE(tracker->orbit_completed);

    double measured_period_days = tracker->time_at_completion / SECONDS_PER_DAY;
    double period_error_days = fabs(measured_period_days - EXPECTED_PERIOD_DAYS);

    INFO("Expected Titan period: " << EXPECTED_PERIOD_DAYS << " days");
    INFO("Measured Titan period: " << measured_period_days << " days");
    INFO("Period error: " << period_error_days << " days");

    REQUIRE(period_error_days < 3.0);

    destroy_orbit_tracker(tracker);
    destroy_simulation(sim);
}

TEST_CASE("Combined solar system with all moons - parent stability", "[moon][integration]") {
    const double TIME_STEP = 60.0;
    const double SECONDS_PER_DAY = 86400.0;
    const double MAX_SIMULATION_DAYS = 60.0;

    SimulationState* sim = create_simulation(20, 0, 0, TIME_STEP);

    REQUIRE(load_system_config(sim, "tests/test_moon_orbits.toml"));

    struct ParentChange {
        double time_days;
        int body_index;
        int old_parent;
        int new_parent;
    };

    std::vector<ParentChange> parent_changes;

    int initial_parents[14];
    for (int i = 0; i < 14; i++) {
        initial_parents[i] = sim->bodies[i].parent_index;
    }

    double max_time = MAX_SIMULATION_DAYS * SECONDS_PER_DAY;

    while (sim->time < max_time) {
        update_simulation(sim);

        for (int i = 0; i < sim->body_count; i++) {
            if (sim->bodies[i].parent_index != initial_parents[i]) {
                ParentChange change;
                change.time_days = sim->time / SECONDS_PER_DAY;
                change.body_index = i;
                change.old_parent = initial_parents[i];
                change.new_parent = sim->bodies[i].parent_index;
                parent_changes.push_back(change);
                initial_parents[i] = sim->bodies[i].parent_index;
            }
        }
    }

    INFO("Total parent changes detected: " << parent_changes.size());
    for (const auto& change : parent_changes) {
        INFO("Body " << sim->bodies[change.body_index].name
             << " (index " << change.body_index << ") changed parent "
             << "from " << change.old_parent << " to " << change.new_parent
             << " at day " << change.time_days);
    }

    REQUIRE(parent_changes.size() == 0);

    destroy_simulation(sim);
}