From fc20f5aa924452ce3feda2cea771ec5522820a55 Mon Sep 17 00:00:00 2001 From: Matthieu Tristram Date: Mon, 29 Apr 2019 17:57:47 +0200 Subject: [PATCH] Version 0 Compared to cib_tSZ_cross.py --- src/camel/Astro/1halo_2halo.cc | 79 +++---- src/camel/Astro/1halo_2halo.h | 20 +- src/camel/exec/writeCIB.cc | 412 +++++++++++++++++++++++++++++++++ 3 files changed, 463 insertions(+), 48 deletions(-) create mode 100644 src/camel/exec/writeCIB.cc diff --git a/src/camel/Astro/1halo_2halo.cc b/src/camel/Astro/1halo_2halo.cc index 80fe690..a2472a7 100644 --- a/src/camel/Astro/1halo_2halo.cc +++ b/src/camel/Astro/1halo_2halo.cc @@ -12,37 +12,42 @@ #include 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_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; } } diff --git a/src/camel/Astro/1halo_2halo.h b/src/camel/Astro/1halo_2halo.h index 9433e35..ca22512 100644 --- a/src/camel/Astro/1halo_2halo.h +++ b/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_; diff --git a/src/camel/exec/writeCIB.cc b/src/camel/exec/writeCIB.cc new file mode 100644 index 0000000..bd6a884 --- /dev/null +++ b/src/camel/exec/writeCIB.cc @@ -0,0 +1,412 @@ +#include "Parser.hh" +#include "Timer.hh" +#include "Class/MnClassEngine.hh" +#include"cxxsupport/cxxutils.h" +#include"cxxsupport/fitshandle.h" +#include +#include +#include +#include +#include +using namespace std; + +#include +#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&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&v){ + + double step=(kmax-kmin)/(npks-1); + v.resize(npks); + for (int i=0;i("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 kval; + string kvfile=parser.params.find("kval_file",""); + //if file + if (kvfile.size()!=0) { + assert_present(kvfile); + cout <<"using k values from file=" <>kk) { + kval.push_back(kk); + } + } + kmax=kval.back(); + } + else{ + kmax = parser.params.find("kmax",1); + const int npks = parser.params.find("n_pks",1000); + const double kmin = parser.params.find("kmin",8.e-4); + bool lin=parser.params.find("linspace",false); + lin? linspace(kmin,kmax,npks,kval) : logspace(kmin,kmax,npks,kval) ; + + } + + //modify class param accordingly + bool h_rescale=parser.params.find("h_rescale",false); + string label=(h_rescale? "P_k_max_h/Mpc" : "P_k_max_1/Mpc"); + cout << label << " set to " << kmax <("precisionFile","")); + } + + + cout << "CLASS \t--> precision parameters taken from file=" << pre << endl; + MnClassEngine* e=new MnClassEngine(parser.vars(), + classparms, + pre); + + ClassEngine* klass=e->getClass(); + + //check cosmology + std::cout << "Omega_c = " << klass->Omega_c() << std::endl; + std::cout << "Omega_b = " << klass->Omega_b() << std::endl; + std::cout << "Omega_m = " << klass->Omega_m() << std::endl; + std::cout << "Omega_L = " << klass->Omega_L() << std::endl; + std::cout << "Omega_g = " << klass->Omega_g() << std::endl; + std::cout << "Omega_nu = " << klass->Omega_n() << std::endl; + std::cout << " H0 = " << (klass->get_H0() * CLHEP::km/CLHEP::s/(1e6*CLHEP::parsec)) /CLHEP::km*CLHEP::s << " km/s" << std::endl; + std::cout << " ns = " << klass->ns() << std::endl; + std::cout << " sigma8 = " << klass->get_sigma8(0) << std::endl; + std::cout << " tau = " << klass->get_tau_reio() << std::endl; + std::cout << " z_reio = " << klass->get_z_reio() << std::endl; + std::cout << " Tcmb = " << klass->Tcmb() << std::endl; + std::cout << " Neff = " << klass->Neff() << std::endl; + const double kB_evK = 8.617330337217213e-05; + std::cout << " Mnu = "; for( int i=0; iMnu().size(); i++) std::cout << klass->Mnu()[i]*kB_evK << ", "; std::cout << std::endl; + + + double _hubble_distance = CLHEP::c_light/(klass->get_H0() * CLHEP::km/CLHEP::s/(1e6*CLHEP::parsec)); + std::cout << " hubble_distance = " << _hubble_distance/(1e6*CLHEP::parsec) << " Mpc" << std::endl; + + vector zvalues; + for (double z = 0; z <= 10; z+=0.5) { + zvalues.push_back(z); + std::cout << "Comoving distance for z = " << z << ": " << klass->com_distance(z) << " Mpc" << std::endl; + } + + cout << "Omega_m=" << klass->Omega_m() << "\tomega_m=" << klass->Omega_m()*klass->get_h()*klass->get_h() << endl; + cout << "sigma8(0)=" << klass->get_sigma8(0) << "\tsigma8(" << z0 << ")=" << klass->get_sigma8(z0) << endl; + + cout << "Ob(z)" << endl; + for(size_t iz = 0; iz < zvalues.size(); iz++) { + cout << klass->get_Ob(zvalues[iz]) << " "; + } + cout << endl; + + cout << "Om(z)" << endl; + for(size_t iz = 0; iz < zvalues.size(); iz++) { + cout << klass->get_Om(zvalues[iz]) << " "; + } + cout << endl; + + cout << "H(z)/H0" << endl; + for(size_t iz = 0; iz < zvalues.size(); iz++) { + cout << klass->get_H(zvalues[iz]) / klass->get_H0() << " "; + } + cout << endl; + + cout << "H(z)" << endl; //Units in + for(size_t iz = 0; iz < zvalues.size(); iz++) { + cout << klass->get_H(zvalues[iz]) << " "; + } + cout << endl; + + cout << "Critical density" << endl; //Units in Msun/Mpc3 + for (size_t iz = 0; iz < zvalues.size(); iz++) { + cout << klass->critical_density(zvalues[iz]) << " "; + } + cout << endl; + + + + + + + cout << "Computing P(k) at z=" << z0 << endl; + + //test fast Pk + Timer timer; + const size_t nloop=100; + cout <<"Pk timer test..." << endl; + vector pkslow(kval.size()); + for (size_t i=0;iget_Pk(kval[k], z0); + } + } + cout << "std=" << timer.partial() < pkfast(kval.size()); + for (size_t i=0;iget_Pkvec(kval,z0,pkfast); + cout << "fast=" << timer.partial() < residual(kval.size()); + std::transform(pkfast.begin(),pkfast.end(),pkslow.begin(), + residual.begin(),std::minus()); + + double diff=*max_element(residual.begin(),residual.end(),abs_compare); + cout << "max diff=" << diff << endl; + + //fits output + const std::string fileName(argv[2]); + cout << "writing FITS file =" << fileName << endl; + if (file_present(fileName)) { + //std::cout << "the file " + fileName + " already exists: removing"< pkcols; + pkcols.push_back(fitscolumn("knode", "1/Mpc",1,TDOUBLE)); + pkcols.push_back(fitscolumn("pklin", "(Mpc)^3" ,1,TDOUBLE)); + pkcols.push_back(fitscolumn("pklin2", "(Mpc)^3" ,1,TDOUBLE)); + fout.insert_bintab(pkcols); + + vector knode, pknode, pknode2; + + klass->get_PklinNodes(knode,pknode,0); + klass->get_PklinNodes(knode,pknode2,1); + + for (size_t i=0;i pkcols; + pkcols.push_back(fitscolumn("k", "1/Mpc",1,TDOUBLE)); + pkcols.push_back(fitscolumn("pklin", "(Mpc)^3" ,1,TDOUBLE)); + + fout.insert_bintab(pkcols); + + //optional rescale + double h=klass->get_H0()/100.; + //warning h_rescaling + if (h_rescale) + std::transform(kval.begin(),kval.end(),kval.begin(),bind1st(multiplies(),h)); + cout <<"computing Pk from " << kval[0] << " to " << kval.back() << endl; + + for (size_t i=0;iget_Pklin(kval[i],z0),i); + } + cout <<"Growth factor D(z=" << z0 << ")="<< klass->get_growthD(z0) << endl; + fout.add_key("D",klass->get_growthD(z0),"growth factor"); + fout.add_key("f",klass->get_f(z0),"growth rate"); + + } + + fout.close(); + + + //CIB model + cib::array2d k_array; + cib::read_from_file("data/k_array.txt", k_array); + cib::dump_shape(k_array); + cib::array2d Pk_int; + cib::read_from_file("data/Pk_int.txt", Pk_int); + cib::dump_shape(Pk_int); + cib::array1d z; + cib::read_from_file("data/z.txt", z); + cib::dump_shape(z); + + const double z_c = 1.5; + cib::array1d sig_z_arr; + for (size_t i = 0; i < z.size(); i++) sig_z_arr.push_back(std::max(0., z_c-z[i])); + + cib::array1d mass; + for (double m = 6; m <= 15; m+=0.1) mass.push_back(std::pow(10, m)); + cib::dump_shape(mass); + + cib::array2d hmf_red; + cib::read_from_file("data/hmf_red.txt", hmf_red); + cib::dump_shape(hmf_red); + + cib::array2d snu; + cib::read_from_file("data/snu.txt", snu); + cib::dump_shape(snu); + + double dell[8] = {187., 320., 502., 684., 890., 1158., 1505., 1956}; + cib::array1d ell( &dell[0], &dell[0]+8); + cib::dump(ell); + + cib::array3d u_int; + cib::read_from_file("data/u_int.txt", u_int); + cib::dump_shape(u_int); + // cib::dump(u_int[0][5]); + + cib::array2d bmz_int; + cib::read_from_file("data/bmz_int.txt", bmz_int); + + //parameter + const std::string law = "lognormal_sigevol"; + double Meffmax = 8753289339381.791; + double etamax = 0.4028353504978569; + double sigmaMh = 1.807080723258688; + double tau = 1.2040244128818796; + + double dcc[7] = {1.076, 1.017, 1.119, 1.097, 1.068, 0.995, 0.96}; + cib::array1d cc_pl(&dcc[0],&dcc[0]+7); + cib::array1d fc_pl(7, 1); + + double dnu[7] = {100e9, 143e9, 217e9, 353e9, 545e9, 857e9, 3000e9}; // GHz + cib::array1d nu(&dnu[0],&dnu[0]+7); + + cib::array2d m500c; + cib::read_from_file("data/m500c.txt", m500c); + + const double delta_h_tsz = 500; + const double B = 1.28; + cib::array1d x; + 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::cosmology *cosmo = new cib::cosmology( *klass, z); +// cout << "dchi_dz" << endl; +// for (size_t iz = 0; iz < zvalues.size(); iz++) { +// cout << cosmo->dchidz[iz] << " "; +// } +// cout << endl; +// cout << "dVc_dz" << endl; +// for (size_t iz = 0; iz < zvalues.size(); iz++) { +// cout << cosmo->dVcdz[iz] << " "; +// } +// cout << endl; +// cout << "rho_crit" << endl; +// for (size_t iz = 0; iz < z.size(); iz++) { +// cout << cosmo->critical_densities[iz] << " "; +// } +// cout << endl; + + + cib::cl_cib_tsz_cross cl_cib(Pk_int, + z, + sig_z_arr, + mass, + hmf_red, + snu, + ell, + u_int, + *klass, + bmz_int, + Meffmax, + etamax, + sigmaMh, + tau, + cc_pl, fc_pl, + law, + nu, + m500c, + delta_h_tsz, + B, + x); + + cib::array3d cl; + cl_cib.twohalo(cl); + cib::dump_shape(cl); + cib::dump(cl[0][0]); + + delete e; + + return 0; +} -- GitLab