#include <numbers>
#include <print>

#include "serif/eos/models/polytropic.hpp"
#include "serif/equations/angular_momentum.hpp"
#include "serif/equations/barotropic_density.hpp"
#include "serif/equations/surface_potential.hpp"
#include "serif/equations/mass_normalization.hpp"
#include "serif/equations/primal_gravity.hpp"
#include "serif/equations/registry.hpp"
#include "serif/model/description.hpp"
#include "serif/model/eos_reference.hpp"
#include "serif/geometry/surface/isobaric.hpp"

namespace
{
    namespace model = serif::model;
    namespace equations = serif::equations;
    namespace dimensions = serif::dimensions;
    namespace surface = serif::geometry::surface;
}

/** @brief Inspect an explicitly chosen rotating equilibrium formulation without creating a mesh or solving it. */
int main()
{
    // This list is the physical formulation. The EOS supplies relations; it does not choose equations.
    const auto equilibrium = [](const auto& physics)
    {
        const auto eos = model::borrow_eos(physics.equation_of_state());
        const auto& rotation = physics.template specification<model::FixedAngularMomentum>();

        return model::EquilibriumEquations
        {
            .interior_equations = equations::EquationRegistry
            {
                equations::BarotropicDensityBalance{eos, dimensions::PressureValue{0.0}, rotation.axis(), rotation.center()},
                equations::PrimalPoissonBalance{4.0 * std::numbers::pi},
                equations::MassNormalization{physics.template specification<model::FixedTotalMass>().target().value()},
                equations::AngularMomentumConstraint{rotation.target().value(), rotation.axis(), rotation.center()}
            },
            .shape_force = equations::SurfacePotentialBalance{rotation.axis(), rotation.center()},
            .surface_condition = surface::Isobaric{dimensions::PressureValue{0.0}}
        };
    };

    // Mass and angular momentum are prescribed. Density, potential, surface level, angular velocity and shape are solved.
    const auto star = model::StellarModel
    {
        serif::eos::models::Polytrope{3.0, 0.25},
        equilibrium,
        model::FixedTotalMass{dimensions::MassValue{1.0}},
        model::FixedAngularMomentum{dimensions::AngularMomentumValue{0.02}}
    };

    // In this example G = 1 in the chosen physical units, and scalar Poisson balance supplies gravity.
    // Rigid rotation uses the same axis and origin in the density equation and on the pressure surface.
    // Resolving reads immutable physics once. Keep star alive while the resolved equations borrow its EOS.
    const auto chosen = model::resolve_equilibrium_equations(star);
    std::println("Interior equations:");
    chosen.interior_equations.for_each([](const auto& equation)
    {
        std::println("  {}", equation.name);
    });
    std::println("Shape equation: {}", chosen.shape_force.name);
    std::println("Surface pressure: {}", chosen.surface_condition.target_pressure().value());
    std::println("Prescribed mass: {}; angular momentum: {}", star.specification<model::FixedTotalMass>().target().value(), star.specification<model::FixedAngularMomentum>().target().value());
}