#include <cmath>
#include <limits>

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

#include "serif/equations/barotropic_density.hpp"
#include "serif/solver/normalization/scales.hpp"
#include "serif/tests/equations/barotropic_eos.hpp"

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

    using MixedPowerBarotrope = serif::tests::equations::MixedPowerBarotrope;
    using Balance = serif::equations::BarotropicDensityBalance<MixedPowerBarotrope>;
    using Point = serif::equations::EquationPoint<Balance::Operands>;

    /** @brief Check one scalar physical column independently of the other field and geometry columns. */
    template <typename Term>
    void check_field_derivative(const Balance& equation, const Point& point) {
        constexpr double step = 2.0e-6;
        Point positive = point, negative = point;
        positive.fields.get<Term>() += step;
        negative.fields.get<Term>() -= step;
        const double difference = (equation.evaluate(positive).residual - equation.evaluate(negative).residual) / (2.0 * step);
        CHECK_THAT(equation.evaluate(point).field_derivatives.get<Term>(), Catch::Matchers::WithinAbs(difference, 2.0e-8));
    }
}

TEST_CASE("Barotropic density balance preserves finite pressure and primitive gauge invariance", "[equations][barotropic_density][unit]") {
    const Balance equation{MixedPowerBarotrope{}, dimensions::PressureValue{2.0}};
    const Balance shifted{MixedPowerBarotrope{.gauge = 7.5}, dimensions::PressureValue{2.0}};
    Point point;
    point.fields.get<Balance::Density>() = 2.0;
    point.fields.get<Balance::Potential>() = -12.0;
    point.fields.get<Balance::SurfaceLevel>() = -2.0 / 3.0;
    const auto balance = equation.evaluate(point);
    const auto shifted_balance = shifted.evaluate(point);
    CHECK_THAT(balance.residual, Catch::Matchers::WithinAbs(0.0, 3.0e-15));
    CHECK_THAT(shifted_balance.residual, Catch::Matchers::WithinAbs(balance.residual, 4.0e-15));
    CHECK(shifted_balance.field_derivatives.values == balance.field_derivatives.values);
    CHECK(shifted_balance.position_derivative == balance.position_derivative);
    CHECK(equation.surface_pressure().value() == 2.0);

    const Balance vacuum{MixedPowerBarotrope{}};
    CHECK_THAT(vacuum.evaluate(point).residual, Catch::Matchers::WithinAbs(10.0 / 3.0, 3.0e-15));
    point.fields.get<Balance::Density>() = 1.0;
    point.fields.get<Balance::Potential>() = -2.0;
    point.fields.get<Balance::SurfaceLevel>() = -2.0;
    CHECK(equation.evaluate(point).residual == 0.0);
    CHECK(shifted.evaluate(point).residual == 0.0);
}

TEST_CASE("Barotropic density tangent includes every field and physical position", "[equations][barotropic_density][unit]") {
    const Balance equation{MixedPowerBarotrope{}, dimensions::PressureValue{2.0}, {1.0, 2.0, 3.0}, {0.2, -0.1, 0.3}};
    Point point;
    point.fields.get<Balance::Density>() = 1.3;
    point.fields.get<Balance::Potential>() = -3.2;
    point.fields.get<Balance::SurfaceLevel>() = -0.9;
    point.fields.get<Balance::AngularVelocity>() = 0.8;
    point.position = {1.2, 0.7, -0.2};
    // The displacement from the rotation origin is (1, 0.8, -0.5), with projection 1.1/sqrt(14).
    const double centrifugal_potential = 0.5 * (1.0 + 0.64 + 0.25 - 1.21 / 14.0);
    const double expected = 2.6 + (4.0 / 3.0) * 1.3 * 1.3 * 1.3 - 10.0 / 3.0 - 3.2 - 0.64 * centrifugal_potential + 0.9;
    const auto balance = equation.evaluate(point);
    CHECK_THAT(balance.residual, Catch::Matchers::WithinAbs(expected, 4.0e-15));
    CHECK_THAT(balance.field_derivatives.get<Balance::Density>(), Catch::Matchers::WithinAbs(2.0 + 4.0 * 1.3 * 1.3, 2.0e-15));
    check_field_derivative<Balance::Density>(equation, point);
    check_field_derivative<Balance::Potential>(equation, point);
    check_field_derivative<Balance::SurfaceLevel>(equation, point);
    check_field_derivative<Balance::AngularVelocity>(equation, point);
    constexpr double step = 2.0e-6;
    for (int axis = 0; axis < 3; ++axis) {
        Point positive = point, negative = point;
        positive.position[axis] += step;
        negative.position[axis] -= step;
        const double difference = (equation.evaluate(positive).residual - equation.evaluate(negative).residual) / (2.0 * step);
        CHECK_THAT(balance.position_derivative[axis], Catch::Matchers::WithinAbs(difference, 2.0e-8));
    }
    point.fields.get<Balance::AngularVelocity>() = 0.0;
    CHECK(equation.evaluate(point).position_derivative == serif::discretization::forms::VectorValue{});
}

TEST_CASE("Barotropic density balance rejects invalid interior states without clipping", "[equations][barotropic_density][unit]") {
    const Balance equation{MixedPowerBarotrope{}};
    Point point;
    for (const double density : {0.0, -0.1, std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()}) {
        point.fields.get<Balance::Density>() = density;
        CHECK_THROWS_AS(equation.evaluate(point), eos::EOSEvaluationError);
    }
    point.fields.get<Balance::Density>() = 1.0e200;
    CHECK_THROWS_AS(equation.evaluate(point), eos::EOSEvaluationError);
    for (const double pressure : {-0.1, std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()}) {
        CHECK_THROWS_AS((Balance{MixedPowerBarotrope{}, dimensions::PressureValue{pressure}}), eos::EOSEvaluationError);
    }
    CHECK_THROWS_AS((Balance{MixedPowerBarotrope{.gauge = std::numeric_limits<double>::infinity()}}), eos::EOSEvaluationError);
}

TEST_CASE("Pressure potential uses the gravitational specific energy reference scale", "[equations][barotropic_density][unit]") {
    const auto scales = serif::solver::normalization::prepare_stellar_reference_scales(dimensions::MassValue{2.0},
        dimensions::LengthValue{4.0}, dimensions::QuantityValue<dimensions::GravitationalConstant>{3.0});
    REQUIRE(scales.has_value());
    CHECK(serif::solver::normalization::characteristic_scale(*scales, dimensions::BarotropicPressurePotential{}).value() == 1.5);
}