idr.h

#pragma once
#include <vector>
 
#include "background.h"
#include "base_species.h"
#include "perturbations.h"
 
class ThermodynamicsModule;  // forward declaration
 
class IDRSpecies : public BaseSpecies {
 public:
  IDRSpecies(const background& pba,
             double omega0_idr,
             bool has_sibling_idm_dr          = false,
             double T_idr                     = 0.,
             int l_max_idr                    = 0,
             double b_idr                     = 0.,
             int idr_nature                   = 0,
             std::vector<double> alpha_idm_dr = {},
             std::vector<double> beta_idr     = {})
      : BaseSpecies("IDR", EnergyType::Radiation), pba_(pba), Omega0_idr_(omega0_idr),
        has_sibling_idm_dr_(has_sibling_idm_dr), T_idr_(T_idr), l_max_idr_(l_max_idr),
        b_idr_(b_idr), idr_nature_(idr_nature), alpha_idm_dr_(std::move(alpha_idm_dr)),
        beta_idr_(std::move(beta_idr)) {}
 
  bool has_sibling_idm_dr() const {
    return has_sibling_idm_dr_;
  }
 
  double T_idr() const {
    return T_idr_;
  }
  int l_max_idr() const {
    return l_max_idr_;
  }
 
  double b_idr() const {
    return b_idr_;
  }
  int idr_nature() const {
    return idr_nature_;
  }
  const std::vector<double>& alpha_idm_dr() const {
    return alpha_idm_dr_;
  }
  const std::vector<double>& beta_idr() const {
    return beta_idr_;
  }
 
  void SetThermodynamicsModule(const ThermodynamicsModule* thm) override {
    thm_ = thm;
  }
  void SetPerturbs(const perturbs* ppt) override {
    ppt_ = ppt;
  }
 
  double GetOmega0() const override {
    return Omega0_idr_;
  }
 
  double GetRadiationOmega0() const override {
    return Omega0_idr_;
  }
 
  bool IsFreestreaming() const override {
    return true;
  }
 
  void RegisterBackgroundIndices(int& index_bg) override {
    index_bg_rho_ = index_bg++;
  }
 
  void ComputeBackground(double a, const double* /*pvecback_B*/, double* pvecback) override {
    pvecback[index_bg_rho_] = Omega0_idr_ * pba_.H0 * pba_.H0 / (a * a * a * a);
  }
 
  double Rho(const double* pvecback) const override {
    return pvecback[index_bg_rho_];
  }
  double P(const double* pvecback) const override {
    return pvecback[index_bg_rho_] / 3.;
  }
  double PPrime(double a,
                double H,
                const double* /*pvecback_B*/,
                const double* pvecback) const override {
    return a * H * (-4. / 3. * pvecback[index_bg_rho_]);
  }
 
  // ── PerturbLayout ──────────────────────────────────────────────────────────
 
  struct PerturbLayout : BaseSpecies::PerturbLayout {
    int idx_delta = -1;  
    int idx_theta = -1;  
    int idx_shear = -1;  
    int idx_l3    = -1;  
    int l_max     = -1;  
  };
 
  std::unique_ptr<BaseSpecies::PerturbLayout> CreatePerturbLayout() const override {
    return std::make_unique<PerturbLayout>();
  }
 
  double TcaShearIdr(const PerturbLayout& layout,
                     const double* y,
                     const perturb_workspace* ppw) const;
 
  // ── Perturbation index registration ────────────────────────────────────────
 
  void RegisterPerturbationIndices(BaseSpecies::PerturbLayout& layout,
                                   perturb_vector* pv,
                                   const precision* ppr,
                                   int& index_pt,
                                   const perturb_workspace* ppw,
                                   int gauge) override;
 
  // ── PerturbDerivs ──────────────────────────────────────────────────────────
 
  void PerturbDerivs(const BaseSpecies::PerturbLayout& layout,
                     double tau,
                     const double* y,
                     double* dy,
                     const perturb_parameters_and_workspace& ppaw) const override;
 
  // ── Stress-energy observables ───────────────────────────────────────────────
 
  StressEnergyContribution StressEnergy(const BaseSpecies::PerturbLayout& layout,
                                        const perturb_vector* pv,
                                        const double* y,
                                        const double* pvecback,
                                        const perturb_workspace* ppw) const override;
 
  void PerturbSynchronousToNewtonian(const BaseSpecies::PerturbLayout& layout,
                                     double* y,
                                     const PerturbIcContext& ctx) override;
 
  // ── Switch-copy hook (for Phase C / Task 25) ───────────────────────────────
 
  void CopyPerturbationsAcrossSwitch(const BaseSpecies::PerturbLayout& old_layout,
                                     const BaseSpecies::PerturbLayout& new_layout,
                                     const double* old_y,
                                     double* new_y,
                                     const PerturbSwitchContext& ctx) const override;
 
 private:
  const background& pba_;
  double Omega0_idr_;
  bool has_sibling_idm_dr_;
  double T_idr_                    = 0.;
  int l_max_idr_                   = 0;
  const ThermodynamicsModule* thm_ = nullptr;
  const perturbs* ppt_             = nullptr;
  double b_idr_                    = 0.;
  int idr_nature_                  = 0;
  std::vector<double> alpha_idm_dr_;
  std::vector<double> beta_idr_;
};