MeanField/experiments/polytrope_analytic_reference.hpp

#pragma once

// Experiment-only closed-form benchmark. This intentionally does not use the
// production Lane-Emden integration, radial interpolation, EOS, or seed helpers.
#include <cmath>
#include <initializer_list>
#include <numbers>
#include <stdexcept>

namespace experiment::polytrope_validation {
    struct AnalyticValues final {
        double density{};
        double enthalpy{};
        double pressure{};
        double potential{};
        double enclosedMass{};
        // Outward-positive dPhi/dr, not the inward gravitational acceleration.
        double radialPotentialGradient{};
    };

    struct N1Reference final {
        double gravitationalConstant{1.0};
        double mass{1.0};
        double radius{1.0};

        // Validate once before using this reference in a quadrature loop.
        void Validate() const {
            if (!std::isfinite(gravitationalConstant) || gravitationalConstant <= 0.0
                || !std::isfinite(mass) || mass <= 0.0
                || !std::isfinite(radius) || radius <= 0.0) {
                throw std::invalid_argument("The n=1 reference requires finite, positive G, M, and R.");
            }
            for (const double scale : {PolytropicConstant(), CentralDensity(), CentralEnthalpy(),
                                       PressureIntegral(), -BindingEnergy(), MomentOfInertia(),
                                       CentralEnthalpy() / radius, 0.5 * CentralEnthalpy() * CentralDensity()}) {
                if (!std::isfinite(scale) || scale <= 0.0) {
                    throw std::invalid_argument("The n=1 reference scales are not representable as positive finite doubles.");
                }
            }
        }

        [[nodiscard]] double PolytropicConstant() const {
            return 2.0 * gravitationalConstant * radius * radius / std::numbers::pi;
        }

        [[nodiscard]] double CentralDensity() const {
            return (std::numbers::pi / 4.0) * (mass / radius) / radius / radius;
        }

        [[nodiscard]] double CentralEnthalpy() const {
            return gravitationalConstant * mass / radius;
        }

        [[nodiscard]] double PressureIntegral() const {
            return 0.25 * CentralEnthalpy() * mass;
        }

        // W = (1/2) integral rho Phi dV, with Phi tending to zero at infinity.
        [[nodiscard]] double BindingEnergy() const {
            return -0.75 * CentralEnthalpy() * mass;
        }

        // Axial moment of inertia, not integral rho r^2 dV.
        [[nodiscard]] double MomentOfInertia() const {
            constexpr double coefficient = (2.0 / 3.0)
                * (1.0 - 6.0 / (std::numbers::pi * std::numbers::pi));
            return coefficient * mass * radius * radius;
        }

        [[nodiscard]] double DimensionlessTheta(const double physicalRadius) const {
            CheckRadius(physicalRadius);
            if (physicalRadius >= radius) return 0.0;
            const double fraction = physicalRadius / radius;
            const double argument = std::numbers::pi * fraction;
            if (argument < 0.25) {
                const double squared = argument * argument;
                return 1.0 + squared * (-1.0 / 6.0 + squared * (1.0 / 120.0
                    + squared * (-1.0 / 5040.0 + squared * (1.0 / 362880.0
                    + squared * (-1.0 / 39916800.0 + squared / 6227020800.0)))));
            }
            if (fraction > 0.5) {
                // sin(pi-delta) avoids the nonzero floating-point sin(pi)
                // floor. Form the small surface distance before dividing.
                const double surfaceDistance = (radius - physicalRadius) / radius;
                return std::sin(std::numbers::pi * surfaceDistance) / argument;
            }
            return std::sin(argument) / argument;
        }

        // This normalization also accepts numerical potentials: do not clamp
        // its result to the stellar theta range or fit an additive constant.
        [[nodiscard]] double NormalizedPotential(const double potential) const {
            return -potential / CentralEnthalpy() - 1.0;
        }

        [[nodiscard]] AnalyticValues AtRadius(const double physicalRadius) const {
            CheckRadius(physicalRadius);
            const double centralEnthalpy = CentralEnthalpy();
            if (physicalRadius >= radius) {
                const double surfaceFraction = radius / physicalRadius;
                return {
                    .density = 0.0,
                    .enthalpy = 0.0,
                    .pressure = 0.0,
                    .potential = -centralEnthalpy * surfaceFraction,
                    .enclosedMass = mass,
                    .radialPotentialGradient = (centralEnthalpy / radius) * surfaceFraction * surfaceFraction
                };
            }
            const double fraction = physicalRadius / radius;
            const double argument = std::numbers::pi * fraction;
            const double theta = DimensionlessTheta(physicalRadius);
            double massFraction = 0.0;
            double gradientFraction = 0.0;
            if (argument < 0.25) {
                // sin(x)-x*cos(x) = x^3 [1/3-x^2/30+x^4/840-...].
                // Evaluate g separately from m/r^2 to remain regular even
                // when the representable enclosed mass underflows at r~0.
                const double squared = argument * argument;
                const double factor = 1.0 / 3.0 + squared * (-1.0 / 30.0
                    + squared * (1.0 / 840.0 + squared * (-1.0 / 45360.0
                    + squared * (1.0 / 3991680.0 - squared / 518918400.0))));
                massFraction = std::numbers::pi * std::numbers::pi
                    * fraction * fraction * fraction * factor;
                gradientFraction = std::numbers::pi * std::numbers::pi * fraction * factor;
            } else {
                double numerator = 0.0;
                if (fraction > 0.5) {
                    const double delta = std::numbers::pi * ((radius - physicalRadius) / radius);
                    numerator = std::sin(delta) + argument * std::cos(delta);
                } else {
                    numerator = std::sin(argument) - argument * std::cos(argument);
                }
                massFraction = numerator / std::numbers::pi;
                gradientFraction = std::numbers::pi * numerator / (argument * argument);
            }
            const double centralDensity = CentralDensity();
            return {
                .density = centralDensity * theta,
                .enthalpy = centralEnthalpy * theta,
                .pressure = 0.5 * centralEnthalpy * centralDensity * theta * theta,
                .potential = -centralEnthalpy * (1.0 + theta),
                .enclosedMass = mass * massFraction,
                .radialPotentialGradient = (centralEnthalpy / radius) * gradientFraction
            };
        }

    private:
        static void CheckRadius(const double physicalRadius) {
            if (!std::isfinite(physicalRadius) || physicalRadius < 0.0) {
                throw std::invalid_argument("The analytic reference radius must be finite and nonnegative.");
            }
        }
    };
} // namespace experiment::polytrope_validation