// Ported from MeanField tests/physics/polytropic_eos_characterization.cpp.

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

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

#include "serif/eos/evaluation.hpp"
#include "serif/eos/exceptions.hpp"
#include "serif/eos/models/polytropic.hpp"

#include "serif/tests/test_tags.hpp"

namespace eos = serif::eos;
namespace dimensions = serif::dimensions;

TEST_CASE(
    "Polytropic EOS Pressure To Specific Enthalpy Relation Is Characterized",
    tags::unit & tags::eos
) {
    constexpr std::array<double, 3> polytropic_indices{1.0, 1.5, 3.0};
    constexpr std::array<double, 3> polytropic_constants{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 polytropic_index : polytropic_indices) {
        for (const double polytropic_constant : polytropic_constants) {
            const eos::models::Polytrope equation_of_state(polytropic_index, polytropic_constant);

            for (const double pressure : pressures) {
                CAPTURE(polytropic_index, polytropic_constant, pressure);

                const double index_plus_one = polytropic_index + 1.0;
                const double expected_enthalpy = index_plus_one *
                                                 std::pow(polytropic_constant, polytropic_index / index_plus_one) *
                                                 std::pow(pressure, 1.0 / index_plus_one);

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

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

                    const double recovered_pressure = eos::evaluate<dimensions::Pressure>(
                        equation_of_state, dimensions::SpecificEnthalpyValue{enthalpy}
                    )
                                                          .value();

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

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

    for (const double invalid_index : std::array<double, 4>{0.999, infinity, -infinity, quiet_nan}) {
        CAPTURE(invalid_index);
        CHECK_THROWS_AS(eos::models::Polytrope(invalid_index, 1.0), eos::EOSEvaluationError);
    }

    for (const double invalid_constant : std::array<double, 5>{0.0, -0.1, infinity, -infinity, quiet_nan}) {
        CAPTURE(invalid_constant);
        CHECK_THROWS_AS(eos::models::Polytrope(3.0, invalid_constant), eos::EOSEvaluationError);
    }

    const eos::models::Polytrope equation_of_state(3.0, 0.75);

    constexpr double negative_density = -0.1;
    CHECK_THROWS_AS(
        eos::evaluate<dimensions::Pressure>(equation_of_state, dimensions::DensityValue{negative_density}),
        eos::EOSEvaluationError
    );
    CHECK_THROWS_AS(
        eos::evaluate<dimensions::SpecificEnthalpy>(equation_of_state, dimensions::DensityValue{negative_density}),
        eos::EOSEvaluationError
    );
    CHECK_THROWS_AS(
        (eos::partial_derivative<dimensions::Pressure, dimensions::MassDensity>(
            equation_of_state, dimensions::DensityValue{negative_density}
        )),
        eos::EOSEvaluationError
    );

    CHECK_THROWS_AS(
        eos::evaluate<dimensions::SpecificEnthalpy>(equation_of_state, dimensions::PressureValue{-0.1}),
        eos::EOSEvaluationError
    );

    constexpr double exterior_enthalpy = -0.1;
    CHECK(
        eos::evaluate<dimensions::MassDensity>(equation_of_state, dimensions::SpecificEnthalpyValue{exterior_enthalpy})
            .value() == 0.0
    );
    CHECK(
        eos::evaluate<dimensions::Pressure>(equation_of_state, dimensions::SpecificEnthalpyValue{exterior_enthalpy})
            .value() == 0.0
    );
    CHECK(
        (eos::partial_derivative<dimensions::MassDensity, dimensions::SpecificEnthalpy>(
             equation_of_state, dimensions::SpecificEnthalpyValue{exterior_enthalpy}
         )
                .value() == 0.0)
    );
    CHECK(
        (eos::partial_derivative<dimensions::Pressure, dimensions::SpecificEnthalpy>(
             equation_of_state, dimensions::SpecificEnthalpyValue{exterior_enthalpy}
         )
                .value() == 0.0)
    );

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

        CHECK_THROWS_AS(
            eos::evaluate<dimensions::Pressure>(equation_of_state, dimensions::DensityValue{nonfinite_value}),
            eos::EOSEvaluationError
        );
        CHECK_THROWS_AS(
            eos::evaluate<dimensions::SpecificEnthalpy>(equation_of_state, dimensions::DensityValue{nonfinite_value}),
            eos::EOSEvaluationError
        );
        CHECK_THROWS_AS(
            (eos::partial_derivative<dimensions::Pressure, dimensions::MassDensity>(
                equation_of_state, dimensions::DensityValue{nonfinite_value}
            )),
            eos::EOSEvaluationError
        );
        CHECK_THROWS_AS(
            eos::evaluate<dimensions::SpecificEnthalpy>(equation_of_state, dimensions::PressureValue{nonfinite_value}),
            eos::EOSEvaluationError
        );
        CHECK_THROWS_AS(
            eos::evaluate<dimensions::MassDensity>(
                equation_of_state, dimensions::SpecificEnthalpyValue{nonfinite_value}
            ),
            eos::EOSEvaluationError
        );
        CHECK_THROWS_AS(
            eos::evaluate<dimensions::Pressure>(equation_of_state, dimensions::SpecificEnthalpyValue{nonfinite_value}),
            eos::EOSEvaluationError
        );
        CHECK_THROWS_AS(
            (eos::partial_derivative<dimensions::MassDensity, dimensions::SpecificEnthalpy>(
                equation_of_state, dimensions::SpecificEnthalpyValue{nonfinite_value}
            )),
            eos::EOSEvaluationError
        );
        CHECK_THROWS_AS(
            (eos::partial_derivative<dimensions::Pressure, dimensions::SpecificEnthalpy>(
                equation_of_state, dimensions::SpecificEnthalpyValue{nonfinite_value}
            )),
            eos::EOSEvaluationError
        );
    }
}