Commit 9f41934f authored by Xavier Garrido's avatar Xavier Garrido
Browse files

compute comoving distances from class and update them

parent 0ec525f9
......@@ -42,9 +42,18 @@ namespace cib {
_ell_ = {187., 320., 502., 684., 890., 1158., 1505., 1956};
_nell_ = _ell_.size();
// Halo mass sampling
for (double m = 6; m <= 15; m+=0.1) _mh_.push_back(std::pow(10, m));
// 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));
}
cib::read_from_file("data/hmf_red.txt", _hmfmz_);
cib::read_from_file("data/snu.txt", _snu_eff_);
_initialized_ = true;
}
......@@ -55,6 +64,10 @@ 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));
}
}
void cib_interface::compute_cls(cib::array3d & cls_) const
......@@ -63,12 +76,11 @@ namespace cib {
throw std::logic_error("cib_interface::update: CIB model not initialized !");
}
std::cout << "cib_interface::compute_cls" << std::endl;
cib::zeros(cls_, _nfreq_, _nfreq_, _nell_);
cib::array3d dj_cen; // [nfreq][mhalo][z]
_compute_djc_dlnMh_(dj_cen);
cib::zeros(cls_, _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++) {
......@@ -88,22 +100,22 @@ namespace cib {
void cib_interface::_compute_djc_dlnMh_(cib::array3d & dj_cen_) const
{
const double fsub = 0.134;
const double KC = 1.0e-10; // Kennicutt constant for Chabroer IMF
cib::array1d mhalo;
for (size_t i = 0; i < _mh_.size(); i++) mhalo.push_back(_mh_[i]*(1-fsub));
cib::array2d sfr;
_compute_sfr_(mhalo, sfr);
// cib::zeros(dj_cen_, this->_snu_eff_.size(), this->_mh_.size(), this->_z_.size());
// for (size_t i = 0; i < this->_snu_eff_.size(); i++) {
// 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];
// 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];
// }
// }
// }
cib::zeros(dj_cen_, _nfreq_, _mh_.size(), _z_.size());
for (size_t i = 0; i < _nfreq_; i++) {
for (size_t j = 0; j < _mh_.size(); j++) {
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];
}
}
}
}
......@@ -111,16 +123,16 @@ namespace cib {
{
cib::array2d sfr_mhdot;
_compute_sfr_mhdot_(mhalo_, sfr_mhdot);
// cib::array2d mhdot;
// _compute_Mdot_(mhalo_, mhdot);
// 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++) {
// sfr_[i][j] = mhdot[i][j] * ob[i]/om[i] * sfr_mhdot[i][j];
// }
// }
cib::array2d mhdot;
_compute_Mdot_(mhalo_, mhdot);
cib::zeros(sfr_, mhalo_.size(), _z_.size());
for (size_t i = 0; i < mhalo_.size(); i++) {
for (size_t j = 0; j < _z_.size(); j++) {
const double ob = _klass_->get_Ob(_z_[i]);
const double om = _klass_->get_Om(_z_[i]);
sfr_[i][j] = mhdot[i][j] * ob/om * sfr_mhdot[i][j];
}
}
}
void cib_interface::_compute_sfr_mhdot_(const cib::array1d & mhalo_, cib::array2d & sfr_mhdot_) const
......@@ -147,21 +159,24 @@ namespace cib {
} // end of halo mass loop
}
// void cl_cib_tsz_cross::Mdot(const cib::array1d & mhalo_, cib::array2d & mdot_) const
// {
// cib::zeros(mdot_, mhalo_.size(), this->_z_.size());
// const bool use_mean = true;
// if (use_mean) {
// for (size_t i = 0; i < mhalo_.size(); i++) {
// 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);
// mdot_[i][j] = b*a;
// }
// }
// }
// }
void cib_interface::_compute_Mdot_(const cib::array1d & mhalo_, cib::array2d & mdot_) const
{
cib::zeros(mdot_, mhalo_.size(), _z_.size());
const double Om0 = _klass_->Omega_m();
const double Ode0 = _klass_->Omega_L();
const bool use_mean = true;
if (use_mean) {
for (size_t i = 0; i < mhalo_.size(); i++) {
const double b = std::pow(mhalo_[i]/1e12, 1.1);
for (size_t j = 0; j < _z_.size(); j++) {
const double a = 46.1*(1 + 1.11*_z_[j]) * std::sqrt(Om0 * std::pow(1 + _z_[j], 3) + Ode0);
mdot_[i][j] = b*a;
}
}
}
}
......
......@@ -46,6 +46,8 @@ namespace cib {
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;
private:
/// Initalized flag
......@@ -76,6 +78,13 @@ namespace cib {
/// Redshift values
cib::array1d _z_;
/// Comoving distances in Mpc
cib::array1d _comoving_distances_;
/// ?
cib::array2d _hmfmz_;
cib::array2d _snu_eff_;
};
} // 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