perturb_source_context.h

#pragma once
#include <string>
#include <vector>
 
#include "common.h"  // for file_format, true, false, class_store_* macros
 
// Forward declarations
struct perturb_workspace;
struct perturb_vector;
struct precision;
class PerturbationsModule;
 
// ── existing PerturbScalarContext (unchanged) ────────────────────────────────
 
struct PerturbScalarContext {
  double k = 0., k2 = 0.;
  double a = 0., a2 = 0., a_prime_over_a = 0.;
  double metric_continuity = 0.;
  double metric_euler      = 0.;
  double metric_shear      = 0.;
  double metric_ufa_class  = 0.;
  double cotKgen = 0., s2_squared = 1.;
  double delta_g = 0., theta_g = 0.;
  double delta_b = 0., theta_b = 0.;
  int b_idx_theta = -1;
  int idr_nature = 0;
  double R = 0.;
  double cb2 = 0.;
  int gauge = 0;
};
 
struct SourceRequestContext {
  bool wants_density;   // ppt->has_density_transfers
  bool wants_velocity;  // ppt->has_velocity_transfers
  int gauge;            // possible_gauges; CDM registers theta only when != synchronous
};
 
// ── PerturbColumnWriter ──────────────────────────────────────────────────────
class PerturbColumnWriter {
 public:
  // Title mode (used from perturb_output_titles and perturb_prepare_k_output)
  explicit PerturbColumnWriter(std::string& titles) : titles_(&titles) {}
 
  // Data mode for WriteOutputColumns: value = tk[tp_index]
  PerturbColumnWriter(double* dataptr, const double* tk, int& storeidx)
      : dataptr_(dataptr), tk_(tk), storeidx_(&storeidx) {}
 
  // Data mode for PrintVariables: value passed directly
  PerturbColumnWriter(double* dataptr, int& storeidx) : dataptr_(dataptr), storeidx_(&storeidx) {}
 
  bool IsTitleMode() const {
    return titles_ != nullptr;
  }
 
  // For WriteOutputColumns: species pass a tp_index; writer looks up tk[tp_index]
  void Add(const char* title, int tp_index, bool active);
  void Add(const std::string& title, int tp_index, bool active) {
    Add(title.c_str(), tp_index, active);
  }
 
  // For PrintVariables: species pass the value directly
  void Add(const char* title, double value, bool active);
  void Add(const std::string& title, double value, bool active) {
    Add(title.c_str(), value, active);
  }
 
 private:
  std::string* titles_ = nullptr;
  double* dataptr_     = nullptr;
  const double* tk_    = nullptr;
  int* storeidx_       = nullptr;
};
 
// ── PerturbSourceContext ─────────────────────────────────────────────────────
struct PerturbSourceContext {
  PerturbationsModule* p_mod = nullptr;
  perturb_workspace* ppw     = nullptr;
  int index_md               = 0;
  int index_ic               = 0;
  int index_k                = 0;
  int index_tau              = 0;
  double k                   = 0.;
  double a                   = 0.;
  double a2                  = 0.;
  double a_prime_over_a      = 0.;
  // Pre-computed N-body gauge corrections (0 when not has_Nbody_gauge_transfers)
  double theta_over_k2 = 0.;
  double theta_shift   = 0.;
};
 
// ── PerturbIcContext ─────────────────────────────────────────────────────────
struct PerturbIcContext {
  // Density fractions
  double fracnu = 0., fracg = 0., fracb = 0., fraccdm = 0., fracidm_drmd = 0.;
  double rho_m_over_rho_r = 0., om = 0.;
  // (k*tau)^2, (k*tau)^3
  double ktau_two = 0., ktau_three = 0.;
  // Curvature factor: 1 - 3K/k^2
  double s2_squared = 1.;
  // IC seed values pre-computed by module (photon delta/theta set by module formulas)
  double delta_g_ic = 0., theta_g_ic = 0.;
  // Shared relativistic IC (= delta_g_ic for adiabatic, computed by module for isocurvatures)
  double delta_ur = 0., theta_ur = 0., shear_ur = 0., l3_ur = 0., delta_dr = 0.;
  // Synchronous metric IC
  double eta = 0.;
  // DO NOT read in ApplyInitialConditions() — filled by module AFTER the species loop
  // for the Newtonian gauge transformation. Zero during the species dispatch.
  double alpha = 0., alpha_prime = 0.;
  // Kinematics
  double k = 0., tau = 0., a = 0., a_prime_over_a = 0.;
  int index_ic = 0;
  int gauge    = 0;
  // Pointers
  perturb_workspace* ppw           = nullptr;
  const precision* ppr             = nullptr;
  const PerturbationsModule* p_mod = nullptr;
};
 
// ── PerturbSwitchContext ─────────────────────────────────────────────────────
struct PerturbSwitchContext {
  double k               = 0.;
  double a               = 0.;
  const double* pvecback = nullptr;
};