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

get updates from Tristram

parent d1ac5cd8
This diff is collapsed.
#ifndef CIB_CL_CIB_TSZ_CROSS
#define CIB_CL_CIB_TSZ_CROSS
// Standard library
#include <string>
#include <vector>
#include"cxxsupport/cxxutils.h"
#include"cxxsupport/fitshandle.h"
// Utilities
#include "utils.h"
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:
//constuctor
cosmology( ClassEngine & klass,
const cib::array1d & z_);
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;
// const double Ob0;
// const double Ogamma0;
// const double Onu0;
// const double n;
// const double sigma8;
// const double tau;
// const double z_reion;
// const double t0;
// const double Tcmb0;
// const double Neff;
// const array1d m_nu;
void dump() const;
double dchi_dz(double z);
private:
// Default constructor can not instantiated
};
class cl_cib_tsz_cross {
public:
// Constructor
cl_cib_tsz_cross(const cib::array2d & power_,
const cib::array1d & z_,
const cib::array1d & sig_z_arr_,
const cib::array1d & mh_,
const cib::array2d & hmfmz_,
const cib::array2d & snu_eff_,
const cib::array1d & ell_,
const cib::array3d & unfw_,
ClassEngine & klass,
const cib::array2d & biasmz_,
const double Meffmax_,
const double etamax_,
const double sigmaMh_,
const double tau_,
const cib::array1d & cc_,
const cib::array1d & fc_,
const std::string & law_,
const cib::array1d & nu_,
const cib::array2d & m500c_,
const double delta_h_tsz_,
const double B_,
const cib::array1d & x_
);
void onehalo_int(cib::array3d & Cl_1h_) const;
void twohalo_int(cib::array3d & Cl_2h_) const;
void C_ell_1h(cib::array3d & cl_) const;
void C_ell_2h(cib::array3d & cl_) const;
void onehalo(cib::array3d & cl_) const;
void twohalo(cib::array3d & cl_) const;
void g_nu(cib::array1d & gnu_) const;
void y_ell(cib::array3d & y_ell_) const;
void C(cib::array2d & c_) const;
void P_e(cib::array3d & pe_) const;
void r_delta(cib::array2d & r_delta_) const;
void ell_delta(cib::array2d & ell_delta_) const;
void J_nu(cib::array3d & Jnu_) const;
void djsub_dlnMh(cib::array3d & dj_sub_) const;
void djc_dlnMh(cib::array3d & dj_cen_) const;
void sfr(const cib::array1d & mhalo_, cib::array2d & sfr_) const;
void sfr_mhdot(const cib::array1d & mhalo_, cib::array2d & sfr_mhdot_) const;
void Mdot(const cib::array1d & mhalo_, cib::array2d & mdot_) const;
void msub(const double mhalo_, cib::array1d & msub_) const;
void subhmf(const double mhalo_, const cib::array1d & ms_, cib::array1d & subhmf_) const;
private:
// Default constructor can not instantiated
cl_cib_tsz_cross();
// Data members
// cib::array2d _kk_; // [ell][z]
cib::array2d _power_; // [ell][z]
cib::array1d _z_; // [z]
// double _z_c_;
cib::array1d _sig_z_arr_; // [z]
cib::array1d _mh_; // [mhalo]
cib::array2d _hmfmz_; // [mhalo][z]
cib::array2d _snu_eff_; // [nfreq][z]
cib::array1d _ell_; // [ell]
cib::array3d _unfw_; // [mhalo][ell][z]
cib::cosmology *_cosmo_;
cib::array2d _biasmz_; // [mhalo][z]
double _Meffmax_;
double _etamax_;
double _sigmaMh_;
// double _alpha_;
// double _delta_;
// double _beta_;
// double _tau_;
cib::array1d _cc_; // [nfreq]
cib::array1d _fc_; // [nfreq]
std::string _law_;
cib::array1d _sigpow_; // [z]
cib::array1d _nu_; // [nfreq]
size_t _nfreq_;
cib::array2d _m500c_; // [mhalo][z]
double _delta_h_tsz_;
double _B_;
cib::array1d _x_; // [?]
cib::array1d _dVcdz_; // [z]
cib::array1d _dchidz_;
double Kcmb_MJy[7] = {244.1, 371.74, 483.69, 287.45, 58.04, 2.27, 1.};
};
}// end of namespace cib
#endif
// Ourself
#include "utils.h"
// Standard libraries
#include <limits>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <functional>
namespace cib {
double _basic_simps(const cib::array1d & y_, const double start_, const double stop_, const double dx_, const cib::array1d * x_)
{
double result = 0.0;
if (x_ == 0) {
for (size_t i = start_; i < stop_; i+=2) {
result += dx_/3 * (y_[i] + 4*y_[i+1] + y_[i+2]);
}
} else {
const cib::array1d & x = *x_;
for (size_t i = start_; i < stop_; i+=2) {
const double h0 = x[i+1] - x[i];
const double h1 = x[i+2] - x[i+1];
const double hsum = h0 + h1;
const double hprod = h0 * h1;
const double hdiv = h0 / h1;
result += hsum/6 * (y_[i]*(2-1/hdiv) + y_[i+1]*hsum*hsum/hprod + y_[i+2]*(2-hdiv));
}
}
return result;
}
double _integration_simps(const cib::array1d & y_, const double dx_, const cib::array1d * x_)
{
double result = 0.0;
const size_t N = y_.size();
if (N % 2 == 0) {
double last_dx = dx_;
double first_dx = dx_;
double val = 0.0;
if (x_ != 0) {
const cib::array1d & x = *x_;
last_dx = x[N-1] - x[N-2];
first_dx = x[1] - x[0];
}
val += 0.5 * last_dx * (y_[N-1] + y_[N-2]);
result += _basic_simps(y_, 0, N-3, dx_, x_);
val += 0.5 * first_dx * (y_[0] + y_[1]);
result += _basic_simps(y_, 1, N-2, dx_, x_);
result = (result + val)/2;
} else {
result = _basic_simps(y_, 0, N-2, dx_, x_);
}
return result;
}
double integration_simps(const cib::array1d & y_, const double dx_)
{
return _integration_simps(y_, dx_, 0);
}
double integration_simps(const cib::array1d & y_, const cib::array1d * x_)
{
return _integration_simps(y_, 1.0, x_);
}
double NaN()
{
return std::numeric_limits<double>::quiet_NaN();
}
void read_from_file(std::string filename, cib::array1d & array1d_) {
std::ifstream fin(filename.c_str());
std::string line;
while (std::getline(fin, line)) {
std::istringstream iss(line);
double val = 0;
while (iss >> val) {
array1d_.push_back(val);
}
}
}
void read_from_file(std::string filename, cib::array2d & array2d_) {
std::ifstream fin(filename.c_str());
std::string line;
while (std::getline(fin, line)) {
{
// New array2d line
cib::array1d dummy;
array2d_.push_back(dummy);
}
cib::array1d & newline = array2d_.back();
std::istringstream iss(line);
double val = 0;
while (iss >> val) {
newline.push_back(val);
}
}
}
void read_from_file(std::string filename, cib::array3d & array3d_) {
std::ifstream fin(filename.c_str());
std::string line;
while (std::getline(fin, line)) {
if (line.find("#") != std::string::npos) {
cib::array2d dummy;
array3d_.push_back(dummy);
continue;
}
cib::array2d & newmatrix = array3d_.back();
{
// New array2d line
cib::array1d dummy;
newmatrix.push_back(dummy);
}
cib::array1d & newline = newmatrix.back();
std::istringstream iss(line);
double val = 0;
while (iss >> val) {
newline.push_back(val);
}
}
}
void dump_shape(const cib::array1d & array1d_) {
std::cout << "Array1d shape = (" << array1d_.size() << ")" << std::endl;
}
void dump_shape(const cib::array2d & array2d_) {
std::cout << "Array2d shape = (" << array2d_.size() << ","
<< array2d_.front().size() << ")" << std::endl;
}
void dump_shape(const cib::array3d & array3d_) {
std::cout << "Array3d shape = (" << array3d_.size() << ","
<< array3d_.front().size() << ","
<< array3d_.front().front().size() << ")" << std::endl;
}
void dump(const cib::array2d & array2d_, std::ostream & out_) {
dump_shape(array2d_);
// std::cout << std::scientific;
for (size_t il = 0; il < array2d_.size(); il++) {
const cib::array1d & line = array2d_.at(il);
for (size_t ic = 0; ic < line.size(); ic++) {
out_ << line.at(ic) << " ";
}
out_ << std::endl;
}
}
void dump(const cib::array1d & array1d_, std::ostream & out_) {
dump_shape(array1d_);
// std::cout << std::scientific;
for (size_t i = 0; i < array1d_.size(); i++) {
out_ << array1d_.at(i) << " ";
}
out_ << std::endl;
}
void _allocate(array1d & array1d_, const size_t n1d_, const double val_)
{
array1d_.resize(n1d_, val_);
}
void zeros(array1d & array1d_, const size_t n1d_)
{
_allocate(array1d_, n1d_, 0.0);
}
void zeros(array2d & array2d_, const size_t n2d_, const size_t n1d_)
{
array2d_.resize(n2d_);
for (size_t i = 0; i < n2d_; i++) {
zeros(array2d_[i], n1d_);
}
}
void zeros(array3d & array3d_, const size_t n3d_, const size_t n2d_, const size_t n1d_)
{
array3d_.resize(n3d_);
for (size_t i = 0; i < n3d_; i++) {
zeros(array3d_[i], n2d_, n1d_);
}
}
void ones(array1d & array1d_, const size_t n1d_)
{
_allocate(array1d_, n1d_, 1.0);
}
void ones(array2d & array2d_, const size_t n2d_, const size_t n1d_)
{
array2d_.resize(n2d_);
for (size_t i = 0; i < n2d_; i++) {
ones(array2d_[i], n1d_);
}
}
void ones(array3d & array3d_, const size_t n3d_, const size_t n2d_, const size_t n1d_)
{
array3d_.resize(n3d_);
for (size_t i = 0; i < n3d_; i++) {
ones(array3d_[i], n2d_, n1d_);
}
}
}// end of namespace cib
#ifndef CIB_UTILS
#define CIB_UTILS
// Standard libraries
#include <string>
#include <vector>
#include <iostream>
namespace cib {
// Singleton
template <typename T>
class singleton {
public:
static T& get_instance() {
static T instance;
return instance;
}
protected:
// derived class can call ctor and dtor
singleton() {}
~singleton() {}
private:
// no one should do copies
singleton(const singleton&);
singleton& operator=(const singleton&);
};
double NaN();
// Define aliases to vector/matrix objects
typedef std::vector<double> array1d;
typedef std::vector<std::vector<double> > array2d;
typedef std::vector<std::vector<std::vector<double> > > array3d;
// Array utilities
void read_from_file(std::string filename, cib::array1d & array1d);
void read_from_file(std::string filename, cib::array2d & array2d);
void read_from_file(std::string filename, cib::array3d & array3d);
void dump_shape(const cib::array1d & array1d);
void dump_shape(const cib::array2d & array2d);
void dump_shape(const cib::array3d & array3d);
void dump(const cib::array1d & array1d, std::ostream & out_ = std::clog);
void dump(const cib::array2d & array2d, std::ostream & out_ = std::clog);
void zeros(array1d & array1d_, const size_t n1d_);
void zeros(array2d & array2d_, const size_t n2d_, const size_t n1d_);
void zeros(array3d & array3d_, const size_t n3d_, const size_t n2d_, const size_t n1d_);
void ones(array1d & array1d_, const size_t n1d_);
void ones(array2d & array2d_, const size_t n2d_, const size_t n1d_);
void ones(array3d & array3d_, const size_t n3d_, const size_t n2d_, const size_t n1d_);
double integration_simps(const cib::array1d & y_,
const double dx_ = 1.0);
double integration_simps(const cib::array1d & y_,
const cib::array1d * x_ = 0);
// double _basic_simps(const cib::array1d & y_,
// const double start_,
// const double stop_,
// const double dx_ = 1.0,
// const cib::array1d * x_ = 0);
// void _allocate(array1d & array1d_, const size_t n1d_, const double val_);
} //end of namespace sib
#endif
......@@ -653,6 +653,7 @@ double ClassEngine::get_Fz(double z)
}
double ClassEngine::get_H(double z)
//km.s-1.Mpc-1
{
double tau;
int index;
......@@ -704,7 +705,7 @@ double ClassEngine::get_Da(double z)
}
double ClassEngine::get_DMod(double z){
double ClassEngine::get_DMod(double z){
// returns distance modulus at redshift z
double tau;
int index;
......@@ -719,6 +720,61 @@ double ClassEngine::get_DMod(double z){
}
double ClassEngine::critical_density( double z)
//Units are Msun/Mpc3
{
double tau;
int index;
//transform redshift in conformal time
background_tau_of_z(&ba,z,&tau);
//pvecback must be allocated
background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
const double G = 4.30091e-3; //pc.Msun-1.(km/s)^2
double to_Msun_per_Mpc3 = 3.*_c_*_c_/8./M_PI/G; //Msun/Mpc.(m/s)^2
double rhoc = pvecback[ba.index_bg_rho_crit]; //units are 3c^2/8piG
return( rhoc * to_Msun_per_Mpc3);
}
double ClassEngine::get_Om(double z)
{
double tau;
int index;
//transform redshift in conformal time
background_tau_of_z(&ba,z,&tau);
//pvecback must be allocated
background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
double om = pvecback[ba.index_bg_Omega_m];
return(om);
}
double ClassEngine::get_Ob(double z)
{
double tau;
int index;
//transform redshift in conformal time
background_tau_of_z(&ba,z,&tau);
//pvecback must be allocated
background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
double rhob = pvecback[ba.index_bg_rho_b];
double rhoc = pvecback[ba.index_bg_rho_crit];
return(rhob/rhoc);
}
bool
ClassEngine::get_Pkvec(const std::vector<double>& knodes, const std::vector<double>& zvec, std::vector<double>& pks)
{
......@@ -729,14 +785,17 @@ ClassEngine::get_Pkvec(const std::vector<double>& knodes, const std::vector<doub
}
bool
ClassEngine::get_PklinNodes(std::vector<double>& knodes,std::vector<double>& pknodes){
ClassEngine::get_PklinNodes(std::vector<double>& knodes,std::vector<double>& pknodes,int eta){
knodes.resize(sp.ln_k_size);
pknodes.resize(sp.ln_k_size);
int index_mode=0;
int index_ic1_ic2=0;
int index_eta = sp.ln_tau_size-1;
int index_eta = eta == 0 ? sp.ln_tau_size-1 : 0;
// for (int i=0; i < sp.ln_tau_size; i++)
// cout << i << " ln(tau):" << sp.ln_tau[i] << " tau:" << std::exp(sp.ln_tau[i]) << endl;
for (int index_k=0; index_k < sp.ln_k_size; index_k++){
knodes[index_k]=std::exp(sp.ln_k[index_k]);
......
......@@ -27,6 +27,7 @@ extern "C"{
#include<utility>
#include<ostream>
#include<sstream>
#include <CLHEP/Units/PhysicalConstants.h>
using std::string;
......@@ -149,9 +150,12 @@ public:
double get_Da(double z);
double get_sigma8(double z);
double get_f(double z);
double get_Ob(double z);
double get_Om(double z);
double critical_density( double z);
double get_Fz(double z);
double get_H(double z);
double get_H(double z); //Km.s-1.Mpc-1
double get_DMod(double z);
......@@ -161,7 +165,7 @@ public:
//nodes used by class to perfrom spline interpolation (at z=0)
//knodes/pknodes resized
bool get_PklinNodes(std::vector<double>& knodes,std::vector<double>& Pknodes);