CIB Cl done

src/camel/Astro/cib_interface.cc

@@ -54,6 +54,8 @@ namespace cib {
     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_);
+    cib::read_from_file("data/Pk_int.txt", _power_);
+    cib::read_from_file("data/bmz_int.txt", _biasmz_);
 
     _initialized_ = true;
   }
@@ -93,6 +95,9 @@ namespace cib {
     _compute_djc_dlnMh_(dj_cen);
     cib::array3d dj_sub; // [nfreq][mhalo][z]
     _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,6 +116,7 @@ namespace cib {
     // cib::array1d g_v;
     // this->g_nu(g_v);
 
+
     cib::zeros(cl_cib_, _nfreq_, _nfreq_, _nell_);
     for (size_t l = 0; l < _nell_; l++) {
       for (size_t f1 = 0; f1 < _nfreq_; f1++) {
@@ -119,7 +125,10 @@ namespace cib {
           // CIB
           cib::array1d intg_mh;
           cib::zeros(intg_mh, _z_.size());
+          cib::array1d rest2;
+          cib::zeros(rest2, _z_.size());
           for (size_t z = 0; z < _z_.size(); z++) {
+            // One halo
             cib::array1d rest1;
             cib::zeros(rest1, _mh_.size());
             for (size_t m = 0; m < _mh_.size(); m++) {
@@ -130,9 +139,17 @@ namespace cib {
             const double dm = std::log10(_mh_[1]/_mh_[0]);
             intg_mh[z] = cib::integration_simps(rest1, dm);
             intg_mh[z] *= geo[z];
+            // Two halos
+            rest2[z] = Jv[f2][z][l]*Jv[f1][z][l]*geo[z]*_power_[l][z];
+          }
+          {
+            const double intg_z = cib::integration_simps(intg_mh, &(_z_));
+            cl_cib_[f1][f2][l] += fcxcc[f1]*intg_z*fcxcc[f2];
+          }
+          {
+            const double intg_z = cib::integration_simps(rest2, &(_z_));
+            cl_cib_[f1][f2][l] += fcxcc[f1]*intg_z*fcxcc[f2];
           }
-          const double intg_z = cib::integration_simps(intg_mh, &(_z_));
-          cl_cib_[f1][f2][l] = fcxcc[f1]*intg_z*fcxcc[f2];
 
 
           // tSZ
@@ -265,6 +282,29 @@ namespace cib {
     }
   }
 
+  void cib_interface::_compute_J_nu_(cib::array3d & Jnu_,
+                                     const cib::array3d & dj_cen_,
+                                     const cib::array3d & dj_sub_) const
+  {
+    cib::zeros(Jnu_, _nfreq_, _z_.size(), _ell_.size());
+    const cib::array3d & u = _unfw_; // [mhalo][ell][z]
+    for (size_t l = 0; l < _ell_.size(); l++) {
+      for (size_t f = 0; f < _nfreq_; f++) {
+        cib::array1d intg_mh;
+        cib::zeros(intg_mh, _z_.size());
+        for (size_t z = 0; z < _z_.size(); z++) {
+          cib::array1d rest1;
+          cib::zeros(rest1, _mh_.size());
+          for (size_t m = 0; m < _mh_.size(); m++) {
+            rest1[m] = (dj_cen_[f][m][z] + dj_sub_[f][m][z]*u[m][l][z])*_biasmz_[m][z];
+          }
+          const double dm = std::log10(_mh_[1]/_mh_[0]);
+          Jnu_[f][z][l] = cib::integration_simps(rest1, dm);
+        }
+      }
+    }
+  }
+
   void cib_interface::_compute_msub_(const double mhalo_, cib::array1d & msub_) const
   {
     const double log10msub_min = 5;

src/camel/Astro/cib_interface.h

@@ -56,6 +56,7 @@ namespace cib {
 
     void _compute_subhmf_(const double mhalo_, const cib::array1d & ms_, cib::array1d & subhmf_) const;
 
+    void _compute_J_nu_(cib::array3d & Jnu_, const cib::array3d & dj_cen_, const cib::array3d & dj_sub_) const;
 
   private:
 
@@ -76,6 +77,9 @@ namespace cib {
     cib::array1d _ell_;
     size_t _nell_;
 
+    // Pk
+    cib::array2d _power_;
+
     /// Color corrections
     cib::array1d _cc_;
 
@@ -99,6 +103,7 @@ namespace cib {
     cib::array2d _hmfmz_;
     cib::array2d _snu_eff_;
     cib::array3d _unfw_;
+    cib::array2d _biasmz_;
   };
 
 } // end of namespace