Commit 3196e334 authored by Xavier Garrido's avatar Xavier Garrido
Browse files

compute one halo for CIB

parent 9f41934f
......@@ -47,13 +47,13 @@ namespace cib {
// Redshift & comoving distances
cib::read_from_file("data/z.txt", _z_);
_comoving_distances_.reserve(_z_.size());
for (auto z: _z_) {
_comoving_distances_.push_back(_klass_->com_distance(z));
}
_comoving_distances_.resize(_z_.size());
_dVcdz_.resize(_z_.size());
_dchidz_.resize(_z_.size());
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_);
_initialized_ = true;
}
......@@ -64,13 +64,25 @@ namespace cib {
throw std::logic_error("cib_interface::update: CIB model not initialized !");
}
_parameters_ = parameters_;
// Update comoving distances for each z given cosmo. parameters
for (auto z: _z_) {
_comoving_distances_.push_back(_klass_->com_distance(z));
const double Om0 = _klass_->Omega_m();
const double Ode0 = _klass_->Omega_L();
const double H0 = _klass_->get_H0(); //Km.s-1.Mpc-1
for (size_t i = 0; i < _z_.size(); i++) {
const double z = _z_[i];
_comoving_distances_[i] = _klass_->com_distance(z);
const double ez = std::sqrt(Om0*std::pow(1+z, 3) + Ode0);
_dVcdz_[i] = _c_*std::pow(_comoving_distances_[i], 2) / (H0 * ez);
_dchidz_[i] = _c_*1e-3/ (H0 * ez);
}
}
void cib_interface::compute_cls(cib::array3d & cls_) const
void cib_interface::compute_cls(cib::array3d & cl_cib_,
cib::array3d & cl_tsz_,
cib::array3d & cl_cibxtsz_) const
{
if (! _initialized_) {
throw std::logic_error("cib_interface::update: CIB model not initialized !");
......@@ -79,13 +91,49 @@ namespace cib {
cib::array3d dj_cen; // [nfreq][mhalo][z]
_compute_djc_dlnMh_(dj_cen);
cib::array3d dj_sub; // [nfreq][mhalo][z]
_compute_djsub_dlnMh_(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::zeros(cls_, _nfreq_, _nfreq_, _nell_);
// cib::array3d y_ll;
// this->y_ell(y_ll);
// 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++) {
for (size_t f2 = 0; f2 < _nfreq_; f2++) {
// CIB
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[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];
}
const double dm = std::log10(_mh_[1]/_mh_[0]);
intg_mh[z] = cib::integration_simps(rest1, dm);
intg_mh[z] *= geo[z];
}
const double intg_z = cib::integration_simps(intg_mh, &(_z_));
cl_cib_[f1][f2][l] = fcxcc[f1]*intg_z*fcxcc[f2];
// tSZ
......@@ -97,6 +145,45 @@ namespace cib {
} // end of ell loop
}
// Private methods
void cib_interface::_compute_djsub_dlnMh_(cib::array3d & dj_sub_) const
{
const double fsub = 0.134;
const double KC = 1.0e-10; // Kennicutt constant for Chabroer IMF
cib::zeros(dj_sub_, _nfreq_, _mh_.size(), _z_.size());
for (size_t i = 0; i < _mh_.size(); i++) {
const double mhsub = _mh_[i]*(1-fsub);
cib::array1d ms;
_compute_msub_(mhsub, ms);
cib::array2d sfrI;
_compute_sfr_(ms, sfrI);
cib::array2d sfrII;
_compute_sfr_(cib::array1d(1, mhsub), sfrII);
cib::array1d subhmf;
_compute_subhmf_(_mh_[i], ms, subhmf);
cib::array1d integral;
cib::zeros(integral, ms.size());
cib::array1d intgn;
cib::zeros(intgn, _z_.size());
for (size_t j = 0; j < _z_.size(); j++) {
// Integrate over mass
for (size_t k = 0; k < ms.size(); k++) {
integral[k] = std::min(sfrI[k][j], sfrII[0][j]*ms[k]/mhsub)*subhmf[k]/KC;
}
const double dlnmsub = std::log10(ms[1]/ms[0]);
intgn[j] = cib::integration_simps(integral, dlnmsub);
const double cd = _comoving_distances_[j];
for (size_t m = 0; m < _nfreq_; m++) {
dj_sub_[m][i][j] = _snu_eff_[m][j]*_hmfmz_[i][j]*(1+_z_[j])*intgn[j]*std::pow(cd, 2);
}
}
} // end of loop on halo mass
}
void cib_interface::_compute_djc_dlnMh_(cib::array3d & dj_cen_) const
{
const double fsub = 0.134;
......@@ -111,8 +198,8 @@ namespace cib {
for (size_t k = 0; k < _z_.size(); k++) {
const double z = _z_[k];
const double cd = _comoving_distances_[k];
const double rest = this->_hmfmz_[j][k]*sfr[j][k]*(1+z)*std::pow(cd, 2)/KC;
dj_cen_[i][j][k] = rest*this->_snu_eff_[i][k];
const double rest = _hmfmz_[j][k]*sfr[j][k]*(1+z)*std::pow(cd, 2)/KC;
dj_cen_[i][j][k] = rest*_snu_eff_[i][k];
}
}
}
......@@ -178,6 +265,28 @@ namespace cib {
}
}
void cib_interface::_compute_msub_(const double mhalo_, cib::array1d & msub_) const
{
const double log10msub_min = 5;
if (std::log10(mhalo_) <= log10msub_min) {
std::stringstream iss;
iss << "Halo mass " << std::log10(mhalo_) << " should be greater than subhalo mass " << log10msub_min;
throw std::logic_error(iss.str());
}
const double logmh = std::log10(mhalo_);
for (double logm = log10msub_min; logm < logmh; logm += 0.1) {
msub_.push_back(std::pow(10, logm));
}
}
void cib_interface::_compute_subhmf_(const double mhalo_, const cib::array1d & ms_, cib::array1d & subhmf_) const
{
cib::zeros(subhmf_, ms_.size());
for (size_t i = 0; i < ms_.size(); i++) {
subhmf_[i] = 0.13*std::pow(ms_[i]/mhalo_, -0.7)*std::exp(-9.9*std::pow(ms_[i]/mhalo_, 2.5))*std::log(10);
}
}
} // end of namespace cib
......@@ -35,19 +35,28 @@ namespace cib {
void update(const cib::cib_parameters & parameters_);
/// Compute and retrieve Cls
void compute_cls(cib::array3d & cls_) const;
void compute_cls(cib::array3d & cl_cib_,
cib::array3d & cl_tsz_,
cib::array3d & cl_cibxtsz_) const;
private:
/// Compute djc_dlnMh
void _compute_djc_dlnMh_(cib::array3d & dj_cen_) const;
void _compute_djsub_dlnMh_(cib::array3d & dj_sub_) const;
void _compute_sfr_(const cib::array1d & mhalo_, cib::array2d & sfr_) const;
void _compute_sfr_mhdot_(const cib::array1d & mhalo_, cib::array2d & sfr_mhdot_) const;
void _compute_Mdot_(const cib::array1d & mhalo_, cib::array2d & sfr_mhdot_) const;
void _compute_msub_(const double mhalo_, cib::array1d & msub_) const;
void _compute_subhmf_(const double mhalo_, const cib::array1d & ms_, cib::array1d & subhmf_) const;
private:
/// Initalized flag
......@@ -82,9 +91,14 @@ namespace cib {
/// Comoving distances in Mpc
cib::array1d _comoving_distances_;
///
cib::array1d _dVcdz_;
cib::array1d _dchidz_;
/// ?
cib::array2d _hmfmz_;
cib::array2d _snu_eff_;
cib::array3d _unfw_;
};
} // end of namespace
......
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