#include <catch2/catch_test_macros.hpp>
#include "../src/physics.h"
#include "../src/bodies.h"
#include "../src/config_loader.h"
#include "../src/test_utilities.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("Comet orbital elements and SOI transitions during Mars encounter", "[comet][orbital-elements][soi]") {
    const double TIME_STEP = 60.0;
    const double SECONDS_PER_DAY = 86400.0;
    const double DAYS_TO_SIMULATE = 5000.0;
    const double AU = 1.496e11;

    const double EXPECTED_SMA = 2.5;
    const double EXPECTED_ECC = 0.7;

    SimulationState* sim = create_simulation(10, TIME_STEP);

    REQUIRE(load_system_config(sim, "configs/test_simple.txt"));

    const int COMET_INDEX = 3;
    const int MARS_INDEX = 2;
    const int SUN_INDEX = 0;

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

    OrbitalElements initial = calculate_orbital_elements(
        &sim->bodies[COMET_INDEX],
        &sim->bodies[SUN_INDEX],
        &sim->bodies[MARS_INDEX],
        sim->time
    );
    snapshots.push_back(initial);

    double checkpoint_days[] = {0, 365, 722, 1444, 1533, 1600, 1700, 2000, 5000};
    size_t next_checkpoint = 1;

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

        int current_parent = sim->bodies[COMET_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_sun_au = vec3_magnitude(sim->bodies[COMET_INDEX].position) / AU;
            change.distance_to_mars_au = vec3_distance(sim->bodies[COMET_INDEX].position,
                                                       sim->bodies[MARS_INDEX].position) / AU;
            parent_changes.push_back(change);
            previous_parent = current_parent;
        }

        double current_days = sim->time / SECONDS_PER_DAY;
        if (next_checkpoint < sizeof(checkpoint_days)/sizeof(checkpoint_days[0]) &&
            current_days >= checkpoint_days[next_checkpoint]) {

            printf("\n day: %f, sim->time: %f, comet parent index: %d",
					current_days,
					sim->time,
					sim->bodies[COMET_INDEX].parent_index);

            OrbitalElements elem = calculate_orbital_elements(
                &sim->bodies[COMET_INDEX],
                &sim->bodies[SUN_INDEX],
                &sim->bodies[MARS_INDEX],
                sim->time
            );
            snapshots.push_back(elem);
            next_checkpoint++;
        }
    }

#if 1
    printf("\n=== Comet Orbital Elements Over Time ===\n");
    printf("Expected: a = %.1f AU, e = %.1f\n\n", EXPECTED_SMA, EXPECTED_ECC);

    for (const auto& elem : snapshots) {
        printf("Day %.0f:\n", elem.time_days);
        printf("  Semi-major axis: %.6f AU (expected: %.1f AU)\n", elem.semi_major_axis_au, EXPECTED_SMA);
        printf("  Eccentricity: %.6f (expected: %.1f)\n", elem.eccentricity, EXPECTED_ECC);
        printf("  Distance to Sun: %.4f AU\n", elem.distance_to_sun_au);
        printf("  Distance to Mars: %.6f AU\n", elem.distance_to_ref_body_au);
        printf("  Velocity: %.2f km/s\n", elem.velocity_magnitude / 1000.0);
        printf("  Specific energy: %.3e J/kg\n", elem.specific_energy);

        double sma_error = fabs(elem.semi_major_axis_au - EXPECTED_SMA);
        double ecc_error = fabs(elem.eccentricity - EXPECTED_ECC);

        if (elem.time_days >= 1530 && elem.time_days <= 1540) {
            printf("  *** CLOSEST APPROACH TO MARS ***\n");
        }

        printf("  Error: Δa = %.6f AU, Δe = %.6f\n\n", sma_error, ecc_error);
    }

    printf("\n=== Parent Index Changes ===\n");

    if (parent_changes.empty()) {
        printf("No parent changes detected\n\n");
    } else {
        for (const auto& change : parent_changes) {
            printf("Time: %.0f days (%.0f seconds)\n", change.time_days, change.time_seconds);
            printf("  Parent: %d -> %d", change.old_parent, change.new_parent);
            if (change.new_parent == MARS_INDEX) {
                printf(" (Sun -> Mars)");
            } else if (change.old_parent == MARS_INDEX) {
                printf(" (Mars -> Sun)");
            }
            printf("\n");
            printf("  Distance to Mars: %.6f AU\n", change.distance_to_mars_au);
            printf("  Distance to Sun: %.6f AU\n\n", change.distance_to_sun_au);
        }
    }
#endif
    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;
            REQUIRE(fabs(change.time_days - 1550.0) < 50.0);
            INFO("Mars encounter at day " << change.time_days << ", distance " << change.distance_to_mars_au << " AU");
        }
    }
    REQUIRE(found_mars_transition);

    destroy_simulation(sim);
}