Version 0

Compared to cib_tSZ_cross.py

src/camel/Astro/1halo_2halo.cc

@@ -12,37 +12,42 @@
 #include <CLHEP/Units/PhysicalConstants.h>
 
 namespace cib {
-
   
   double cosmology::dchi_dz(double z)
+  //Units in Mpc
   {
-    double dchidz = CLHEP::c_light/(H0*std::sqrt(Om0*std::pow(1+z, 3) + Ode0));
+    double dchidz = (_c_*1e-3)/(cib::cosmology::H0*std::sqrt(cib::cosmology::Om0*std::pow(1+z, 3) + cib::cosmology::Ode0));
     
     return( dchidz);
   }
 
+
   cosmology::cosmology( ClassEngine & klass, 
 			const cib::array1d & z_)
   {
-    H0 = klass.get_H0() * CLHEP::km/CLHEP::s/(1e6*CLHEP::parsec);
+    H0 = klass.get_H0();  //Km.s-1.Mpc-1
     Om0 = klass.Omega_m();
     Ode0 = klass.Omega_L();
+    Tcmb = klass.Tcmb();
     
     for( size_t i=0; i<z_.size(); i++) {
       comoving_distances.push_back( klass.com_distance(z_[i]));
+
+      // const in critical density Msun.Mpc-3
       critical_densities.push_back( klass.critical_density(z_[i]));
+
+      double ez = std::sqrt(cib::cosmology::Om0*std::pow(1+z_[i], 3) + cib::cosmology::Ode0);
+      
       Ob.push_back( klass.get_Ob(z_[i]));
       Om.push_back( klass.get_Om(z_[i]));
       dchidz.push_back( dchi_dz(z_[i]));
-      Ez.push_back( std::sqrt(Om0*std::pow(1+z_[i], 3) + Ode0));
-    }    
+      Ez.push_back( ez);
+      dVcdz.push_back( _c_*std::pow(klass.com_distance(z_[i]),2) / (cib::cosmology::H0 * ez));
+    }
     
   }
   
 
-  const double KC = 1.0e-10; // Kennicutt constant for Chabroer IMF
-  const double fsub = 0.134;
-
   cl_cib_tsz_cross::cl_cib_tsz_cross(const cib::array2d & power_,
                                      const cib::array1d & z_,
                                      const cib::array1d & sig_z_arr_,
@@ -51,7 +56,7 @@ namespace cib {
                                      const cib::array2d & snu_eff_,
                                      const cib::array1d & ell_,
                                      const cib::array3d & unfw_,
-                                     const cib::cosmology & cosmo_,
+                                     ClassEngine & klass,
                                      const cib::array2d & biasmz_,
                                      const double Meffmax_,
                                      const double etamax_,
@@ -70,7 +75,6 @@ namespace cib {
 //     this->_kk_ = kk_;
     this->_power_ = power_;
     this->_z_ = z_;
-//     this->_z_c_ = z_c_;
     this->_sig_z_arr_ = sig_z_arr_;
     this->_mh_ = mh_;
     this->_hmfmz_ = hmfmz_;
@@ -81,10 +85,6 @@ namespace cib {
     this->_Meffmax_ = Meffmax_;
     this->_etamax_ = etamax_;
     this->_sigmaMh_ = sigmaMh_;
-//     this->_alpha_ = alpha_;
-//     this->_delta_ = delta_;
-//     this->_beta_ = beta_;
-//     this->_tau_ = tau_;
     this->_cc_ = cc_;
     this->_fc_ = fc_;
     this->_law_ = law_;
@@ -102,13 +102,12 @@ namespace cib {
     this->_delta_h_tsz_ = delta_h_tsz_;
     this->_B_ = B_;
     this->_x_ = x_;
-    this->_cosmo_ = cib::cosmology(cosmo_);
-    this->_dchidz_ = cosmo_.dchidz;
+    this->_cosmo_ = new cib::cosmology( klass, z_);
+    this->_dchidz_ = this->_cosmo_->dchidz;
     for (size_t i = 0; i < this->_z_.size(); i++) {
-      this->_dchidz_[i] /= std::pow(this->_cosmo_.comoving_distances[i]*(1+this->_z_[i]), 2);
-      this->_dchidz_[i] *= 1e6*CLHEP::parsec;
+      this->_dchidz_[i] /= std::pow(this->_cosmo_->comoving_distances[i]*(1+this->_z_[i]), 2);
     }
-    cib::read_from_file("data/dVcdz.txt", this->_dVcdz_);
+    this->_dVcdz_ = this->_cosmo_->dVcdz;
   }
 
   void cl_cib_tsz_cross::onehalo_int(cib::array3d & Cl_1h_) const
@@ -209,7 +208,7 @@ namespace cib {
             integral2[z] = intgn1 * a_z[z];
           }
           const double fin = cib::integration_simps(integral2, &(this->_z_));
-          const double T_cmb = 2.725;
+	  const double T_cmb = this->_cosmo_->Tcmb;
           cl_[f1][f2][l] = T_cmb*T_cmb*gnu[f1]*gnu[f2]*fin;
           // {
           //   std::ostringstream oss;
@@ -246,7 +245,7 @@ namespace cib {
             integral2[z] = std::pow(intgn1, 2) * a_z[z] * this->_power_[l][z];
           }
           const double fin = cib::integration_simps(integral2, &(this->_z_));
-          const double T_cmb = 2.725;
+          const double T_cmb = this->_cosmo_->Tcmb;
           cl_[f1][f2][l] = T_cmb*T_cmb*gnu[f1]*gnu[f2]*fin;
           // {
           //   std::ostringstream oss;
@@ -262,7 +261,6 @@ namespace cib {
   void cl_cib_tsz_cross::onehalo(cib::array3d & cl_) const
   {
     cib::zeros(cl_, this->_nfreq_, this->_nfreq_, this->_ell_.size());
-    const double T_cmb = 2.725;
     cib::array3d dj_cen; // [nfreq][mhalo][z]
     this->djc_dlnMh(dj_cen);
     cib::array3d dj_sub; // [nfreq][mhalo][z]
@@ -280,7 +278,7 @@ namespace cib {
           cib::array1d intg_mh;
           cib::zeros(intg_mh, this->_z_.size());
           for (size_t z = 0; z < this->_z_.size(); z++) {
-            const double cd_Mpc = this->_cosmo_.comoving_distances[z]/(1e6*CLHEP::parsec);
+            const double cd_Mpc = this->_cosmo_->comoving_distances[z];
             const double cosm = (1+this->_z_[z])*cd_Mpc*cd_Mpc;
             cib::array1d rest1;
             cib::zeros(rest1, this->_mh_.size());
@@ -293,6 +291,7 @@ namespace cib {
             intg_mh[z] *= geo[z];
           }
           const double intg_z = cib::integration_simps(intg_mh, &(this->_z_));
+	  const double T_cmb = this->_cosmo_->Tcmb;
           cl_[f1][f2][l] = T_cmb*intg_z*Kcmb_MJy[f2]*1e6;
           // {
           //   std::ostringstream oss;
@@ -308,7 +307,6 @@ namespace cib {
   void cl_cib_tsz_cross::twohalo(cib::array3d & cl_) const
   {
     cib::zeros(cl_, this->_nfreq_, this->_nfreq_, this->_ell_.size());
-    const double T_cmb = 2.725;
     cib::array3d dj_cen; // [nfreq][mhalo][z]
     this->djc_dlnMh(dj_cen);
     cib::array3d dj_sub; // [nfreq][mhalo][z]
@@ -325,7 +323,7 @@ namespace cib {
           cib::array1d intg_mh;
           cib::zeros(intg_mh, this->_z_.size());
           for (size_t z = 0; z < this->_z_.size(); z++) {
-            const double cd_Mpc = this->_cosmo_.comoving_distances[z]/(1e6*CLHEP::parsec);
+            const double cd_Mpc = this->_cosmo_->comoving_distances[z];
             const double cosm = (1+this->_z_[z])*cd_Mpc*cd_Mpc;
             cib::array1d a1, a2;
             cib::zeros(a1, this->_mh_.size());
@@ -343,6 +341,7 @@ namespace cib {
             intg_mh[z] = geo*intgn_mh1*intgn_mh2;
           }
           const double intg_z = cib::integration_simps(intg_mh, &(this->_z_));
+	  const double T_cmb = this->_cosmo_->Tcmb;
           cl_[f1][f2][l] = T_cmb*intg_z*Kcmb_MJy[f2]*1e6;
           // {
           //   std::ostringstream oss;
@@ -358,9 +357,7 @@ namespace cib {
   void cl_cib_tsz_cross::g_nu(cib::array1d & gnu_) const
   {
     cib::zeros(gnu_, this->_nu_.size());
-    const double h_p = 6.62607004e-34;
-    const double k_B = 1.38064852e-23;
-    const double T_cmb = 2.725;
+    const double T_cmb = this->_cosmo_->Tcmb;
     for (size_t i = 0; i < this->_nu_.size(); i++) {
       const double x = h_p*this->_nu_[i]/k_B/T_cmb;
       gnu_[i] = x*((std::exp(x) + 1)/(std::exp(x) - 1)) - 4;
@@ -376,15 +373,10 @@ namespace cib {
     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->_cosmo_.comoving_distances[z]/(1e6*CLHEP::parsec);
+          const double cd_Mpc = this->_cosmo_->comoving_distances[z];
           const double l500 = cd_Mpc/(1+this->_z_[z])/r500[m][z];
           cib::array1d integral;
           cib::zeros(integral, this->_x_.size());
@@ -407,9 +399,9 @@ namespace cib {
     cib::zeros(c_, this->_mh_.size(), this->_z_.size());
     for (size_t m = 0; m < this->_mh_.size(); m++) {
       for (size_t z = 0; z < this->_z_.size(); z++) {
-        const double h = this->_cosmo_.H0/CLHEP::km*CLHEP::s*(1e6*CLHEP::pc)/100;
+        const double h = this->_cosmo_->H0/100;
         const double M_tilde = this->_m500c_[m][z]*this->_B_;
-        const double a = 1.65*std::pow(h/0.7, 2)*std::pow(this->_cosmo_.Ez[z], 8./3);
+        const double a = 1.65*std::pow(h/0.7, 2)*std::pow(this->_cosmo_->Ez[z], 8./3);
         const double b = std::pow(h/0.7*M_tilde/3e14, 2./3+0.12);
         c_[m][z] = a*b;
       }
@@ -446,12 +438,11 @@ namespace cib {
   void cl_cib_tsz_cross::r_delta(cib::array2d & r_delta_) const
   {
     cib::zeros(r_delta_, this->_mh_.size(), this->_z_.size());
-    cib::array1d rho_crit = this->_cosmo_.critical_densities;
-    const double Msun = 1.98848e+36; // no units to make sure the conversion has the correct unit
+    cib::array1d rho_crit = this->_cosmo_->critical_densities;
     for (size_t m = 0; m < this->_mh_.size(); m++) {
       for (size_t z = 0; z < this->_z_.size(); z++) {
         const double r3 = 3*this->_m500c_[m][z]/(4*CLHEP::pi*this->_delta_h_tsz_*rho_crit[z]);
-        r_delta_[m][z] = std::pow(r3*(Msun/std::pow(1e6*CLHEP::pc, 3)), 1./3.);
+        r_delta_[m][z] = std::pow(r3, 1./3.);  //Mpc
       }
     }
   }
@@ -511,7 +502,7 @@ namespace cib {
           integral[k] = std::min(sfrI[k][j], sfrII[0][j]*ms[k]/mhsub)*subhmf[k]/KC;
         }
         intgn[j] = cib::integration_simps(integral, dlnmsub);
-        const double cd_Mpc = this->_cosmo_.comoving_distances[j]/(1e6*CLHEP::parsec);
+        const double cd_Mpc = this->_cosmo_->comoving_distances[j];
         for (size_t m = 0; m < this->_snu_eff_.size(); m++) {
           dj_sub_[m][i][j] = this->_snu_eff_[m][j]*this->_hmfmz_[i][j]* \
             (1+this->_z_[j])*intgn[j]*std::pow(cd_Mpc, 2);
@@ -533,7 +524,7 @@ namespace cib {
       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]/(1e6*CLHEP::parsec);
+          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];
         }
@@ -547,8 +538,8 @@ namespace cib {
     this->sfr_mhdot(mhalo_, sfr_mhdot);
     cib::array2d mhdot;
     this->Mdot(mhalo_, mhdot);
-    cib::array1d ob = this->_cosmo_.Ob;
-    cib::array1d om = this->_cosmo_.Om;
+    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++) {
@@ -583,7 +574,7 @@ namespace cib {
         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);
+            std::sqrt(this->_cosmo_->Om0 * std::pow(1 + this->_z_[j], 3) + this->_cosmo_->Ode0);
           mdot_[i][j] = b*a;
         }
       }

src/camel/Astro/1halo_2halo.h

@@ -12,6 +12,16 @@
 
 namespace cib {
 
+  const double KC = 1.0e-10; // Kennicutt constant for Chabroer IMF
+  const double fsub = 0.134;
+  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
+  const double h_p = 6.62607004e-34;
+  const double k_B = 1.38064852e-23;
+  const double Msun = 1.98847e30;  //in kg
+
   class cosmology {
   public:
 
@@ -22,13 +32,16 @@ namespace cib {
     double H0;
     double Om0;
     double Ode0;
-    
+
     array1d comoving_distances;
     array1d critical_densities;
     array1d Ob;
     array1d Om;
     array1d Ez;
     array1d dchidz;
+    array1d dVcdz;
+
+    double Tcmb;
     
     //useless ?
 //     const double Oc0;
@@ -48,7 +61,6 @@ namespace cib {
     
     double dchi_dz(double z);
 
-    cosmology();
 
   private:
     // Default constructor can not instantiated
@@ -68,7 +80,7 @@ namespace cib {
 		     const cib::array2d & snu_eff_,
 		     const cib::array1d & ell_,
 		     const cib::array3d & unfw_,
-		     const cib::cosmology & cosmo_,
+		     ClassEngine & klass, 
 		     const cib::array2d & biasmz_,
 		     const double Meffmax_,
 		     const double etamax_,
@@ -139,7 +151,7 @@ namespace cib {
     cib::array2d _snu_eff_;   // [nfreq][z]
     cib::array1d _ell_;       // [ell]
     cib::array3d _unfw_;      // [mhalo][ell][z]
-    cib::cosmology _cosmo_;
+    cib::cosmology *_cosmo_;
     cib::array2d _biasmz_;    // [mhalo][z]
     double _Meffmax_;
     double _etamax_;

src/camel/exec/writeCIB.cc

+#include "Parser.hh"
+#include "Timer.hh"
+#include "Class/MnClassEngine.hh"
+#include"cxxsupport/cxxutils.h"
+#include"cxxsupport/fitshandle.h"
+#include <cmath>
+#include <iostream>
+#include<fstream>
+#include<numeric>
+#include<string>
+using namespace std;
+
+#include <CLHEP/Units/SystemOfUnits.h>
+#include "Astro/1halo_2halo.h"
+
+/*
+Utlities to check the CLASS Pk spectrum: fixed nodes location + inetrpolated values
+
+author : M. Tristram, Apr.2019
+*/
+
+
+static bool abs_compare(double a, double b)
+{
+  return (std::abs(a) < std::abs(b));
+};
+
+
+static void logspace(double kmin,double kmax, int npks,vector<double>&v){
+
+  double lmax=std::log(kmax);
+  double lmin=std::log(kmin);
+  double step=(lmax-lmin)/(npks-1);
+
+  v.resize(npks);
+  for (int i=0;i<npks;i++){
+    double lval=lmin+step*i;
+    v[i]=std::exp(lval);
+  }
+}
+
+
+static void linspace(double kmin,double kmax, int npks,vector<double>&v){
+
+  double step=(kmax-kmin)/(npks-1);
+  v.resize(npks);
+  for (int i=0;i<npks;i++){
+    double kval=kmin+step*i;
+    v[i]=kval;
+  }
+}
+
+
+
+
+
+int main(int argc,char** argv){
+
+  planck_assert(argc==3 || argc==4,"Usage: writeCIB parfile(in) file.fits(out)");
+
+  assert_present(argv[1]);
+
+  float z0 = 10.; /* z for Pk */
+  
+  //decodage arguments
+  Parser parser(argv[1]);
+
+  //choose engine
+  string engine= parser.params.find<string>("engine","CLASS");
+  if (engine=="pico") {
+    throw  Message_error("pico cannot be used for Pk !");
+  }
+
+  ClassParams classparms(parser.classparms);
+  classparms.add("output","mPk"); //polar +lens+clphi
+  
+  classparms.add("z_pk",z0); //polar +lens+clphi  
+
+  //no lensing
+  try{
+    size_t i=classparms.findIndex("lensing");
+    cout << "setting lensing to false" << endl;
+    classparms.updateVal(i,str(false));
+  }
+  catch (Message_error&){
+    classparms.add("lensing",false);
+  }
+  
+  //define kmax and kvalues either from file either given by k_end
+
+    double kmax=-1;
+    vector<double>  kval;
+    string kvfile=parser.params.find<string>("kval_file","");
+    //if file
+    if (kvfile.size()!=0) {
+      assert_present(kvfile);
+      cout <<"using k values from file=" <<kvfile <<endl;
+      //decode 1st column
+      ifstream f(kvfile.c_str());
+      string line;
+      double kk;
+      while (getline(f,line)) {
+	//skip comments
+	if (line.find("#",0) !=std:: string::npos) continue;
+	std::istringstream str(line);
+	if (str>>kk) {
+	  kval.push_back(kk);
+	}
+      }
+      kmax=kval.back();
+    }
+    else{
+      kmax = parser.params.find<double>("kmax",1);
+      const int npks = parser.params.find<int>("n_pks",1000);
+      const double kmin = parser.params.find<double>("kmin",8.e-4);
+      bool lin=parser.params.find<bool>("linspace",false);
+      lin? linspace(kmin,kmax,npks,kval) : logspace(kmin,kmax,npks,kval) ; 
+      
+    }
+
+    //modify class param accordingly
+    bool h_rescale=parser.params.find<bool>("h_rescale",false);