Program Listing for File Galaxy.hpp

↰ Return to documentation for file (src/core/Galaxy.hpp)

#ifndef GALAXY_HPP
#define GALAXY_HPP

#include "../io/SimControls.hpp"
#include "../pdfs/PDF.hpp"
#include "../pdfs/PDFReflect.hpp"
#include "../pdfs/PDFSegmentPowerlaw.hpp"
#include "../specsyn/Specsyn.hpp"
#include "../utils/PDFIntegrator.hpp"
#include "Cluster.hpp"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iterator>
#include <memory>
#include <numbers>
#include <optional>
#include <variant>
#include <vector>

namespace extinct
{
    class Extinct;
} // namespace extinct

namespace nebular
{
    class Nebular;
} // namespace nebular

namespace core
{

    class Galaxy
    {
    public:

        using SFRVar = std::variant<pdfs::PDF, std::reference_wrapper<const pdfs::PDF>>;

        struct FieldStar
        {
            double mass_;
            double feh_;
            double formTime_;
            double deathTime_;
            double aV_;
        };

        explicit Galaxy(const io::SimControls& controls);

        // Observers

        [[nodiscard]] auto sfr() const -> const pdfs::PDF&
        {
            if (const auto* owned = std::get_if<pdfs::PDF>(&sfr_)) { return *owned; }
            return std::get<std::reference_wrapper<const pdfs::PDF>>(sfr_).get();
        }

        [[nodiscard]] auto curTime() const { return curTime_; }

        [[nodiscard]] auto clusters() -> auto& { return clusters_; }

        [[nodiscard]] auto disruptedClusters() -> auto& { return disruptedClusters_; }

        [[nodiscard]] auto fieldStars() const -> const auto& { return fieldStars_; }

        [[nodiscard]] auto deadFieldStars() const -> const auto& { return deadFieldStars_; }

        [[nodiscard]] auto spec() -> const auto&
        {
            if (!specCurrent_) { computeSpec(); specCurrent_ = true; }
            return spec_;
        }

        [[nodiscard]] auto specExtinct() -> const auto&
        {
            if (!specCurrent_) { computeSpec(); specCurrent_ = true; }
            return specExtinct_;
        }

        [[nodiscard]] auto phot() -> const auto&
        {
            if (!photCurrent_) { computePhot(); photCurrent_ = true; }
            return phot_;
        }

        [[nodiscard]] auto photExtinct() -> const auto&
        {
            if (!photCurrent_) { computePhot(); photCurrent_ = true; }
            return photExtinct_;
        }

        [[nodiscard]] auto specNeb() -> const auto&
        {
            if (!specCurrent_) { computeSpec(); specCurrent_ = true; }
            return specNeb_;
        }

        [[nodiscard]] auto specNebExtinct() -> const auto&
        {
            if (!specCurrent_) { computeSpec(); specCurrent_ = true; }
            return specNebExtinct_;
        }

        [[nodiscard]] auto lineLum() -> const auto&
        {
            if (!specCurrent_) { computeSpec(); specCurrent_ = true; }
            return lineLum_;
        }

        [[nodiscard]] auto lineLumExtinct() -> const auto&
        {
            if (!specCurrent_) { computeSpec(); specCurrent_ = true; }
            return lineLumExtinct_;
        }

        [[nodiscard]] auto photNeb() -> const auto&
        {
            if (!photCurrent_) { computePhot(); photCurrent_ = true; }
            return photNeb_;
        }

        [[nodiscard]] auto photNebExtinct() -> const auto&
        {
            if (!photCurrent_) { computePhot(); photCurrent_ = true; }
            return photNebExtinct_;
        }

        [[nodiscard]] auto lbol()
        {
            if (!lbolCurrent_) { computeLbol(); lbolCurrent_ = true; }
            return lbol_;
        }

        [[nodiscard]] auto yields() const -> const auto&
        {
            return yields_;
        }

        [[nodiscard]] auto cumSNe() const -> double
        {
            double nSN = static_cast<double>(fieldStarStochSNe_) + fieldStarNonStochSNe_;
            for (const auto& cluster : clusters_) { nSN += cluster.cumSNe(); }
            for (const auto& cluster : disruptedClusters_) { nSN += cluster.cumSNe(); }
            return nSN;
        }

        [[nodiscard]] auto mDotWind() -> double
        {
            if (!windsCurrent_) { computeWinds(); windsCurrent_ = true; }
            return mDotWind_;
        }

        [[nodiscard]] auto pDotWind() -> double
        {
            if (!windsCurrent_) { computeWinds(); windsCurrent_ = true; }
            return pDotWind_;
        }

        [[nodiscard]] auto eDotWind() -> double
        {
            if (!windsCurrent_) { computeWinds(); windsCurrent_ = true; }
            return eDotWind_;
        }

        [[nodiscard]] auto yieldsRate(double t, double feh) const -> std::vector<double>;

        [[nodiscard]] auto yieldsRate(double t) const -> std::vector<double>;

        struct NoDeathRateHook
        {
            void operator()(double /*t*/, std::vector<double>& /*rate*/) const {}
        };

        template <class F, class G = NoDeathRateHook>
        [[nodiscard]] auto continuousDeathQuantity(const F& f, const std::size_t n, const double tLast,
            const double tNow, const G& g = G{}) const -> std::vector<double>
        {
            const auto& sc = controls_.get();
            std::vector<double> result(n, 0.0);
            if (tNow <= tLast || sc.minStochMass() <= 0.0 || sc.fCluster() >= 1.0) { return result; }

            const double dt = tNow - tLast;
            const pdfs::PDF step(std::make_unique<pdfs::PDFSegmentPowerlaw>(tLast, tNow, 0.0));
            const auto rateAt = [&](const double t) -> std::vector<double>
            { return deathRate(t, f, n, g, std::nullopt); };
            result = utils::integrateScaled(step, rateAt, n, tLast, tNow,
                { tLast + (0.25 * dt), tLast + (0.5 * dt), tLast + (0.75 * dt), tNow },
                sc.intMaxIter(), sc.intAbsTol(), sc.intRelTol());
            const double factor = (1.0 - sc.fCluster()) * dt / step.integral(tLast, tNow);
            for (double& r : result) { r *= factor; }
            return result;
        }

        [[nodiscard]] auto targetMass() const { return targetMass_; }

        [[nodiscard]] auto actualMass() const { return actualMass_; }

        void advance(double t);

    private:

        template <class F>
        [[nodiscard]] auto integrateCtsIntegrand(const double age, const double feh, const F& f,
            const std::size_t nInt) const -> std::vector<double>
        {
            const auto& sc = controls_.get();
            const auto& imf = sc.imf();
            const double logAge = std::max(std::log10(age), sc.tracks()->logTMin());
            const auto isochrone = sc.constFeH() ?
                sc.tracks2D()->getIsochrone(logAge) : sc.tracks()->getIsochrone(logAge, feh);

            std::vector<double> result(nInt, 0.0);
            for (const auto& seg : isochrone)
            {
                const double a = std::max(imf.getMin(), seg->xMin());
                const double b = std::min(sc.minStochMass(), seg->xMax());
                if (a >= b) { continue; } // empty intersection with [imf.getMin(), minStochMass()]

                const auto fAt = [&f, &seg, feh](const double m) -> std::vector<double>
                {
                    const auto v = f(m, *seg, feh);
                    std::vector<double> out(std::begin(v), std::end(v));
                    return out;
                };
                const auto segResult = utils::integrateScaled(imf, fAt, nInt, a, b,
                    { a, std::sqrt(a * b), b }, sc.intMaxIter(), sc.intAbsTol(), sc.intRelTol());
                for (std::size_t k = 0; k < nInt; ++k) { result.at(k) += segResult.at(k); }
            }
            return result;
        }

        template <class F>
        [[nodiscard]] auto integrateCts(const F& f, const std::size_t nInt) const -> std::vector<double>
        {
            const auto& sc = controls_.get();
            const double massInt = sc.nonStochIMFMass();
            if (!(curTime_ > 0.0) || !(massInt > 0.0)) { return std::vector<double>(nInt, 0.0); } // NOLINT(modernize-return-braced-init-list) -- braces would select vector's initializer_list constructor, giving {double(nInt), 0.0}

            // See Specsyn::specCtsHelper()'s own comment for why this
            // reflects sfr() about curTime_ / 2 (so this layer's own
            // coordinate is age directly), why ageMin is
            // min(1e4 yr, 1e-3 curTime_), and why it is log-transformed
            const pdfs::PDFReflect sfrAge(sfr(), 0.5 * curTime_);
            const double ageMin = std::min(1e4, 1e-3 * curTime_);
            constexpr int nAgeScale = 5;
            std::vector<double> ageScalePoints;
            ageScalePoints.reserve(nAgeScale);
            for (int i = 0; i < nAgeScale; ++i)
            {
                ageScalePoints.push_back(ageMin * std::pow(curTime_ / ageMin,
                    static_cast<double>(i) / (nAgeScale - 1)));
            }

            const auto ageIntegral = [&](const double feh) -> std::vector<double>
            {
                const auto atAge = [&, feh](const double age) -> std::vector<double>
                { return integrateCtsIntegrand(age, feh, f, nInt); };
                return utils::integrateScaled(sfrAge, atAge, nInt, ageMin, curTime_, ageScalePoints,
                    sc.intMaxIter(), sc.intAbsTol(), sc.intRelTol(), true);
            };

            const auto& fehDist = sc.fehDist();
            std::vector<double> result;
            if (fehDist.getMin() == fehDist.getMax())
            {
                result = ageIntegral(fehDist.getMin());
            }
            else
            {
                result = utils::integrateScaled(fehDist, ageIntegral, nInt, fehDist.getMin(),
                    fehDist.getMax(), { fehDist.expectationValue() }, sc.intMaxIter(),
                    sc.intAbsTol(), sc.intRelTol());
                const double fehNorm = fehDist.integral(fehDist.getMin(), fehDist.getMax());
                for (double& r : result) { r /= fehNorm; }
            }

            // imf() is normalized by number, so result is per star
            // formed; convert to this population's own total -- see
            // SimControls::nonStochIMFMass()'s own comment
            const double scale = (1.0 - sc.fCluster()) * (1.0 - sc.fracStochMass()) / massInt;
            for (double& r : result) { r *= scale; }
            return result;
        }

        template <class F>
        [[nodiscard]] auto deathIntegrand(const double age, const double feh, const F& f,
            const std::size_t n) const -> std::vector<double>
        {
            const auto& sc = controls_.get();
            std::vector<double> result(n, 0.0);
            const auto massDeriv = sc.constFeH() ?
                sc.tracks2D()->massAndDerivFromLifetime(std::log10(age)) :
                sc.tracks()->massAndDerivFromLifetime(std::log10(age), feh);
            for (const auto& [m, dmDlogT] : massDeriv)
            {
                if (m > sc.minStochMass()) { continue; } // stochastically sampled, handled separately
                const double dmDt = dmDlogT / (age * std::numbers::ln10); // chain rule
                const double weight = std::abs(dmDt) * sc.imf()(m) / sc.imf().expectationValue();
                const std::vector<double> y = f(m, feh);
                for (std::size_t k = 0; k < result.size(); ++k)
                {
                    result[k] += y[k] * weight; // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- f returns n values by its own contract, and result has size n
                }
            }
            return result;
        }

        template <class F>
        [[nodiscard]] auto fehAveragedDeathIntegrand(const double age, const F& f,
            const std::size_t n) const -> std::vector<double>
        {
            const auto& sc = controls_.get();
            const auto& fehDist = sc.fehDist();
            if (fehDist.getMin() == fehDist.getMax()) { return deathIntegrand(age, fehDist.getMin(), f, n); }
            const auto integrandAt = [&](const double feh) -> std::vector<double>
            { return deathIntegrand(age, feh, f, n); };
            auto result = utils::integrateScaled(fehDist, integrandAt, n, fehDist.getMin(),
                fehDist.getMax(), { fehDist.expectationValue() }, sc.intMaxIter(), sc.intAbsTol(),
                sc.intRelTol());
            const double fehNorm = fehDist.integral(fehDist.getMin(), fehDist.getMax());
            for (double& r : result) { r /= fehNorm; }
            return result;
        }

        template <class F, class G>
        [[nodiscard]] auto deathRate(const double t, const F& f, const std::size_t n, const G& g,
            const std::optional<double> fehFixed) const -> std::vector<double>
        {
            const auto& sc = controls_.get();
            const pdfs::PDFReflect sfrAge(sfr(), 0.5 * t);
            const auto integrandAt = [&](const double age) -> std::vector<double>
            {
                return fehFixed.has_value() ? deathIntegrand(age, *fehFixed, f, n) :
                    fehAveragedDeathIntegrand(age, f, n);
            };
            auto result = utils::integrateScaled(sfrAge, integrandAt, n, 0.0, t,
                { t, 0.3 * t, 0.1 * t, 0.03 * t, 0.01 * t }, sc.intMaxIter(), sc.intAbsTol(),
                sc.intRelTol());
            g(t, result);
            return result;
        }

        double curTime_ = 0.0;
        std::vector<Cluster> clusters_;
        std::vector<Cluster> disruptedClusters_;
        std::vector<FieldStar> fieldStars_;
        std::vector<std::optional<specsyn::Specsyn::StarData>> fieldStarProps_;
        std::vector<FieldStar> deadFieldStars_;
        SFRVar sfr_;

        std::vector<double> spec_;
        std::vector<double> specExtinct_;
        std::vector<double> specNeb_;
        std::vector<double> specNebExtinct_;
        std::vector<double> lineLum_;
        std::vector<double> lineLumExtinct_;
        std::vector<double> phot_;
        std::vector<double> photExtinct_;
        std::vector<double> photNeb_;
        std::vector<double> photNebExtinct_;
        double lbol_ = 0.0;
        double targetMass_ = 0.0;
        double actualMass_ = 0.0;
        bool specCurrent_ = true;

        bool photCurrent_ = true;

        bool lbolCurrent_ = true;

        bool windsCurrent_ = true;

        double lbolCts_ = 0.0;

        bool lbolCtsCurrent_ = false;

        std::vector<double> yields_;
        std::vector<double> fieldYields_;
        double lastYieldTime_ = 0.0;

        unsigned long fieldStarStochSNe_ = 0;
        double fieldStarNonStochSNe_ = 0.0;
        double mDotWind_ = 0.0;
        double pDotWind_ = 0.0;
        double eDotWind_ = 0.0;
        double lastFeedbackTime_ = 0.0;

        std::reference_wrapper<const io::SimControls> controls_;

        void computeSpec();

        void addClusterSpecNeb(const extinct::Extinct* ext, const nebular::Nebular* neb);

        void addContinuousSpec(const extinct::Extinct* ext, const nebular::Nebular* neb);

        void addContinuousNebSpec(const extinct::Extinct* ext, const nebular::Nebular* neb,
            const std::vector<double>& contSpec);

        void computePhot();

        void computeLbol();

        void computeLbolCts();

        [[nodiscard]] auto getFieldStarProps() const -> std::vector<std::optional<specsyn::Specsyn::StarData>>;

        void computeYields();

        void computeFeedback();

        void computeWinds();

    };

} // namespace core

#endif // GALAXY_HPP