[WIP] working on new interface to CIB model

src/camel/Astro/cib_interface.cc

+// Ourselves
+#include "cib_interface.h"
+
+// Standard libraries
+#include <iostream>
+#include <cmath>
+
+namespace cib {
+
+  cib_interface::cib_interface()
+  {
+    _initialized_ = false;
+  }
+
+  cib_interface::~cib_interface()
+  {
+  }
+
+  void cib_interface::initialize(// File to be parsed)
+)
+  {
+    // law == "lognormal_sigevol"
+    _parameters_.Meffmax = 8753289339381.791;
+    _parameters_.etamax = 0.4028353504978569;
+    _parameters_.sigmaMh = 1.807080723258688;
+    //_parameters_.alpha = _parameters_.delta = _parameters_.beta = 0.0;
+    _parameters_.tau = 1.2040244128818796;
+
+    _cc_ = {1.076, 1.017, 1.119, 1.097, 1.068, 0.995, 0.960};
+    _fc_ = {1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000};
+    _nu_ = {100e9, 143e9, 217e9, 353e9, 545e9, 857e9, 3000e9}; // Hz
+    _nfreq_ = _nu_.size();
+    _ell_ = {187.,  320.,  502.,  684.,   890., 1158., 1505., 1956};
+    _nell_ = _ell_.size();
+
+    for (double m = 6; m <= 15; m+=0.1) _mh_.push_back(std::pow(10, m));
+
+    cib::read_from_file("data/z.txt", _z_);
+
+    _initialized_ = true;
+  }
+
+  void cib_interface::update(const cib::cib_parameters & parameters_)
+  {
+    if (! _initialized_) {
+      throw std::logic_error("cib_interface::update: CIB model not initialized !");
+    }
+    _parameters_ = parameters_;
+  }
+
+  void cib_interface::compute_cls(cib::array3d & cls_) const
+  {
+    if (! _initialized_) {
+      throw std::logic_error("cib_interface::update: CIB model not initialized !");
+    }
+    std::cout << "cib_interface::compute_cls" << std::endl;
+    cib::zeros(cls_, _nfreq_, _nfreq_, _nell_);
+
+    cib::array3d dj_cen; // [nfreq][mhalo][z]
+    _compute_djc_dlnMh_(dj_cen);
+
+
+    for (size_t l = 0; l < _nell_; l++) {
+      for (size_t f1 = 0; f1 < _nfreq_; f1++) {
+        for (size_t f2 = 0; f2 < _nfreq_; f2++) {
+
+          // CIB
+
+          // tSZ
+
+          // CIBxtSZ
+
+
+        } // end of f2 loop
+      } // end of f1 loop
+    } // end of ell loop
+  }
+
+  void cib_interface::_compute_djc_dlnMh_(cib::array3d & dj_cen_) const
+  {
+    const double fsub = 0.134;
+    cib::array1d mhalo;
+    for (size_t i = 0; i < _mh_.size(); i++) mhalo.push_back(_mh_[i]*(1-fsub));
+    cib::array2d sfr;
+    _compute_sfr_(mhalo, sfr);
+
+    // cib::zeros(dj_cen_, this->_snu_eff_.size(), this->_mh_.size(), this->_z_.size());
+    // for (size_t i = 0; i < this->_snu_eff_.size(); i++) {
+    //   for (size_t j = 0; j < this->_mh_.size(); j++) {
+    //     for (size_t k = 0; k < this->_z_.size(); k++) {
+    //       const double z = this->_z_[k];
+    //       const double cd_Mpc = this->_cosmo_->comoving_distances[k];
+    //       const double rest = this->_hmfmz_[j][k]*sfr[j][k]*(1+z)*std::pow(cd_Mpc, 2)/KC;
+    //       dj_cen_[i][j][k] = rest*this->_snu_eff_[i][k];
+    //     }
+    //   }
+    // }
+  }
+
+
+  void cib_interface::_compute_sfr_(const cib::array1d & mhalo_, cib::array2d & sfr_) const
+  {
+    cib::array2d sfr_mhdot;
+    _compute_sfr_mhdot_(mhalo_, sfr_mhdot);
+    cib::dump(sfr_mhdot);
+    // cib::array2d mhdot;
+    // _compute_Mdot_(mhalo_, mhdot);
+    // cib::array1d ob = this->_cosmo_->Ob;
+    // cib::array1d om = this->_cosmo_->Om;
+    // cib::zeros(sfr_, mhalo_.size(), this->_z_.size());
+    // for (size_t i = 0; i < mhalo_.size(); i++) {
+    //   for (size_t j = 0; j < this->_z_.size(); j++) {
+    //     sfr_[i][j] = mhdot[i][j] * ob[i]/om[i] * sfr_mhdot[i][j];
+    //   }
+    // }
+  }
+
+  void cib_interface::_compute_sfr_mhdot_(const cib::array1d & mhalo_, cib::array2d & sfr_mhdot_) const
+  {
+    const double Meffmax = _parameters_.Meffmax;
+    const double etamax = _parameters_.etamax;
+    const double sigmaMh = _parameters_.sigmaMh;
+    const double tau = _parameters_.tau;
+
+    sfr_mhdot_.resize(mhalo_.size());
+    for (size_t i = 0; i < mhalo_.size(); i++) {
+      if (mhalo_[i] < Meffmax) {
+        const double val = etamax * std::exp(-std::pow(std::log(mhalo_[i]/Meffmax)/sigmaMh, 2)/2);
+        sfr_mhdot_[i].resize(_z_.size(), val);
+      } else {
+        sfr_mhdot_[i].resize(_z_.size());
+        for (size_t j = 0; j < _z_.size(); j++) {
+          const double z_c = 1.5;
+          const double sigpow = sigmaMh - std::max(0., z_c-_z_[j])*tau;
+          const double val = etamax * std::exp(-std::pow(std::log(mhalo_[i]/Meffmax)/sigpow, 2)/2);
+          sfr_mhdot_[i][j] = val;
+        } // end of z loop
+      }
+    } // end of halo mass loop
+  }
+
+  // void cl_cib_tsz_cross::Mdot(const cib::array1d & mhalo_, cib::array2d & mdot_) const
+  // {
+  //   cib::zeros(mdot_, mhalo_.size(), this->_z_.size());
+  //   const bool use_mean = true;
+  //   if (use_mean) {
+  //     for (size_t i = 0; i < mhalo_.size(); i++) {
+  //       const double b = std::pow(mhalo_[i]/1e12, 1.1);
+  //       for (size_t j = 0; j < this->_z_.size(); j++) {
+  //         const double a = 46.1*(1 + 1.11*this->_z_[j]) *       \
+  //           std::sqrt(this->_cosmo_->Om0 * std::pow(1 + this->_z_[j], 3) + this->_cosmo_->Ode0);
+  //         mdot_[i][j] = b*a;
+  //       }
+  //     }
+  //   }
+  // }
+
+
+
+} // end of namespace cib

src/camel/Astro/cib_interface.h

+#ifndef CIB_INTERFACE
+#define CIB_INTERFACE
+
+#include "Astro/utils.h"
+
+namespace cib {
+
+  struct cib_parameters {
+    double Meffmax;
+    double etamax;
+    double sigmaMh;
+    // double alpha;
+    // double beta;
+    // double delta;
+    double tau;
+  };
+
+  class cib_interface {
+
+  public:
+    /// Constructor
+    cib_interface();
+
+    /// Destructor
+    ~cib_interface();
+
+    /// Initialize method
+    void initialize();
+
+    /// Update CIB parameters
+    void update(const cib::cib_parameters & parameters_);
+
+    /// Compute and retrieve Cls
+    void compute_cls(cib::array3d & cls_) const;
+
+  private:
+
+    /// Compute djc_dlnMh
+    void _compute_djc_dlnMh_(cib::array3d & dj_cen_) const;
+
+    void _compute_sfr_(const cib::array1d & mhalo_, cib::array2d & sfr_) const;
+
+    void _compute_sfr_mhdot_(const cib::array1d & mhalo_, cib::array2d & sfr_mhdot_) const;
+
+  private:
+
+    /// Initalized flag
+    bool _initialized_;
+
+    /// CIB parameter model
+    cib_parameters _parameters_;
+
+    /// Frequencies
+    cib::array1d _nu_;
+    size_t _nfreq_;
+
+    /// Multipole
+    cib::array1d _ell_;
+    size_t _nell_;
+
+    /// Color corrections
+    cib::array1d _cc_;
+
+    /// Calibration factors
+    cib::array1d _fc_;
+
+    /// Mass of halo
+    cib::array1d _mh_;
+
+    /// Redshift values
+    cib::array1d _z_;
+  };
+
+} // end of namespace
+#endif