#include <array>
#include <cmath>
#include <limits>

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>

import mean_field;
import test_helpers;

TEST_CASE(
    "Polytropic EOS Satisfies Its Analytic Identities",
    tags::barotrope_eos_unit
) {
    using namespace mean_field::eos;

    constexpr double polytropic_index    = 3.0;
    constexpr double polytropic_constant = 1.5;

    const Polytrope barotrope(polytropic_index, polytropic_constant);

    using densityV  = DensityValue;
    using pressureV = PressureValue;
    using enthalpyV = SpecificEnthalpyValue;

    constexpr std::array<densityV, 5> densities{
        densityV{1.0e-6}, densityV{1.0e-3}, densityV{0.1}, densityV{0.7}, densityV{2.0}
    };

    for (const densityV density : densities) {
        const pressureV pressure               = evaluate<quantity::Pressure>(barotrope, density);

        const enthalpyV enthalpy               = evaluate<quantity::SpecificEnthalpy>(barotrope, density);
        const densityV reconstructed_density   = evaluate<quantity::Density>(barotrope, enthalpy);

        const pressureV reconstructed_pressure = evaluate<quantity::Pressure>(barotrope, enthalpy);

        const enthalpyV reconstructed_enthalpy = evaluate<quantity::SpecificEnthalpy>(barotrope, pressure);

        CHECK_THAT(reconstructed_density.value(), Catch::Matchers::WithinRel(density.value(), 2.0e-14));

        CHECK_THAT(reconstructed_pressure.value(), Catch::Matchers::WithinRel(pressure.value(), 2.0e-14));

        CHECK_THAT(reconstructed_enthalpy.value(), Catch::Matchers::WithinRel(enthalpy.value(), 2.0e-14));

        CHECK_THAT(
            pressure.value(),
            Catch::Matchers::WithinRel(density.value() * enthalpy.value() / (polytropic_index + 1.0), 2.0e-14)
        );

        CHECK_THAT(
            (mean_field::eos::partialDerivative<
                 mean_field::eos::quantity::Pressure, mean_field::eos::quantity::SpecificEnthalpy>(
                 barotrope, mean_field::eos::SpecificEnthalpyValue{enthalpy}
            )
                 .value()),
            Catch::Matchers::WithinRel(density.value(), 2.0e-14)
        );

        CHECK_THAT(
            (mean_field::eos::partialDerivative<
                 mean_field::eos::quantity::Pressure, mean_field::eos::quantity::Density>(
                 barotrope, mean_field::eos::DensityValue{density}
            )
                 .value()),
            Catch::Matchers::WithinRel(enthalpy.value() / polytropic_index, 2.0e-14)
        );
    }
}

TEST_CASE(
    "Polytropic EOS Derivatives Match Centered Differences",
    tags::barotrope_eos_jacobian
) {
    using namespace mean_field::eos;

    const Polytrope barotrope(3.0, 1.5);

    using densityV  = DensityValue;
    using pressureV = PressureValue;
    using enthalpyV = SpecificEnthalpyValue;

    constexpr std::array<enthalpyV, 4> enthalpies{enthalpyV{0.05}, enthalpyV{0.2}, enthalpyV{0.7}, enthalpyV{1.4}};

    for (const enthalpyV enthalpy : enthalpies) {
        const enthalpyV step                = enthalpyV{1.0e-6} * std::max(1.0, enthalpy.value());

        const densityV density_difference   = (evaluate<quantity::Density>(barotrope, enthalpy + step) -
                                               evaluate<quantity::Density>(barotrope, enthalpy - step)) /
                                              (2.0 * step.value());

        const pressureV pressure_difference = (evaluate<quantity::Pressure>(barotrope, enthalpy + step) -
                                               evaluate<quantity::Pressure>(barotrope, enthalpy - step)) /
                                              (2.0 * step.value());

        CHECK_THAT(
            density_difference.value(),
            Catch::Matchers::WithinRel(
                mean_field::eos::partialDerivative<
                    mean_field::eos::quantity::Density, mean_field::eos::quantity::SpecificEnthalpy>(
                    barotrope, enthalpy
                )
                    .value(),
                5.0e-10
            )
        );

        CHECK_THAT(
            pressure_difference.value(),
            Catch::Matchers::WithinRel(
                mean_field::eos::partialDerivative<
                    mean_field::eos::quantity::Pressure, mean_field::eos::quantity::SpecificEnthalpy>(
                    barotrope, enthalpy
                )
                    .value(),
                5.0e-10
            )
        );
    }
}

TEST_CASE(
    "Polytropic EOS Has An Exact Zero Density Surface",
    tags::barotrope_eos_unit
) {
    const mean_field::eos::Polytrope barotrope(3.0, 1.5);

    CHECK(
        mean_field::eos::evaluate<mean_field::eos::quantity::Density>(
            barotrope, mean_field::eos::SpecificEnthalpyValue{-1.0}
        )
            .value() == 0.0
    );
    CHECK(
        mean_field::eos::evaluate<mean_field::eos::quantity::Density>(
            barotrope, mean_field::eos::SpecificEnthalpyValue{0.0}
        )
            .value() == 0.0
    );

    CHECK(
        mean_field::eos::evaluate<mean_field::eos::quantity::Pressure>(
            barotrope, mean_field::eos::SpecificEnthalpyValue{-1.0}
        )
            .value() == 0.0
    );
    CHECK(
        mean_field::eos::evaluate<mean_field::eos::quantity::Pressure>(
            barotrope, mean_field::eos::SpecificEnthalpyValue{0.0}
        )
            .value() == 0.0
    );

    CHECK(
        (mean_field::eos::partialDerivative<
             mean_field::eos::quantity::Density, mean_field::eos::quantity::SpecificEnthalpy>(
             barotrope, mean_field::eos::SpecificEnthalpyValue{-1.0}
         )
             .value() == 0.0)
    );

    CHECK(
        (mean_field::eos::partialDerivative<
             mean_field::eos::quantity::Density, mean_field::eos::quantity::SpecificEnthalpy>(
             barotrope, mean_field::eos::SpecificEnthalpyValue{0.0}
         )
             .value() == 0.0)
    );

    CHECK(
        (mean_field::eos::partialDerivative<
             mean_field::eos::quantity::Pressure, mean_field::eos::quantity::SpecificEnthalpy>(
             barotrope, mean_field::eos::SpecificEnthalpyValue{0.0}
         )
             .value() == 0.0)
    );
}

TEST_CASE(
    "Polytropic EOS Rejects Invalid Material Parameters",
    tags::barotrope_eos_unit
) {
    CHECK_THROWS_AS(mean_field::eos::Polytrope(0.5, 1.0), std::invalid_argument);

    CHECK_THROWS_AS(mean_field::eos::Polytrope(3.0, 0.0), std::invalid_argument);

    CHECK_THROWS_AS(mean_field::eos::Polytrope(std::numeric_limits<double>::infinity(), 1.0), std::invalid_argument);

    const mean_field::eos::Polytrope barotrope(3.0, 1.0);

    CHECK_THROWS_AS(
        mean_field::eos::evaluate<mean_field::eos::quantity::Pressure>(barotrope, mean_field::eos::DensityValue{-1.0}),
        std::domain_error
    );

    CHECK_THROWS_AS(
        mean_field::eos::evaluate<mean_field::eos::quantity::SpecificEnthalpy>(
            barotrope, mean_field::eos::DensityValue{-1.0}
        ),
        std::domain_error
    );

    CHECK_THROWS_AS(
        mean_field::eos::evaluate<mean_field::eos::quantity::SpecificEnthalpy>(
            barotrope, mean_field::eos::PressureValue{-1.0}
        ),
        std::domain_error
    );
}