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

using Catch::Matchers::WithinAbs;

SCENARIO("Parabolic orbit - escape trajectory and initial conditions",
         "[parabolic][energy][escape][initial]") {
    // Fixture constants
    const double TIME_STEP = 60.0;
    const double DAYS_TO_SIMULATE = 300.0;
    const double SECONDS_PER_DAY = 86400.0;
    const double AU = 1.496e11;

    // Tolerance constants (precise per observed errors)
    const double V_ESCAPE_TOL = 1e-6;    // velocity match to escape velocity
    const double ECC_TOL = 1e-4;         // eccentricity = 1.0
    const double ENERGY_REL_TOL = 1e-10; // energy relative error
    const double DIST_TOL = 1.0;         // final distance (m) - Python/C++ match to 0.27m
    const double VEL_TOL = 0.001;        // final velocity (m/s) - Python/C++ match to 0.2mm/s
    const double DRIFT_TOL = 1e-12;      // energy drift percent
    const double VEL_DECREASE_TOL = 0.9; // velocity decrease ratio

    SimulationState* sim = create_simulation(10, 0, 0, TIME_STEP);
    REQUIRE(load_system_config(sim, "tests/test_parabolic_orbit.toml"));

    const int COMET_INDEX = 1;
    const int SUN_INDEX = 0;
    CelestialBody* comet = &sim->bodies[COMET_INDEX];
    CelestialBody* sun = &sim->bodies[SUN_INDEX];

    // Initial state
    const double initial_distance = vec3_magnitude(comet->global_position);
    const double initial_velocity = vec3_magnitude(comet->global_velocity);
    const double initial_kinetic = calculate_kinetic_energy(comet);
    const double initial_potential = calculate_potential_energy_pair(comet, sun);
    const double initial_total_energy = initial_kinetic + initial_potential;

    INFO("Initial distance: " << initial_distance / AU << " AU");
    INFO("Initial velocity: " << initial_velocity / 1000.0 << " km/s");
    INFO("Initial kinetic energy: " << initial_kinetic);
    INFO("Initial potential energy: " << initial_potential);
    INFO("Initial total energy: " << initial_total_energy);

    SECTION("velocity matches escape velocity") {
        const double distance = vec3_distance(comet->global_position, sun->global_position);
        const double escape_velocity = sqrt(2.0 * G * sun->mass / distance);
        const double circular_velocity = sqrt(G * sun->mass / distance);

        INFO("Distance: " << distance / AU << " AU");
        INFO("Actual velocity: " << initial_velocity / 1000.0 << " km/s");
        INFO("Escape velocity: " << escape_velocity / 1000.0 << " km/s");
        INFO("Circular velocity: " << circular_velocity / 1000.0 << " km/s");

        const double velocity_error = fabs(initial_velocity - escape_velocity) / escape_velocity;
        INFO("Velocity error from escape velocity: " << velocity_error * 100.0 << "%");
        REQUIRE_THAT(velocity_error, WithinAbs(0.0, V_ESCAPE_TOL));
    }

    SECTION("eccentricity equals 1.0") {
        INFO("Eccentricity: " << comet->orbit.eccentricity);
        REQUIRE_THAT(comet->orbit.eccentricity, WithinAbs(1.0, ECC_TOL));
    }

    SECTION("total energy near zero (relative to KE)") {
        // For a parabolic orbit, total energy should be zero. Due to
        // floating-point cancellation of two large terms (~8.87e22), the
        // absolute value is ~1.68e7 J, but the relative error is ~2e-16.
        const double relative_error = fabs(initial_total_energy) / initial_kinetic;
        INFO("Initial total energy: " << initial_total_energy << " J");
        INFO("Relative error: " << relative_error);
        REQUIRE_THAT(relative_error, WithinAbs(0.0, ENERGY_REL_TOL));
    }

    // Record velocities for trend analysis (every 1000 steps)
    std::vector<double> velocities;
    velocities.push_back(initial_velocity);

    const double max_time = DAYS_TO_SIMULATE * SECONDS_PER_DAY;
    int step_count = 0;
    while (sim->time < max_time) {
        if (step_count % 1000 == 0) {
            velocities.push_back(vec3_magnitude(comet->global_velocity));
        }
        update_simulation(sim);
        step_count++;
    }

    // Final state
    const double final_distance = vec3_magnitude(comet->global_position);
    const double final_velocity = vec3_magnitude(comet->global_velocity);
    const double final_kinetic = calculate_kinetic_energy(comet);
    const double final_potential = calculate_potential_energy_pair(comet, sun);
    const double final_total_energy = final_kinetic + final_potential;

    INFO("Final distance: " << final_distance / AU << " AU");
    INFO("Final velocity: " << final_velocity / 1000.0 << " km/s");
    INFO("Final kinetic energy: " << final_kinetic);
    INFO("Final potential energy: " << final_potential);
    INFO("Final total energy: " << final_total_energy);

    // Precalculated expected values from scripts/precalc_parabolic_orbit.py
    const double expected_distance = 372192353748.3338;  // 2.487917 AU
    const double expected_velocity = 26708.624837;       // 26.708625 km/s

    SECTION("final distance matches escape trajectory") {
        REQUIRE_THAT(final_distance, WithinAbs(expected_distance, DIST_TOL));
    }

    SECTION("final velocity matches escape trajectory") {
        REQUIRE_THAT(final_velocity, WithinAbs(expected_velocity, VEL_TOL));
    }

    SECTION("energy drift near zero") {
        const double energy_drift = fabs(final_total_energy - initial_total_energy);
        const double avg_kinetic = (initial_kinetic + final_kinetic) / 2.0;
        const double drift_pct = (energy_drift / avg_kinetic) * 100.0;

        INFO("Energy drift: " << energy_drift << " J");
        INFO("Energy drift percent: " << drift_pct << "%");
        REQUIRE_THAT(drift_pct, WithinAbs(0.0, DRIFT_TOL));
    }

    SECTION("velocity monotonically decreases (escape trajectory)") {
        int velocity_decreases = 0;
        for (size_t i = 1; i < velocities.size(); i++) {
            if (velocities[i] < velocities[i - 1]) {
                velocity_decreases++;
            }
        }
        const int total_checks = static_cast<int>(velocities.size()) - 1;
        const double decrease_ratio = static_cast<double>(velocity_decreases) / total_checks;

        INFO("Velocity decreases: " << velocity_decreases << " / " << total_checks);
        INFO("Decrease ratio: " << decrease_ratio);
        REQUIRE_THAT(decrease_ratio, WithinAbs(1.0, 1.0 - VEL_DECREASE_TOL));
    }

    destroy_simulation(sim);
}