#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 <cstdio>

struct TestBody {
    const char* name;
    int body_index;
    int parent_index;
    double expected_period_days;
};

struct StabilityResult {
    const char* name;
    double max_stable_dt;
    double period_days;
};

const double SECONDS_PER_DAY = 86400.0;
const double MIN_DT = 30.0;
const double MAX_DT = 600.0;
const double ENERGY_TOLERANCE = 1.0;
const int NUM_ORBITS = 100;

double calculate_orbital_period(CelestialBody* body, CelestialBody* parent) {
    Vec3 relative_pos = vec3_sub(body->global_position, parent->global_position);
    double r = vec3_magnitude(relative_pos);
    Vec3 relative_vel = vec3_sub(body->global_velocity, parent->global_velocity);
    double v = vec3_magnitude(relative_vel);

    double specific_energy = (v * v) / 2.0 - G * parent->mass / r;
    double semi_major_axis = -G * parent->mass / (2.0 * specific_energy);
    double period_seconds = 2.0 * M_PI * sqrt(pow(semi_major_axis, 3.0) / (G * parent->mass));

    return period_seconds;
}

bool is_dt_stable(SimulationState* sim, const TestBody& test_body, double dt, int num_orbits) {
    SimulationState* test_sim = create_simulation(sim->max_bodies, sim->max_craft, sim->max_maneuvers, dt);

    REQUIRE(load_system_config(test_sim, "tests/informational/test_time_step_stability.toml"));

    int body_index = test_body.body_index;
    int parent_index = test_body.parent_index;

    double initial_energy = calculate_system_total_energy(test_sim);

    Vec3 initial_pos_relative = vec3_sub(
        test_sim->bodies[body_index].global_position,
        test_sim->bodies[parent_index].global_position
    );
    double initial_distance = vec3_magnitude(initial_pos_relative);

    double period = calculate_orbital_period(&test_sim->bodies[body_index], &test_sim->bodies[parent_index]);
    double max_time = period * num_orbits;

    bool completed = true;
    while (test_sim->time < max_time) {
        update_simulation(test_sim);

        if (test_sim->bodies[body_index].parent_index != parent_index) {
            completed = false;
            break;
        }
    }

    double final_energy = calculate_system_total_energy(test_sim);
    double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0;

    Vec3 final_pos_relative = vec3_sub(
        test_sim->bodies[body_index].global_position,
        test_sim->bodies[parent_index].global_position
    );
    double final_distance = vec3_magnitude(final_pos_relative);
    double distance_drift_percent = fabs((final_distance - initial_distance) / initial_distance) * 100.0;

    bool stable = completed && (energy_drift_percent < ENERGY_TOLERANCE) && (distance_drift_percent < 5.0);

    destroy_simulation(test_sim);

    return stable;
}

double find_max_stable_dt(SimulationState* sim, const TestBody& test_body) {
    double low = MIN_DT;
    double high = MAX_DT;
    double max_stable = low;

    printf("Testing %s (period ~%.2f days):\n", test_body.name, test_body.expected_period_days);

    for (int iter = 0; iter < 10; iter++) {
        double mid = (low + high) / 2.0;
        bool stable = is_dt_stable(sim, test_body, mid, NUM_ORBITS);

        if (stable) {
            max_stable = mid;
            low = mid;
            printf("  dt=%.0fs: STABLE\n", mid);
        } else {
            high = mid;
            printf("  dt=%.0fs: UNSTABLE\n", mid);
        }

        if (high - low < 5.0) break;
    }

    printf("  Maximum stable dt: %.0f seconds\n\n", max_stable);
    return max_stable;
}

void print_summary(const StabilityResult* results, int num_results, double min_stable_dt, double default_dt) {
    printf("\n");
    printf("===============================================================================\n");
    printf("                    TIME STEP STABILITY TEST RESULTS\n");
    printf("===============================================================================\n\n");

    printf("STABILITY CRITERIA:\n");
    printf("  - Energy drift < %.1f%% over %d orbits\n", ENERGY_TOLERANCE, NUM_ORBITS);
    printf("  - Distance drift < 5.0%%\n");
    printf("  - No SOI transitions (parent changes)\n\n");

    printf("PER-BODY RESULTS:\n");
    printf("+----------------------+----------------+------------------+----------------+\n");
    printf("| Body                 | Period (days)  | Max Stable dt (s) | Stability Status |\n");
    printf("+----------------------+----------------+------------------+----------------+\n");

    for (int i = 0; i < num_results; i++) {
        const StabilityResult& r = results[i];
        double ratio = default_dt / r.max_stable_dt;
        const char* status = ratio < 0.5 ? "Very Stable" : ratio < 0.8 ? "Stable" : "Limited Margin";
        printf("| %-20s | %14.2f | %16.0f | %-14s |\n", r.name, r.period_days, r.max_stable_dt, status);
    }

    printf("+----------------------+----------------+------------------+----------------+\n\n");

    printf("OVERALL ANALYSIS:\n");
    printf("  Minimum stable time step: %.0f seconds\n", min_stable_dt);
    printf("  Recommended safe dt:     %.0f seconds (0.7x safety margin)\n", min_stable_dt * 0.7);
    printf("  Current default dt:        %.0f seconds\n", default_dt);
    printf("  Current dt stability:     %.0fx\n", default_dt / min_stable_dt);
    printf("\n");

    if (default_dt < min_stable_dt * 0.7) {
        printf("STATUS: Current time step (60s) is VERY STABLE with good margin.\n");
        printf("        Can be increased significantly if needed.\n\n");
    } else if (default_dt < min_stable_dt) {
        printf("STATUS: Current time step (60s) is STABLE with adequate margin.\n");
        printf("        Moderate increases possible.\n\n");
    } else {
        printf("STATUS: Current time step (60s) is near stability limit.\n");
        printf("        Consider reducing for safety.\n\n");
    }

    printf("RECOMMENDATIONS:\n");
    printf("  - For MESSENGER-like close orbits: Keep dt <= %.0f seconds\n", min_stable_dt);
    printf("  - For planetary missions: Current dt=60s is excellent\n");
    printf("  - For Moon-scale orbits: Could use dt=120s+ safely\n");
    printf("\n");
    printf("===============================================================================\n\n");
}

TEST_CASE("Time step stability - Mercury orbiter (MESSENGER-like)", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);

    TestBody mercury_orbiter = {"Mercury_Orbiter", 1, 0, 0.5};
    double max_dt = find_max_stable_dt(sim, mercury_orbiter);

    INFO("Mercury orbiter maximum stable dt: " << max_dt << " seconds");

    REQUIRE(max_dt >= MIN_DT);

    destroy_simulation(sim);
}

TEST_CASE("Time step stability - Io (Jupiter's moon)", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);

    TestBody io = {"Io", 3, 2, 1.77};
    double max_dt = find_max_stable_dt(sim, io);

    INFO("Io maximum stable dt: " << max_dt << " seconds");

    REQUIRE(max_dt >= MIN_DT);

    destroy_simulation(sim);
}

TEST_CASE("Time step stability - Moon (Earth's moon)", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);

    TestBody moon = {"Moon", 5, 4, 27.3};
    double max_dt = find_max_stable_dt(sim, moon);

    INFO("Moon maximum stable dt: " << max_dt << " seconds");

    REQUIRE(max_dt >= MIN_DT);

    destroy_simulation(sim);
}

TEST_CASE("Find minimum stable time step across all bodies", "[timestep][stability]") {
    const double BASE_DT = 60.0;
    SimulationState* sim = create_simulation(10, 0, 0, BASE_DT);

    TestBody bodies[] = {
        {"Mercury_Orbiter", 1, 0, 0.5},
        {"Io", 3, 2, 1.77},
        {"Moon", 5, 4, 27.3}
    };

    StabilityResult results[3];
    double max_dt = MAX_DT;

    printf("\n=== Finding minimum stable dt across all bodies ===\n\n");

    for (int i = 0; i < 3; i++) {
        results[i].name = bodies[i].name;
        results[i].period_days = bodies[i].expected_period_days;
        results[i].max_stable_dt = find_max_stable_dt(sim, bodies[i]);
        if (results[i].max_stable_dt < max_dt) {
            max_dt = results[i].max_stable_dt;
        }
    }

    printf("\n=== RESULTS ===\n");
    printf("Minimum stable time step: %.0f seconds\n", max_dt);
    printf("Recommended safe time step: %.0f seconds (%.0fx safety margin)\n", max_dt * 0.7, 1.0/0.7);

    INFO("Minimum stable dt: " << max_dt << " seconds");

    print_summary(results, 3, max_dt, BASE_DT);

    REQUIRE(max_dt >= MIN_DT);

    destroy_simulation(sim);
}

TEST_CASE("Verify current default dt (60s) stability", "[timestep][stability]") {
    const double DT = 60.0;
    const int NUM_ORBITS = 10;

    SimulationState* sim = create_simulation(10, 0, 0, DT);
    REQUIRE(load_system_config(sim, "tests/informational/test_time_step_stability.toml"));

    struct BodyTest {
        int body_index;
        int parent_index;
        const char* name;
    };

    BodyTest tests[] = {
        {1, 0, "Mercury_Orbiter"},
        {3, 2, "Io"},
        {5, 4, "Moon"}
    };

    for (int t = 0; t < 3; t++) {
        int body_index = tests[t].body_index;
        int parent_index = tests[t].parent_index;
        const char* name = tests[t].name;

        double period = calculate_orbital_period(&sim->bodies[body_index], &sim->bodies[parent_index]);
        double max_time = period * NUM_ORBITS;

        double initial_energy = calculate_system_total_energy(sim);

        INFO("Testing " << name << " with dt=" << DT << "s for " << NUM_ORBITS << " orbits");

        bool completed = true;
        while (sim->time < max_time) {
            update_simulation(sim);

            if (sim->bodies[body_index].parent_index != parent_index) {
                completed = false;
                break;
            }
        }

        double final_energy = calculate_system_total_energy(sim);
        double energy_drift_percent = fabs((final_energy - initial_energy) / initial_energy) * 100.0;

        INFO(name << " completed: " << (completed ? "yes" : "no"));
        INFO(name << " energy drift: " << energy_drift_percent << "%");

        REQUIRE(completed);
        REQUIRE(energy_drift_percent < ENERGY_TOLERANCE);
    }

    destroy_simulation(sim);
}