Commit 9041755c authored by Xavier Garrido's avatar Xavier Garrido
Browse files

done tSZ and tSZxCIB

parent b8f65e4d
......@@ -42,6 +42,15 @@ namespace cib {
_ell_ = {187., 320., 502., 684., 890., 1158., 1505., 1956};
_nell_ = _ell_.size();
cib::zeros(_gnu_, _nfreq_);
const double h_p = 6.62607004e-34;
const double k_B = 1.38064852e-23;
const double T_cmb = 2.725;
for (size_t i = 0; i < _nu_.size(); i++) {
const double x = h_p*_nu_[i]/k_B/T_cmb;
_gnu_[i] = x*((std::exp(x) + 1)/(std::exp(x) - 1)) - 4;
}
// Halo mass sampling
for (double m = 6; m <= 15; m+=0.1) _mh_.push_back(std::pow(10, m));
......@@ -92,7 +101,7 @@ namespace cib {
void cib_interface::compute_cls(cib::array3d & cl_cib_,
cib::array3d & cl_tsz_,
cib::array3d & cl_cibxtsz_) const
cib::array3d & cl_tszxcib_) const
{
if (! _initialized_) {
throw std::logic_error("cib_interface::update: CIB model not initialized !");
......@@ -107,90 +116,94 @@ namespace cib {
_compute_J_nu_(Jv, dj_cen, dj_sub);
const cib::array3d & u = _unfw_; // [mhalo][ell][z]
cib::array1d geo;
cib::zeros(geo, _z_.size());
for (size_t i = 0; i < _z_.size(); i++) {
geo[i] = _dchidz_[i];
geo[i] /= std::pow(_comoving_distances_[i]*(1+_z_[i]), 2);
}
cib::array1d fcxcc;
for (size_t i = 0; i < _nfreq_; i++) {
fcxcc.push_back(_fc_[i]*_cc_[i]);
}
cib::array3d y_ll;
_compute_y_ell_(y_ll);
// cib::array1d g_v;
// this->g_nu(g_v);
const double T_cmb = 2.725;
cib::array1d Kcmb_MJy = {244.1, 371.74, 483.69, 287.45, 58.04, 2.27, 1.};
cib::zeros(cl_cib_, _nfreq_, _nfreq_, _nell_);
cib::zeros(cl_tsz_, _nfreq_, _nfreq_, _nell_);
cib::zeros(cl_cibxtsz_, _nfreq_, _nfreq_, _nell_);
cib::zeros(cl_tszxcib_, _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++) {
// CIB
cib::array1d intg_mh;
cib::zeros(intg_mh, _z_.size());
cib::array1d rest2;
cib::zeros(rest2, _z_.size());
cib::array1d cib_intgz1;
cib::zeros(cib_intgz1, _z_.size());
cib::array1d cib_intgz2;
cib::zeros(cib_intgz2, _z_.size());
cib::array1d tsz_intgz1;
cib::zeros(tsz_intgz1, _z_.size());
cib::array1d tsz_intgz2;
cib::zeros(tsz_intgz2, _z_.size());
cib::array1d tszxcib_intgz1;
cib::zeros(tszxcib_intgz1, _z_.size());
cib::array1d tszxcib_intgz2;
cib::zeros(tszxcib_intgz2, _z_.size());
for (size_t z = 0; z < _z_.size(); z++) {
// One halo
cib::array1d rest1;
cib::zeros(rest1, _mh_.size());
cib::array1d cib_intgmh;
cib::zeros(cib_intgmh, _mh_.size());
cib::array1d tsz_intgmh1;
cib::zeros(tsz_intgmh1, _mh_.size());
cib::array1d tsz_intgmh2;
cib::zeros(tsz_intgmh2, _mh_.size());
cib::array1d tszxcib_intgmh1;
cib::zeros(tszxcib_intgmh1, _mh_.size());
cib::array1d tszxcib_intgmh2;
cib::zeros(tszxcib_intgmh2, _mh_.size());
cib::array1d tszxcib_intgmh3;
cib::zeros(tszxcib_intgmh3, _mh_.size());
for (size_t m = 0; m < _mh_.size(); m++) {
rest1[m] = (dj_cen[f1][m][z]*dj_sub[f2][m][z]*u[m][l][z] + dj_cen[f2][m][z] * \
dj_sub[f1][m][z]*u[m][l][z] + dj_sub[f1][m][z]*dj_sub[f2][m][z] * \
std::pow(u[m][l][z], 2))/_hmfmz_[m][z];
// CIB
cib_intgmh[m] = (dj_cen[f1][m][z]*dj_sub[f2][m][z]*u[m][l][z] + dj_cen[f2][m][z] *
dj_sub[f1][m][z]*u[m][l][z] + dj_sub[f1][m][z]*dj_sub[f2][m][z] *
std::pow(u[m][l][z], 2))/_hmfmz_[m][z];
// tSZ
tsz_intgmh1[m] = _hmfmz_[m][z]*std::pow(y_ll[l][m][z], 2);
tsz_intgmh2[m] = _hmfmz_[m][z]*y_ll[l][m][z]*_biasmz_[m][z];
// tSZxCIB
const double cosm = (1+_z_[z])*std::pow(_comoving_distances_[z], 2);
tszxcib_intgmh1[m] = y_ll[l][m][z]*((dj_cen[f2][m][z]+dj_sub[f2][m][z]*u[m][l][z])*_gnu_[f1] +
(dj_cen[f1][m][z]+dj_sub[f1][m][z]*u[m][l][z])*_gnu_[f2]) / cosm;
const double bhmf = _biasmz_[m][z]*_hmfmz_[m][z];
tszxcib_intgmh2[m] = y_ll[l][m][z]*bhmf;
tszxcib_intgmh3[m] = ((dj_cen[f2][m][z]+dj_sub[f2][m][z]*u[m][l][z])*_gnu_[f1] +
(dj_cen[f1][m][z]+dj_sub[f1][m][z]*u[m][l][z])*_gnu_[f2])/cosm/_hmfmz_[m][z]*bhmf;
}
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 geo = _dchidz_[z] / std::pow(_comoving_distances_[z]*(1+_z_[z]), 2);
cib_intgz1[z] = cib::integration_simps(cib_intgmh, dm) * geo;
cib_intgz2[z] = Jv[f2][z][l]*Jv[f1][z][l]*geo*_power_[l][z];
const double dlogm = std::log10(_m500c_[1][0]/_m500c_[0][0]);
tsz_intgz1[z] = cib::integration_simps(tsz_intgmh1, dlogm) * _dVcdz_[z];
tsz_intgz2[z] = std::pow(cib::integration_simps(tsz_intgmh2, dlogm), 2) * _dVcdz_[z] * _power_[l][z];
tszxcib_intgz1[z] = cib::integration_simps(tszxcib_intgmh1, dm) * _dVcdz_[z];
tszxcib_intgz2[z] = cib::integration_simps(tszxcib_intgmh2, dm)
* cib::integration_simps(tszxcib_intgmh3, dm) * _dVcdz_[z] * _power_[l][z];
}
// CIB
{
const double intg_z = cib::integration_simps(intg_mh, &(_z_));
cl_cib_[f1][f2][l] += fcxcc[f1]*intg_z*fcxcc[f2];
cl_cib_[f1][f2][l] += fcxcc[f1]*cib::integration_simps(cib_intgz1, &(_z_))*fcxcc[f2];
cl_cib_[f1][f2][l] += fcxcc[f1]*cib::integration_simps(cib_intgz2, &(_z_))*fcxcc[f2];
}
// tSZ
{
const double intg_z = cib::integration_simps(rest2, &(_z_));
cl_cib_[f1][f2][l] += fcxcc[f1]*intg_z*fcxcc[f2];
cl_tsz_[f1][f2][l] += T_cmb*T_cmb*_gnu_[f1]*_gnu_[f2]*cib::integration_simps(tsz_intgz1, &(_z_));
cl_tsz_[f1][f2][l] += T_cmb*T_cmb*_gnu_[f1]*_gnu_[f2]*cib::integration_simps(tsz_intgz2, &(_z_));
}
// 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
// tSZxCIB
{
cl_tszxcib_[f1][f2][l] += T_cmb*cib::integration_simps(tszxcib_intgz1, &(_z_))*Kcmb_MJy[f2]*1e6;
cl_tszxcib_[f1][f2][l] += T_cmb*cib::integration_simps(tszxcib_intgz2, &(_z_))*Kcmb_MJy[f2]*1e6;
}
} // end of f2 loop
......@@ -371,41 +384,40 @@ namespace cib {
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
// 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;
// }
// }
// }
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 < _ell_.size(); l++) {
for (size_t m = 0; m < _mh_.size(); m++) {
for (size_t z = 0; z < _z_.size(); z++) {
const double cd = _comoving_distances_[z];
const double l500 = cd/(1+_z_[z])/r500[m][z];
cib::array1d integral;
cib::zeros(integral, _x_.size());
for (size_t x = 0; x < _x_.size(); x++) {
const double Pex2 = Pe[m][z][x]*std::pow(_x_[x], 2);
const double x_ls = _x_[x]/l500;
integral[x] = Pex2 * std::sin(_ell_[l]*x_ls)/(_ell_[l]*x_ls);
}
const double dx = std::log10(_x_[1]/_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 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);
const double r3 = 3*_m500c_[m][z]/(4*CLHEP::pi*delta_h_tsz*rho_crit);
r_delta_[m][z] = std::pow(r3, 1./3.);
}
}
......@@ -444,5 +456,4 @@ namespace cib {
}
} // end of namespace cib
......@@ -37,7 +37,7 @@ namespace cib {
/// Compute and retrieve Cls
void compute_cls(cib::array3d & cl_cib_,
cib::array3d & cl_tsz_,
cib::array3d & cl_cibxtsz_) const;
cib::array3d & cl_tszxcib_) const;
private:
......@@ -78,6 +78,7 @@ namespace cib {
/// Frequencies
cib::array1d _nu_;
size_t _nfreq_;
cib::array1d _gnu_;
/// Multipole
cib::array1d _ell_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment