add P_e

src/camel/Astro/cib_interface.cc

@@ -48,9 +48,15 @@ namespace cib {
     // Redshift & comoving distances
     cib::read_from_file("data/z.txt", _z_);
     _comoving_distances_.resize(_z_.size());
+    _Ez_.resize(_z_.size());
     _dVcdz_.resize(_z_.size());
     _dchidz_.resize(_z_.size());
 
+    const size_t xstep = 50;
+    const double xmin = -6, xmax = 1;
+    const double dx = (xmax - xmin)/(xstep-1);
+    for (size_t i = 0; i < xstep; i++) _x_.push_back(std::pow(10, xmin + i*dx));
+
     cib::read_from_file("data/hmf_red.txt", _hmfmz_);
     cib::read_from_file("data/snu.txt", _snu_eff_);
     cib::read_from_file("data/u_int.txt", _unfw_);
@@ -78,6 +84,7 @@ namespace cib {
       _comoving_distances_[i] = _klass_->com_distance(z);
 
       const double ez = std::sqrt(Om0*std::pow(1+z, 3) + Ode0);
+      _Ez_[i] = ez;
       _dVcdz_[i] = _c_*std::pow(_comoving_distances_[i], 2) / (H0 * ez);
       _dchidz_[i] = _c_*1e-3/ (H0 * ez);
     }
@@ -361,9 +368,8 @@ namespace cib {
     cib::zeros(y_ell_, _ell_.size(), _mh_.size(), _z_.size());
     cib::array2d r500;
     _compute_r_delta_(r500);
-    cib::dump(r500);
-    // cib::array3d Pe;
-    // this->P_e(Pe);
+    cib::array3d Pe;
+    _compute_P_e_(Pe);
 
     // const double Mpc_to_m = 3.086e22;  // Mpc to m
     // const double sig_T = 6.6524587158e-29; // m^2 Thomson cross section
@@ -394,7 +400,6 @@ namespace cib {
   void cib_interface::_compute_r_delta_(cib::array2d & r_delta_) const
   {
     cib::zeros(r_delta_, _mh_.size(), _z_.size());
-    const double Msun = 1.98848e+36; // no units to make sure the conversion has the correct unit
     const double delta_h_tsz = 500;
     const double B = 1.28;
     for (size_t m = 0; m < _mh_.size(); m++) {
@@ -406,6 +411,38 @@ namespace cib {
     }
   }
 
+  void cib_interface::_compute_P_e_(cib::array3d & pe_) const
+  {
+    // values of the constants taken from
+    // https://www.aanda.org/articles/aa/pdf/2013/02/aa20040-12.pdf
+    cib::zeros(pe_, _mh_.size(), _z_.size(), _x_.size());
+    const double gamma_t = 0.31;
+    const double alpha_t = 1.33;
+    const double beta_t  = 4.13;
+    const double P_0_t   = 6.41;
+    const double c_500_t = 1.81;
+
+    const double eV_to_J = 1.6e-19;
+    const double cm_to_m = 0.01;
+    const double B = 1.28;
+    for (size_t m = 0; m < _mh_.size(); m++) {
+      for (size_t z = 0; z < _z_.size(); z++) {
+        const double h = _klass_->get_H0()/100;
+        // const double h = p15.H0/CLHEP::km*CLHEP::s*(1e6*CLHEP::pc)/100;
+        const double M_tilde = _m500c_[m][z]*B;
+        const double a1 = 1.65*std::pow(h/0.7, 2)*std::pow(_Ez_[z], 8./3);
+        const double b1 = std::pow(h/0.7*M_tilde/3e14, 2./3+0.12);
+        const double c1 = a1*b1;
+        for (size_t x = 0; x < this->_x_.size(); x++) {
+          const double c2 = c1*eV_to_J/std::pow(cm_to_m, 3);
+          const double a2 = c2*P_0_t*std::pow(c_500_t*_x_[x], -gamma_t);
+          const double b2 = std::pow(1 + std::pow(c_500_t*_x_[x], alpha_t), (gamma_t-beta_t)/alpha_t);
+          pe_[m][z][x] = a2*b2;
+        }
+      }
+    }
+  }
+
 
 
 } // end of namespace cib

src/camel/Astro/cib_interface.h

@@ -62,6 +62,8 @@ namespace cib {
 
     void _compute_y_ell_(cib::array3d & y_ell_) const;
 
+    void _compute_P_e_(cib::array3d & pe_) const;
+
   private:
 
     /// Initalized flag
@@ -102,6 +104,7 @@ namespace cib {
     ///
     cib::array1d _dVcdz_;
     cib::array1d _dchidz_;
+    cib::array1d _Ez_;
 
     /// ?
     cib::array2d _hmfmz_;
@@ -109,6 +112,7 @@ namespace cib {
     cib::array3d _unfw_;
     cib::array2d _biasmz_;
     cib::array2d _m500c_;
+    cib::array1d _x_;
   };
 
 } // end of namespace