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

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

import mean_field;
import test_helpers;

namespace eos = mean_field::eos;

TEST_CASE(
    "Polytropic EOS Pressure To Specific Enthalpy Relation Is Characterized",
    tags::polytropic_eos_characterization
) {
    constexpr std::array<double, 3> polytropicIndices{1.0, 1.5, 3.0};
    constexpr std::array<double, 3> polytropicConstants{0.25, 0.73, 2.0};
    constexpr std::array<double, 5> pressures{0.0, 1.0e-12, 1.0e-4, 0.3, 5.0};

    for (const double polytropicIndex : polytropicIndices) {
        for (const double polytropicConstant : polytropicConstants) {
            const mean_field::eos::Polytrope equationOfState(polytropicIndex, polytropicConstant);

            for (const double pressure : pressures) {
                CAPTURE(polytropicIndex, polytropicConstant, pressure);

                const double indexPlusOne     = polytropicIndex + 1.0;
                const double expectedEnthalpy = indexPlusOne *
                                                std::pow(polytropicConstant, polytropicIndex / indexPlusOne) *
                                                std::pow(pressure, 1.0 / indexPlusOne);

                const double enthalpy =
                    eos::evaluate<eos::quantity::SpecificEnthalpy>(equationOfState, eos::PressureValue{pressure})
                        .value();

                if (pressure == 0.0) {
                    CHECK(enthalpy == 0.0);
                } else {
                    CHECK_THAT(enthalpy, Catch::Matchers::WithinRel(expectedEnthalpy, 5.0e-14));

                    const double recoveredPressure =
                        eos::evaluate<eos::quantity::Pressure>(equationOfState, eos::SpecificEnthalpyValue{enthalpy})
                            .value();

                    CHECK_THAT(recoveredPressure, Catch::Matchers::WithinRel(pressure, 5.0e-13));
                }
            }
        }
    }
}

TEST_CASE(
    "Polytropic EOS Domain Contract Covers Every Relation",
    tags::polytropic_eos_characterization
) {
    constexpr double infinity = std::numeric_limits<double>::infinity();
    constexpr double quietNaN = std::numeric_limits<double>::quiet_NaN();

    for (const double invalidIndex : std::array<double, 4>{0.999, infinity, -infinity, quietNaN}) {
        CAPTURE(invalidIndex);
        CHECK_THROWS_AS(mean_field::eos::Polytrope(invalidIndex, 1.0), std::invalid_argument);
    }

    for (const double invalidConstant : std::array<double, 5>{0.0, -0.1, infinity, -infinity, quietNaN}) {
        CAPTURE(invalidConstant);
        CHECK_THROWS_AS(mean_field::eos::Polytrope(3.0, invalidConstant), std::invalid_argument);
    }

    const mean_field::eos::Polytrope equationOfState(3.0, 0.75);

    constexpr double negativeDensity = -0.1;
    CHECK_THROWS_AS(
        eos::evaluate<eos::quantity::Pressure>(equationOfState, eos::DensityValue{negativeDensity}), std::domain_error
    );
    CHECK_THROWS_AS(
        eos::evaluate<eos::quantity::SpecificEnthalpy>(equationOfState, eos::DensityValue{negativeDensity}),
        std::domain_error
    );
    CHECK_THROWS_AS(
        (eos::partialDerivative<eos::quantity::Pressure, eos::quantity::Density>(
            equationOfState, eos::DensityValue{negativeDensity}
        )),
        std::domain_error
    );

    CHECK_THROWS_AS(
        eos::evaluate<eos::quantity::SpecificEnthalpy>(equationOfState, eos::PressureValue{-0.1}), std::domain_error
    );

    constexpr double exteriorEnthalpy = -0.1;
    CHECK(
        eos::evaluate<eos::quantity::Density>(equationOfState, eos::SpecificEnthalpyValue{exteriorEnthalpy}).value() ==
        0.0
    );
    CHECK(
        eos::evaluate<eos::quantity::Pressure>(equationOfState, eos::SpecificEnthalpyValue{exteriorEnthalpy}).value() ==
        0.0
    );
    CHECK(
        (eos::partialDerivative<eos::quantity::Density, eos::quantity::SpecificEnthalpy>(
             equationOfState, eos::SpecificEnthalpyValue{exteriorEnthalpy}
         )
             .value() == 0.0)
    );
    CHECK(
        (eos::partialDerivative<eos::quantity::Pressure, eos::quantity::SpecificEnthalpy>(
             equationOfState, eos::SpecificEnthalpyValue{exteriorEnthalpy}
         )
             .value() == 0.0)
    );

    for (const double nonfiniteValue : std::array<double, 3>{infinity, -infinity, quietNaN}) {
        CAPTURE(nonfiniteValue);

        CHECK_THROWS_AS(
            eos::evaluate<eos::quantity::Pressure>(equationOfState, eos::DensityValue{nonfiniteValue}),
            std::domain_error
        );
        CHECK_THROWS_AS(
            eos::evaluate<eos::quantity::SpecificEnthalpy>(equationOfState, eos::DensityValue{nonfiniteValue}),
            std::domain_error
        );
        CHECK_THROWS_AS(
            (eos::partialDerivative<eos::quantity::Pressure, eos::quantity::Density>(
                equationOfState, eos::DensityValue{nonfiniteValue}
            )),
            std::domain_error
        );
        CHECK_THROWS_AS(
            eos::evaluate<eos::quantity::SpecificEnthalpy>(equationOfState, eos::PressureValue{nonfiniteValue}),
            std::domain_error
        );
        CHECK_THROWS_AS(
            eos::evaluate<eos::quantity::Density>(equationOfState, eos::SpecificEnthalpyValue{nonfiniteValue}),
            std::domain_error
        );
        CHECK_THROWS_AS(
            eos::evaluate<eos::quantity::Pressure>(equationOfState, eos::SpecificEnthalpyValue{nonfiniteValue}),
            std::domain_error
        );
        CHECK_THROWS_AS(
            (eos::partialDerivative<eos::quantity::Density, eos::quantity::SpecificEnthalpy>(
                equationOfState, eos::SpecificEnthalpyValue{nonfiniteValue}
            )),
            std::domain_error
        );
        CHECK_THROWS_AS(
            (eos::partialDerivative<eos::quantity::Pressure, eos::quantity::SpecificEnthalpy>(
                equationOfState, eos::SpecificEnthalpyValue{nonfiniteValue}
            )),
            std::domain_error
        );
    }
}