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

struct ParentChange {
    double time_seconds;
    double time_days;
    int old_parent;
    int new_parent;
    double distance_to_mars_au;
    double distance_to_sun_au;
};

TEST_CASE("SOI transition - Sun to Mars", "[soi][transition]") {
    const double TIME_STEP = 60.0;
    const double DAYS_TO_SIMULATE = 2000.0;
    const double SECONDS_PER_DAY = 86400.0;
    const double AU = 1.496e11;

    INFO("Note: SOI transitions use simple model (no hysteresis)");

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

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

    const int SMALL_BODY_INDEX = 2;
    const int MARS_INDEX = 1;
    const int SUN_INDEX = 0;

    std::vector<ParentChange> parent_changes;
    int previous_parent = sim->bodies[SMALL_BODY_INDEX].parent_index;

    INFO("Initial parent index: " << previous_parent);
    INFO("Expected: Sun (0)");

    REQUIRE(previous_parent == SUN_INDEX);

    double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY;
    while (sim->time < max_time) {
        update_simulation(sim);

        int current_parent = sim->bodies[SMALL_BODY_INDEX].parent_index;
        if (current_parent != previous_parent) {
            ParentChange change;
            change.time_seconds = sim->time;
            change.time_days = sim->time / SECONDS_PER_DAY;
            change.old_parent = previous_parent;
            change.new_parent = current_parent;
            change.distance_to_mars_au = vec3_distance(sim->bodies[SMALL_BODY_INDEX].global_position,
                                                       sim->bodies[MARS_INDEX].global_position) / AU;
            change.distance_to_sun_au = vec3_magnitude(sim->bodies[SMALL_BODY_INDEX].global_position) / AU;
            parent_changes.push_back(change);
            previous_parent = current_parent;
        }
    }

    INFO("Total parent changes: " << parent_changes.size());

    REQUIRE(parent_changes.size() > 0);

    bool found_mars_transition = false;

    for (const auto& change : parent_changes) {
        if (change.new_parent == MARS_INDEX || change.old_parent == MARS_INDEX) {
            found_mars_transition = true;
            if (change.new_parent == MARS_INDEX) {
                INFO("Entered Mars SOI at day " << change.time_days
                     << ", distance to Mars: " << change.distance_to_mars_au << " AU"
                     << ", distance to Sun: " << change.distance_to_sun_au << " AU");
            } else if (change.old_parent == MARS_INDEX) {
                INFO("Exited Mars SOI at day " << change.time_days
                     << ", distance to Mars: " << change.distance_to_mars_au << " AU"
                     << ", distance to Sun: " << change.distance_to_sun_au << " AU");
            }
        }
    }

    REQUIRE(found_mars_transition);

    bool found_mars_exit = false;
    for (const auto& change : parent_changes) {
        if (change.old_parent == MARS_INDEX && change.new_parent == SUN_INDEX) {
            found_mars_exit = true;
            INFO("Mars -> Sun exit at day " << change.time_days
                 << ", distance to Mars: " << change.distance_to_mars_au << " AU"
                 << ", distance to Sun: " << change.distance_to_sun_au << " AU");
            break;
        }
    }

    REQUIRE(found_mars_exit);

    destroy_simulation(sim);
}

TEST_CASE("SOI transition - verify SOI radii", "[soi][radii]") {
    const double TIME_STEP = 60.0;
    const double AU = 1.496e11;

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

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

    const int MARS_INDEX = 1;

    double mars_soi_meters = sim->bodies[MARS_INDEX].soi_radius;
    double mars_soi_au = mars_soi_meters / AU;

    INFO("Mars SOI radius: " << mars_soi_meters / 1e9 << " million km");
    INFO("Mars SOI radius: " << mars_soi_au << " AU");

    double expected_mars_soi = 2.244e11 * pow(6.39e23 / 1.989e30, 0.4);
    double expected_mars_soi_au = expected_mars_soi / AU;

    INFO("Expected Mars SOI: " << expected_mars_soi / 1e9 << " million km");
    INFO("Expected Mars SOI: " << expected_mars_soi_au << " AU");

    double soi_error = fabs(mars_soi_au - expected_mars_soi_au) / expected_mars_soi_au;
    INFO("SOI calculation error: " << soi_error * 100.0 << "%");

    REQUIRE(soi_error < 0.01);

    REQUIRE(mars_soi_au > 0.003);
    REQUIRE(mars_soi_au < 0.005);

    destroy_simulation(sim);
}