add r_delta

src/camel/Astro/cib_interface.cc

@@ -56,6 +56,7 @@ namespace cib {
     cib::read_from_file("data/u_int.txt", _unfw_);
     cib::read_from_file("data/Pk_int.txt", _power_);
     cib::read_from_file("data/bmz_int.txt", _biasmz_);
+    cib::read_from_file("data/m500c.txt", _m500c_);
 
     _initialized_ = true;
   }
@@ -97,7 +98,6 @@ namespace cib {
     _compute_djsub_dlnMh_(dj_sub);
     cib::array3d Jv;
     _compute_J_nu_(Jv, dj_cen, dj_sub);
-    cib::dump(Jv[0]);
 
     const cib::array3d & u = _unfw_; // [mhalo][ell][z]
     cib::array1d geo;
@@ -111,13 +111,15 @@ namespace cib {
       fcxcc.push_back(_fc_[i]*_cc_[i]);
     }
 
-    // cib::array3d y_ll;
-    // this->y_ell(y_ll);
+    cib::array3d y_ll;
+    _compute_y_ell_(y_ll);
     // cib::array1d g_v;
     // this->g_nu(g_v);
 
 
     cib::zeros(cl_cib_, _nfreq_, _nfreq_, _nell_);
+    cib::zeros(cl_tsz_, _nfreq_, _nfreq_, _nell_);
+    cib::zeros(cl_cibxtsz_, _nfreq_, _nfreq_, _nell_);
     for (size_t l = 0; l < _nell_; l++) {
       for (size_t f1 = 0; f1 < _nfreq_; f1++) {
         for (size_t f2 = 0; f2 < _nfreq_; f2++) {
@@ -153,6 +155,33 @@ namespace cib {
 
 
           // tSZ
+          // for (size_t l = 0; l < this->_ell_.size(); l++) {
+          //   for (size_t f1 = 0; f1 < this->_nfreq_; f1++) {
+          //     for (size_t f2 = 0; f2 < this->_nfreq_; f2++) {
+          //       cib::array1d integral2;
+          //       cib::zeros(integral2, this->_z_.size());
+          //       for (size_t z = 0; z < this->_z_.size(); z++) {
+          //         cib::array1d integral1;
+          //         cib::zeros(integral1, this->_mh_.size());
+          //         for (size_t m = 0; m < this->_mh_.size(); m++) {
+          //           integral1[m] = this->_hmfmz_[m][z]*std::pow(y_l[l][m][z], 2);
+          //         }
+          //         const double dlogm = std::log10(this->_m500c_[1][0]/this->_m500c_[0][0]);
+          //         const double intgn1 = cib::integration_simps(integral1, dlogm);
+          //         integral2[z] = intgn1 * a_z[z];
+          //       }
+          //       const double fin = cib::integration_simps(integral2, &(this->_z_));
+          //       const double T_cmb = 2.725;
+          //       cl_[f1][f2][l] = T_cmb*T_cmb*gnu[f1]*gnu[f2]*fin;
+          //       // {
+          //       //   std::ostringstream oss;
+          //       //   oss << "./data/C_ell_1h_cpp_" << f1 << ".dat";
+          //       //   std::ofstream fout(oss.str());
+          //       //   cib::dump(cl_[f1], fout);
+          //       // }
+          //     }
+          //   }
+          // }
 
           // CIBxtSZ
 
@@ -327,6 +356,56 @@ namespace cib {
     }
   }
 
+  void cib_interface::_compute_y_ell_(cib::array3d & y_ell_) const
+  {
+    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);
+
+    // const double Mpc_to_m = 3.086e22;  // Mpc to m
+    // const double sig_T = 6.6524587158e-29; // m^2 Thomson cross section
+    // const double m_e = 9.10938356e-31; // kg electron mass
+    // const double c_light = 299792458.0; // m/s light speed
+
+    // for (size_t l = 0; l < this->_ell_.size(); l++) {
+    //   for (size_t m = 0; m < this->_mh_.size(); m++) {
+    //     for (size_t z = 0; z < this->_z_.size(); z++) {
+    //       const double cd_Mpc = this->_comoving_distances_[z]/(1e6*CLHEP::parsec);
+    //       const double l500 = cd_Mpc/(1+this->_z_[z])/r500[m][z];
+    //       cib::array1d integral;
+    //       cib::zeros(integral, this->_x_.size());
+    //       for (size_t x = 0; x < this->_x_.size(); x++) {
+    //         const double Pex2 = Pe[m][z][x]*std::pow(this->_x_[x], 2);
+    //         const double x_ls = this->_x_[x]/l500;
+    //         integral[x] =  Pex2 * std::sin(this->_ell_[l]*x_ls)/(this->_ell_[l]*x_ls);
+    //       }
+    //       const double dx = std::log10(this->_x_[1]/this->_x_[0]);
+    //       const double intgn = cib::integration_simps(integral, dx);
+    //       const double a = (sig_T/(m_e*std::pow(c_light, 2)))*(4*CLHEP::pi*r500[m][z]*Mpc_to_m/(l500*l500));
+    //       y_ell_[l][m][z] = intgn*a;
+    //     }
+    //   }
+    // }
+  }
+
+  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++) {
+      for (size_t z = 0; z < _z_.size(); z++) {
+        const double rho_crit = _klass_->critical_density(_z_[z]);
+        const double r3 = 3*_m500c_[m][z]/B/(4*CLHEP::pi*delta_h_tsz*rho_crit);
+        r_delta_[m][z] = std::pow(r3, 1./3.);
+      }
+    }
+  }
+
 
 
 } // end of namespace cib

src/camel/Astro/cib_interface.h

@@ -58,6 +58,10 @@ namespace cib {
 
     void _compute_J_nu_(cib::array3d & Jnu_, const cib::array3d & dj_cen_, const cib::array3d & dj_sub_) const;
 
+    void _compute_r_delta_(cib::array2d & r_delta_) const;
+
+    void _compute_y_ell_(cib::array3d & y_ell_) const;
+
   private:
 
     /// Initalized flag
@@ -104,6 +108,7 @@ namespace cib {
     cib::array2d _snu_eff_;
     cib::array3d _unfw_;
     cib::array2d _biasmz_;
+    cib::array2d _m500c_;
   };
 
 } // end of namespace