wdm_decay_product.h

#pragma once
#include <optional>
#include <string>
#include <string_view>
#include <vector>
 
#include "../species/ncdm_base_species.h"
#include "background.h"
 
class BackgroundModule;
class DCDMSpecies;
 
class WdmDecayProductSpecies : public NCDMBaseSpecies {
 public:
  static constexpr std::string_view kTypeName = "dcdm_wdm";
  static constexpr double kQKick              = 10.0;  // p/T0, internal convention
 
  static constexpr double kSparsityFloor = 1e-12;
 
  // ── Fiducial-cosmology placement helpers ───────────────────────────────────
  // Placement-only fixed LCDM background — NOT physics. Maps a decay time t (Mpc)
  // to the scale factor a at which those decays inject, so the momentum grid is a
  // function of Gamma alone (invariant across the run's cosmology). See spec
  // docs/superpowers/specs/2026-07-18-dcdm-wdm-adaptive-momentum-grid-design.md.
  static constexpr double kFidOmegaM = 0.31;
  static constexpr double kFidOmegaR = 9.2e-5;            // photons+nu, T_cmb=2.7255, Neff=3.044
  static constexpr double kFidH      = 0.674;             // h
  static double FiducialCosmicTime(double a);             // t(a) [Mpc], radiation+matter
  static double FiducialScaleFactorAtTime(double t_mpc);  // inverse (monotone)
 
  WdmDecayProductSpecies(FileContent* pfc,
                         const std::string& instance_name,
                         const NcdmSettings& settings,
                         const background* pba,
                         const BackgroundModule* bgm);
 
  double Gamma() const {
    return Gamma_;
  }
  double vkick() const {
    return vkick_;
  }
  double epsilon_retention() const {
    return eps_;
  }
  const std::vector<double>& dq() const {
    return dq_;
  }
  const std::vector<double>& u() const {
    return u_;
  }
  const std::vector<double>& u_edge() const {
    return u_edge_;
  }
  double SigmaAt(double u_cut) const;
 
  void SetParent(const DCDMSpecies* dcdm) {
    parent_ = dcdm;
  }
 
  const std::optional<double>& Omega_ini_pending() const {
    return Omega_ini_pending_;
  }
  const std::optional<double>& Omega_combined_pending() const {
    return Omega_combined_pending_;
  }
 
  void FillInjection(double a, double rho_dcdm, double* J, double* dJdlnq) const;
 
  // ── Background (Task 2) ────────────────────────────────────────────────────
  void RegisterBackgroundIndices(int& index_bg) override;
  void RegisterIntegrationIndices(int& index_bi) override;
  void SetBackgroundInitialConditions(const BackgroundICContext& ctx) override;
  void ComputeBackground(double a, const double* pvecback_B, double* pvecback) override;
 
  double Rho(const double* pvecback) const override {
    return pvecback[index_bg_rho_];
  }
  double P(const double* pvecback) const override {
    return pvecback[index_bg_p_];
  }
  double PPrime(double a,
                double H,
                const double* pvecback_B,
                const double* pvecback) const override;
  double RhoDotOverRho(const double* pvecback, double a_prime_over_a) const override;
 
  void WriteBackgroundColumnTitles(BackgroundColumnWriter& w) const override;
  void WriteBackgroundData(const double* pvecback, BackgroundColumnWriter& w) const override;
 
  // ── Decay-product neutrality overrides ─────────────────────────────────────
  // The daughter starts empty: it contributes nothing to the early-time
  // radiation budget, fnu, N_eff, or the flatness budget (the combined sector
  // density is owned by the DCDM_WDM composite).
  double GetOmega0() const override {
    return 0.;
  }
  double NeutrinoOmega0() const override {
    return 0.;
  }
  double NeffContribution(double /*z*/) const override {
    return 0.;
  }
  bool IsFreestreaming() const override {
    return false;
  }
  double BackgroundAIni(double a_proposed, double /*tol*/) const override {
    return a_proposed;  // empty at early times; no relativistic-start constraint
  }
  void CheckUltraRelativisticAtIc(const double* /*pvecback*/, double /*tol*/) const override {}
  bool IsUltraRelativisticAtIc(const double* /*pvecback*/, double /*tol*/) const override {
    return true;
  }
  void PrintNeffInfo() const override;
  void PrintMassInfo() const override;
 
  // ── Perturbations (Task 4) ─────────────────────────────────────────────────
  void RegisterPerturbationIndices(BaseSpecies::PerturbLayout& layout,
                                   perturb_vector* pv,
                                   const precision* ppr,
                                   int& index_pt,
                                   const perturb_workspace* ppw,
                                   int gauge) override;
  void PerturbDerivs(const BaseSpecies::PerturbLayout& layout,
                     double tau,
                     const double* y,
                     double* dy,
                     const perturb_parameters_and_workspace& ppaw) const override;
  void ApplyInitialConditions(const BaseSpecies::PerturbLayout& layout,
                              double* y,
                              const PerturbIcContext& ctx) override;
  StressEnergyContribution StressEnergy(const BaseSpecies::PerturbLayout& layout,
                                        const perturb_vector* pv,
                                        const double* y,
                                        const double* pvecback,
                                        const perturb_workspace* ppw) const override;
  void FillSources(const BaseSpecies::PerturbLayout& layout,
                   const double* y,
                   const double* dy,
                   PerturbSourceContext& ctx) const override;
  void WriteOutputColumns(
      PerturbColumnWriter& writer,
      const PerturbationsModule& mod,
      file_format fmt,
      TransferColumnSection section = TransferColumnSection::all) const override;
 
  void RegisterTensorPerturbationIndices(BaseSpecies::PerturbLayout& /*layout*/,
                                         perturb_vector* /*pv*/,
                                         const precision* /*ppr*/,
                                         int& /*index_pt*/,
                                         const perturb_workspace* /*ppw*/,
                                         int /*gauge*/) override {}
  void PerturbSynchronousToNewtonian(const BaseSpecies::PerturbLayout& /*layout*/,
                                     double* /*y*/,
                                     const PerturbIcContext& /*ctx*/) override {}
 
  // ── Index accessors for the composite ──────────────────────────────────────
  int bg_f_index() const {
    return index_bg_f_;
  }
  int bg_dfdlnq_index() const {
    return index_bg_dfdlnq_;
  }
  int bg_inj_index() const {
    return index_bg_inj_;
  }
  int bg_number_index() const {
    return index_bg_number_;
  }
  int bi_f_index() const {
    return index_bi_f_;
  }
  int bi_dfdlnq_index() const {
    return index_bi_dfdlnq_;
  }
 
 protected:
  double GetDlnf0Dlnq(int /*iq*/, const double* /*pvecback*/) const override {
    return 0.;
  }
 
 private:
  static constexpr double kFSeed = 1e-10;
 
  void BuildInjectionAdaptedGrid();
  void GaussWeights(double u_cut, double sigma, double* w) const;
 
  const background* pba_;
  const DCDMSpecies* parent_ = nullptr;
 
  double Gamma_ = 0.;  // decay rate, CLASS units (Mpc^-1 after both conversions)
  double vkick_ = 0.;  // kick velocity in units of c
  double eps_   = 0.;  // total mass retention (m_d = eps*m_parent/2)
 
  int n_bins_        = 96;
  double q_edge_tol_ = 1e-3;  // decayed-fraction trimmed off each end (grid F-span);
                              // the trimmed tail is clamped into the edge bin, not dropped
  int l_max_input_ = -1;      // -1 → ppr->l_max_ncdm
  std::optional<double> q_min_ratio_floor_;  // optional hard floor on q_lo/q_kick
 
  static constexpr double kSigmaMin = 0.03;  // ≈ 4 background-table samples (7e-3 in ln a)
  static constexpr double kSigmaMax = 0.25;  // bounds the wide-bin placement bias
  double F_lo_ = 0.;  // realized decayed-fraction span of the grid (0,0 = fallback grid)
  double F_hi_ = 0.;
 
  std::vector<double> u_;       // ln q_i (bin centres = q(F_mid[i]))
  std::vector<double> dq_;      // cell widths q_edge[i+1] - q_edge[i]
  std::vector<double> u_edge_;  // N+1 cell edges in u (F-quantile edges)
 
  std::optional<double> Omega_ini_pending_;
  std::optional<double> Omega_combined_pending_;
 
  // Background indices (per-bin arrays are n_bins_ long)
  int index_bg_number_   = -1;
  int index_bg_pseudo_p_ = -1;
  int index_bg_f_        = -1;
  int index_bg_dfdlnq_   = -1;
  int index_bg_inj_      = -1;
  int index_bi_f_        = -1;
  int index_bi_dfdlnq_   = -1;
 
  mutable std::vector<double> scratch_w_;  // erf energy-deposit weights (background thread only)
};