mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
new unified access to cosmological parameters
This commit is contained in:
Oliver Hahn
parent
8686c0365c
commit
bf6dd7e0ee
8 changed files with 197 additions and 171 deletions
|
|
@ -98,6 +98,8 @@ private:
|
|||
double t = t0;
|
||||
const double eps = 1e-10;
|
||||
|
||||
const double Omega_m( cosmo_param_["Omega_m"] ), H0( cosmo_param_["H0"] );
|
||||
|
||||
while (yy[0] < amax)
|
||||
{
|
||||
// RHS of ODEs
|
||||
|
|
@ -111,7 +113,7 @@ private:
|
|||
// d D/dtau
|
||||
dy[1] = Dprime;
|
||||
// d^2 D / dtau^2
|
||||
dy[2] = -a * H_of_a(a) * Dprime + 3.0 / 2.0 * cosmo_param_.Omega_m * std::pow(cosmo_param_.H0, 2) * D / a;
|
||||
dy[2] = -a * H_of_a(a) * Dprime + 3.0 / 2.0 * Omega_m * std::pow(H0, 2) * D / a;
|
||||
return dy;
|
||||
};
|
||||
|
||||
|
|
@ -166,16 +168,16 @@ public:
|
|||
music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
|
||||
|
||||
// set up transfer functions and compute normalisation
|
||||
transfer_function_ = std::move(select_TransferFunction_plugin(cf));
|
||||
transfer_function_ = std::move(select_TransferFunction_plugin(cf, cosmo_param_));
|
||||
transfer_function_->intialise();
|
||||
if( !transfer_function_->tf_isnormalised_ )
|
||||
cosmo_param_.pnorm = this->compute_pnorm_from_sigma8();
|
||||
cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() );
|
||||
else{
|
||||
cosmo_param_.pnorm = 1.0/Dplus_target_/Dplus_target_;
|
||||
cosmo_param_.set("pnorm", 1.0/Dplus_target_/Dplus_target_);
|
||||
auto sigma8 = this->compute_sigma8();
|
||||
music::ilog << "Measured sigma_8 for given PS normalisation is " << sigma8 << std::endl;
|
||||
}
|
||||
cosmo_param_.sqrtpnorm = std::sqrt(cosmo_param_.pnorm);
|
||||
cosmo_param_.set("sqrtpnorm", std::sqrt(cosmo_param_["pnorm"]));
|
||||
|
||||
music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons"
|
||||
<< " : " << (transfer_function_->tf_is_distinct() ? "yes" : "no") << std::endl;
|
||||
|
|
@ -252,7 +254,7 @@ public:
|
|||
<< std::setw(20) << ("delta_m(k,a=ap)")
|
||||
<< std::setw(20) << ("delta_bc(k,a=ap)")
|
||||
<< std::endl;
|
||||
double fb = cosmo_param_.Omega_b / cosmo_param_.Omega_m, fc = 1.0-fb;
|
||||
double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"];
|
||||
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.05)
|
||||
{
|
||||
double dm = this->get_amplitude(k, total) * Dplus_start_ / Dplus_target_;
|
||||
|
|
@ -279,11 +281,11 @@ public:
|
|||
inline double H_of_a(double a) const noexcept
|
||||
{
|
||||
double HH2 = 0.0;
|
||||
HH2 += cosmo_param_.Omega_r / (a * a * a * a);
|
||||
HH2 += cosmo_param_.Omega_m / (a * a * a);
|
||||
HH2 += cosmo_param_.Omega_k / (a * a);
|
||||
HH2 += cosmo_param_.Omega_DE * std::pow(a, -3. * (1. + cosmo_param_.w_0 + cosmo_param_.w_a)) * exp(-3. * (1.0 - a) * cosmo_param_.w_a);
|
||||
return cosmo_param_.H0 * std::sqrt(HH2);
|
||||
HH2 += cosmo_param_["Omega_r"] / (a * a * a * a);
|
||||
HH2 += cosmo_param_["Omega_m"] / (a * a * a);
|
||||
HH2 += cosmo_param_["Omega_k"] / (a * a);
|
||||
HH2 += cosmo_param_["Omega_DE"] * std::pow(a, -3. * (1. + cosmo_param_["w_0"] + cosmo_param_["w_a"])) * exp(-3. * (1.0 - a) * cosmo_param_["w_a"]);
|
||||
return cosmo_param_["H0"] * std::sqrt(HH2);
|
||||
}
|
||||
|
||||
//! Computes the linear theory growth factor D+, normalised to D+(a=1)=1
|
||||
|
|
@ -313,7 +315,7 @@ public:
|
|||
*/
|
||||
real_t get_vfact(real_t a) const noexcept
|
||||
{
|
||||
return f_of_a_(a) * a * H_of_a(a) / cosmo_param_.h;
|
||||
return f_of_a_(a) * a * H_of_a(a) / cosmo_param_["h"];
|
||||
}
|
||||
|
||||
//! Integrand for the sigma_8 normalization of the power spectrum
|
||||
|
|
@ -328,7 +330,7 @@ public:
|
|||
|
||||
double x = k * 8.0;
|
||||
double w = 3.0 * (sin(x) - x * cos(x)) / (x * x * x);
|
||||
static double nspect = (double)pcc->cosmo_param_.nspect;
|
||||
static double nspect = (double)pcc->cosmo_param_["n_s"];
|
||||
double tf = pcc->transfer_function_->compute(k, total);
|
||||
|
||||
//... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
|
||||
|
|
@ -347,7 +349,7 @@ public:
|
|||
|
||||
double x = k * 8.0;
|
||||
double w = 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x);
|
||||
static double nspect = static_cast<double>(pcc->cosmo_param_.nspect);
|
||||
static double nspect = static_cast<double>(pcc->cosmo_param_["n_s"]);
|
||||
double tf = pcc->transfer_function_->compute(k, total0);
|
||||
|
||||
//... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
|
||||
|
|
@ -375,18 +377,18 @@ public:
|
|||
*/
|
||||
inline real_t get_amplitude( const real_t k, const tf_type type) const
|
||||
{
|
||||
return std::pow(k, 0.5 * cosmo_param_.nspect) * transfer_function_->compute(k, type) * cosmo_param_.sqrtpnorm;// * ((type!=deltabc)? 1.0 : 1.0/Dplus_target_);
|
||||
return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, type) * cosmo_param_["sqrtpnorm"];// * ((type!=deltabc)? 1.0 : 1.0/Dplus_target_);
|
||||
}
|
||||
|
||||
inline real_t get_amplitude_phibc( const real_t k ) const
|
||||
{
|
||||
// need to multiply with Dplus_target since sqrtpnorm rescales like that
|
||||
return -std::pow(k, 0.5 * cosmo_param_.nspect-2.0) * transfer_function_->compute(k, deltabc) * cosmo_param_.sqrtpnorm * Dplus_target_;
|
||||
return -std::pow(k, 0.5 * cosmo_param_["n_s"]-2.0) * transfer_function_->compute(k, deltabc) * cosmo_param_["sqrtpnorm"] * Dplus_target_;
|
||||
}
|
||||
inline real_t get_amplitude_rhobc( const real_t k ) const
|
||||
{
|
||||
// need to multiply with Dplus_target since sqrtpnorm rescales like that
|
||||
return std::pow(k, 0.5 * cosmo_param_.nspect) * transfer_function_->compute(k, deltabc) * cosmo_param_.sqrtpnorm * Dplus_target_;
|
||||
return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, deltabc) * cosmo_param_["sqrtpnorm"] * Dplus_target_;
|
||||
}
|
||||
|
||||
//! Computes the normalization for the power spectrum
|
||||
|
|
@ -430,7 +432,7 @@ public:
|
|||
real_t compute_pnorm_from_sigma8(void)
|
||||
{
|
||||
auto measured_sigma8 = this->compute_sigma8();
|
||||
return cosmo_param_.sigma8 * cosmo_param_.sigma8 / (measured_sigma8 * measured_sigma8);
|
||||
return cosmo_param_["sigma_8"] * cosmo_param_["sigma_8"] / (measured_sigma8 * measured_sigma8);
|
||||
}
|
||||
};
|
||||
|
||||
|
|
@ -441,10 +443,10 @@ public:
|
|||
* @param mass mass scale
|
||||
* @returns jeans sound speed
|
||||
*/
|
||||
inline double jeans_sound_speed(double rho, double mass)
|
||||
{
|
||||
const double G = 6.67e-8;
|
||||
return pow(6.0 * mass / M_PI * std::sqrt(rho) * std::pow(G, 1.5), 1.0 / 3.0);
|
||||
}
|
||||
// inline double jeans_sound_speed(double rho, double mass)
|
||||
// {
|
||||
// const double G = 6.67e-8;
|
||||
// return pow(6.0 * mass / M_PI * std::sqrt(rho) * std::pow(G, 1.5), 1.0 / 3.0);
|
||||
// }
|
||||
|
||||
} // namespace cosmology
|
||||
|
|
@ -16,16 +16,18 @@
|
|||
// along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
#pragma once
|
||||
|
||||
#include <map>
|
||||
#include <string>
|
||||
|
||||
#include <physical_constants.hh>
|
||||
#include <config_file.hh>
|
||||
|
||||
namespace cosmology
|
||||
{
|
||||
//! structure for cosmological parameters
|
||||
struct parameters
|
||||
//! singleton structure for cosmological parameters
|
||||
class parameters
|
||||
{
|
||||
|
||||
double
|
||||
/*double
|
||||
Omega_m, //!< baryon+dark matter density
|
||||
Omega_b, //!< baryon matter density
|
||||
Omega_DE, //!< dark energy density (cosmological constant or parameterised)
|
||||
|
|
@ -37,6 +39,7 @@ struct parameters
|
|||
nspect, //!< long-wave spectral index (scale free is nspect=1)
|
||||
sigma8, //!< power spectrum normalization
|
||||
Tcmb, //!< CMB temperature (used to set Omega_r)
|
||||
YHe, //!< Helium fraction
|
||||
Neff, //!< effective number of neutrino species (used to set Omega_r)
|
||||
w_0, //!< dark energy equation of state parameter 1: w = w0 + a * wa
|
||||
w_a, //!< dark energy equation of state parameter 2: w = w0 + a * wa
|
||||
|
|
@ -47,73 +50,126 @@ struct parameters
|
|||
pnorm, //!< actual power spectrum normalisation factor
|
||||
sqrtpnorm, //!< sqrt of power spectrum normalisation factor
|
||||
vfact; //!< velocity<->displacement conversion factor in Zel'dovich approx.
|
||||
*/
|
||||
private:
|
||||
std::map<std::string,double> pmap_;
|
||||
|
||||
public:
|
||||
|
||||
double get( const std::string& key ) const
|
||||
{
|
||||
auto it = pmap_.find(key);
|
||||
if( it == pmap_.end() ){
|
||||
auto errmsg = std::string("Cosmological parameter \'")+key+std::string("\' does not exist in internal list.");
|
||||
music::elog << errmsg << std::endl;
|
||||
throw std::runtime_error(errmsg.c_str());
|
||||
}
|
||||
return it->second;
|
||||
}
|
||||
|
||||
void set( const std::string& key, const double value )
|
||||
{
|
||||
auto it = pmap_.find(key);
|
||||
if( it != pmap_.end() ){
|
||||
pmap_[key] = value;
|
||||
}else{
|
||||
auto errmsg = std::string("Cosmological parameter \'")+key+std::string("\' does not exist in internal list. Needs to be defaulted before it can be set!");
|
||||
music::elog << errmsg << std::endl;
|
||||
throw std::runtime_error(errmsg.c_str());
|
||||
}
|
||||
}
|
||||
|
||||
inline double operator[]( const std::string& key ) const { return this->get( key ); }
|
||||
|
||||
parameters() = delete;
|
||||
|
||||
parameters( const parameters& ) = default;
|
||||
|
||||
explicit parameters(config_file cf)
|
||||
explicit parameters( config_file& cf )
|
||||
{
|
||||
H0 = cf.get_value<double>("cosmology", "H0");
|
||||
h = H0 / 100.0;
|
||||
pmap_["H0"] = cf.get_value<double>("cosmology", "H0");
|
||||
pmap_["h"] = cf.get_value<double>("cosmology", "H0") / 100.0;
|
||||
pmap_["n_s"] = cf.get_value<double>("cosmology", "nspec");
|
||||
pmap_["Omega_b"] = cf.get_value<double>("cosmology", "Omega_b");
|
||||
pmap_["Omega_m"] = cf.get_value<double>("cosmology", "Omega_m");
|
||||
pmap_["Omega_c"] = this->get("Omega_m") - this->get("Omega_b");
|
||||
pmap_["Omega_DE"] = cf.get_value<double>("cosmology", "Omega_L");
|
||||
pmap_["w_0"] = cf.get_value_safe<double>("cosmology", "w0", -1.0);
|
||||
pmap_["w_a"] = cf.get_value_safe<double>("cosmology", "wa", 0.0);
|
||||
pmap_["Tcmb"] = cf.get_value_safe<double>("cosmology", "Tcmb", 2.7255);
|
||||
pmap_["YHe"] = cf.get_value_safe<double>("cosmology", "YHe", 0.248);
|
||||
pmap_["Neff"] = cf.get_value_safe<double>("cosmology", "Neff", 3.046);
|
||||
pmap_["sigma_8"] = cf.get_value_safe<double>("cosmology", "sigma_8",-1.0);
|
||||
pmap_["A_s"] = cf.get_value_safe<double>("cosmology", "A_s",-1.0);
|
||||
pmap_["Omega_gamma"] = 4 * phys_const::sigma_SI / std::pow(phys_const::c_SI, 3) * std::pow(this->get("Tcmb"), 4.0) / phys_const::rhocrit_h2_SI / (this->get("h") * this->get("h"));
|
||||
pmap_["Omega_nu"] = this->get("Neff") * this->get("Omega_gamma") * 7. / 8. * std::pow(4. / 11., 4. / 3.);
|
||||
pmap_["Omega_r"] = this->get("Omega_gamma") + this->get("Omega_nu");
|
||||
// H0 = cf.get_value<double>("cosmology", "H0");
|
||||
// h = H0 / 100.0;
|
||||
|
||||
nspect = cf.get_value<double>("cosmology", "nspec");
|
||||
// nspect = cf.get_value<double>("cosmology", "nspec");
|
||||
|
||||
Omega_b = cf.get_value<double>("cosmology", "Omega_b");
|
||||
// Omega_b = cf.get_value<double>("cosmology", "Omega_b");
|
||||
|
||||
Omega_m = cf.get_value<double>("cosmology", "Omega_m");
|
||||
// Omega_m = cf.get_value<double>("cosmology", "Omega_m");
|
||||
|
||||
Omega_DE = cf.get_value<double>("cosmology", "Omega_L");
|
||||
// Omega_DE = cf.get_value<double>("cosmology", "Omega_L");
|
||||
|
||||
w_0 = cf.get_value_safe<double>("cosmology", "w0", -1.0);
|
||||
// w_0 = cf.get_value_safe<double>("cosmology", "w0", -1.0);
|
||||
|
||||
w_a = cf.get_value_safe<double>("cosmology", "wa", 0.0);
|
||||
// w_a = cf.get_value_safe<double>("cosmology", "wa", 0.0);
|
||||
|
||||
Tcmb = cf.get_value_safe<double>("cosmology", "Tcmb", 2.7255);
|
||||
// Tcmb = cf.get_value_safe<double>("cosmology", "Tcmb", 2.7255);
|
||||
|
||||
Neff = cf.get_value_safe<double>("cosmology", "Neff", 3.046);
|
||||
// YHe = cf.get_value_safe<double>("cosmology", "YHe", 0.248);
|
||||
|
||||
sigma8 = cf.get_value_safe<double>("cosmology", "sigma_8",-1.0);
|
||||
// Neff = cf.get_value_safe<double>("cosmology", "Neff", 3.046);
|
||||
|
||||
// sigma8 = cf.get_value_safe<double>("cosmology", "sigma_8",-1.0);
|
||||
|
||||
// calculate energy density in ultrarelativistic species from Tcmb and Neff
|
||||
double Omega_gamma = 4 * phys_const::sigma_SI / std::pow(phys_const::c_SI, 3) * std::pow(Tcmb, 4.0) / phys_const::rhocrit_h2_SI / (h * h);
|
||||
double Omega_nu = Neff * Omega_gamma * 7. / 8. * std::pow(4. / 11., 4. / 3.);
|
||||
Omega_r = Omega_gamma + Omega_nu;
|
||||
// double Omega_gamma = 4 * phys_const::sigma_SI / std::pow(phys_const::c_SI, 3) * std::pow(Tcmb, 4.0) / phys_const::rhocrit_h2_SI / (h * h);
|
||||
// double Omega_nu = Neff * Omega_gamma * 7. / 8. * std::pow(4. / 11., 4. / 3.);
|
||||
// Omega_r = Omega_gamma + Omega_nu;
|
||||
|
||||
if (cf.get_value_safe<bool>("cosmology", "ZeroRadiation", false))
|
||||
{
|
||||
Omega_r = 0.0;
|
||||
pmap_["Omega_r"] = 0.0;
|
||||
}
|
||||
|
||||
f_b = Omega_b / Omega_m;
|
||||
pmap_["f_b"] = this->get("Omega_b") / this->get("Omega_m");
|
||||
pmap_["f_c"] = 1.0 - this->get("f_b");
|
||||
|
||||
#if 1
|
||||
// assume zero curvature, take difference from dark energy
|
||||
Omega_DE += 1.0 - Omega_m - Omega_DE - Omega_r;
|
||||
Omega_k = 0.0;
|
||||
pmap_["Omega_DE"] += 1.0 - this->get("Omega_m") - this->get("Omega_DE") - this->get("Omega_r");
|
||||
// Omega_DE += 1.0 - Omega_m - Omega_DE - Omega_r;
|
||||
pmap_["Omega_k"] = 0.0;
|
||||
#else
|
||||
// allow for curvature
|
||||
Omega_k = 1.0 - Omega_m - Omega_DE - Omega_r;
|
||||
#endif
|
||||
|
||||
dplus = 0.0;
|
||||
pnorm = 0.0;
|
||||
vfact = 0.0;
|
||||
pmap_["dplus"] = 0.0;
|
||||
pmap_["pnorm"] = 0.0;
|
||||
pmap_["sqrtpnorm"] = 0.0;
|
||||
pmap_["vfact"] = 0.0;
|
||||
|
||||
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
|
||||
music::ilog << "Cosmological parameters are: " << std::endl;
|
||||
music::ilog << " H0 = " << std::setw(16) << H0 << "sigma_8 = " << std::setw(16) << sigma8 << std::endl;
|
||||
music::ilog << " Omega_c = " << std::setw(16) << Omega_m-Omega_b << "Omega_b = " << std::setw(16) << Omega_b << std::endl;
|
||||
music::ilog << " H0 = " << std::setw(16) << this->get("H0") << "sigma_8 = " << std::setw(16) << this->get("sigma_8") << std::endl;
|
||||
music::ilog << " Omega_c = " << std::setw(16) << this->get("Omega_c") << "Omega_b = " << std::setw(16) << this->get("Omega_b") << std::endl;
|
||||
if (!cf.get_value_safe<bool>("cosmology", "ZeroRadiation", false)){
|
||||
music::ilog << " Omega_g = " << std::setw(16) << Omega_gamma << "Omega_nu = " << std::setw(16) << Omega_nu << std::endl;
|
||||
music::ilog << " Omega_g = " << std::setw(16) << this->get("Omega_gamma") << "Omega_nu = " << std::setw(16) << this->get("Omega_nu") << std::endl;
|
||||
}else{
|
||||
music::ilog << " Omega_r = " << std::setw(16) << Omega_r << std::endl;
|
||||
music::ilog << " Omega_r = " << std::setw(16) << this->get("Omega_r") << std::endl;
|
||||
}
|
||||
music::ilog << " Omega_DE = " << std::setw(16) << Omega_DE << "nspect = " << std::setw(16) << nspect << std::endl;
|
||||
music::ilog << " w0 = " << std::setw(16) << w_0 << "w_a = " << std::setw(16) << w_a << std::endl;
|
||||
music::ilog << " Omega_DE = " << std::setw(16) << this->get("Omega_DE") << "nspect = " << std::setw(16) << this->get("n_s") << std::endl;
|
||||
music::ilog << " w0 = " << std::setw(16) << this->get("w_0") << "w_a = " << std::setw(16) << this->get("w_a") << std::endl;
|
||||
|
||||
if( Omega_r > 0.0 )
|
||||
if( this->get("Omega_r") > 0.0 )
|
||||
{
|
||||
music::wlog << "Radiation enabled, using Omega_r=" << Omega_r << " internally."<< std::endl;
|
||||
music::wlog << "Radiation enabled, using Omega_r=" << this->get("Omega_r") << " internally."<< std::endl;
|
||||
music::wlog << "Make sure your sim code supports this..." << std::endl;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -21,6 +21,7 @@
|
|||
#include <memory>
|
||||
#include <general.hh>
|
||||
#include <config_file.hh>
|
||||
#include <cosmology_parameters.hh>
|
||||
|
||||
enum tf_type
|
||||
{
|
||||
|
|
@ -42,8 +43,8 @@ enum tf_type
|
|||
class TransferFunction_plugin
|
||||
{
|
||||
public:
|
||||
// Cosmology cosmo_; //!< cosmological parameter, read from config_file
|
||||
config_file *pcf_; //!< pointer to config_file from which to read parameters
|
||||
const cosmology::parameters& cosmo_params_; //!< cosmological parameters are stored here
|
||||
bool tf_distinct_; //!< bool if density transfer function is distinct for baryons and DM
|
||||
bool tf_withvel_; //!< bool if also have velocity transfer functions
|
||||
bool tf_withtotal0_; //!< have the z=0 spectrum for normalisation purposes
|
||||
|
|
@ -52,8 +53,9 @@ class TransferFunction_plugin
|
|||
|
||||
public:
|
||||
//! constructor
|
||||
TransferFunction_plugin(config_file &cf)
|
||||
: pcf_(&cf), tf_distinct_(false), tf_withvel_(false), tf_withtotal0_(false), tf_velunits_(false), tf_isnormalised_(false)
|
||||
TransferFunction_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
|
||||
: pcf_(&cf), cosmo_params_( cosmo_params ), tf_distinct_(false), tf_withvel_(false),
|
||||
tf_withtotal0_(false), tf_velunits_(false), tf_isnormalised_(false)
|
||||
{ }
|
||||
|
||||
//! destructor
|
||||
|
|
@ -100,7 +102,7 @@ class TransferFunction_plugin
|
|||
struct TransferFunction_plugin_creator
|
||||
{
|
||||
//! create an instance of a transfer function plug-in
|
||||
virtual std::unique_ptr<TransferFunction_plugin> create(config_file &cf) const = 0;
|
||||
virtual std::unique_ptr<TransferFunction_plugin> create(config_file &cf, const cosmology::parameters& cp) const = 0;
|
||||
|
||||
//! destroy an instance of a plug-in
|
||||
virtual ~TransferFunction_plugin_creator() {}
|
||||
|
|
@ -121,12 +123,12 @@ struct TransferFunction_plugin_creator_concrete : public TransferFunction_plugin
|
|||
}
|
||||
|
||||
//! create an instance of the plug-in
|
||||
std::unique_ptr<TransferFunction_plugin> create(config_file &cf) const
|
||||
std::unique_ptr<TransferFunction_plugin> create(config_file &cf, const cosmology::parameters& cp) const
|
||||
{
|
||||
return std::make_unique<Derived>(cf);
|
||||
return std::make_unique<Derived>(cf,cp);
|
||||
}
|
||||
};
|
||||
|
||||
// typedef TransferFunction_plugin TransferFunction;
|
||||
|
||||
std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf);
|
||||
std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf, const cosmology::parameters& cp);
|
||||
|
|
|
|||
|
|
@ -490,7 +490,7 @@ int run( config_file& the_config )
|
|||
music::ilog << std::endl
|
||||
<< ">>> Computing ICs for species \'" << cosmo_species_name[this_species] << "\' <<<\n" << std::endl;
|
||||
|
||||
const real_t C_species = (this_species == cosmo_species::baryon)? (1.0-the_cosmo_calc->cosmo_param_.f_b) : -the_cosmo_calc->cosmo_param_.f_b;
|
||||
const real_t C_species = (this_species == cosmo_species::baryon)? (1.0-the_cosmo_calc->cosmo_param_["f_b"]) : -the_cosmo_calc->cosmo_param_["f_b"];
|
||||
|
||||
// main loop block
|
||||
{
|
||||
|
|
|
|||
|
|
@ -20,7 +20,7 @@
|
|||
|
||||
#include <vector>
|
||||
|
||||
#include "transfer_function_plugin.hh"
|
||||
#include <transfer_function_plugin.hh>
|
||||
|
||||
const double tiny = 1e-30;
|
||||
|
||||
|
|
@ -28,6 +28,9 @@ class transfer_CAMB_file_plugin : public TransferFunction_plugin
|
|||
{
|
||||
|
||||
private:
|
||||
|
||||
using TransferFunction_plugin::cosmo_params_;
|
||||
|
||||
std::string m_filename_Pk, m_filename_Tk;
|
||||
std::vector<double> m_tab_k, m_tab_Tk_tot, m_tab_Tk_cdm, m_tab_Tk_baryon;
|
||||
std::vector<double> m_tab_Tvk_tot, m_tab_Tvk_cdm, m_tab_Tvk_baryon;
|
||||
|
|
@ -36,7 +39,7 @@ private:
|
|||
gsl_spline *spline_tot, *spline_cdm, *spline_baryon;
|
||||
gsl_spline *spline_vtot, *spline_vcdm, *spline_vbaryon;
|
||||
|
||||
double m_kmin, m_kmax, m_Omega_b, m_Omega_m, m_zstart;
|
||||
double m_kmin, m_kmax;
|
||||
unsigned m_nlines;
|
||||
|
||||
bool m_linbaryoninterp;
|
||||
|
|
@ -105,10 +108,10 @@ private:
|
|||
throw std::runtime_error("error reading transfer function file \'" + m_filename_Tk + "\'");
|
||||
}
|
||||
|
||||
if (m_Omega_b < 1e-6)
|
||||
if (cosmo_params_["Omega_b"] < 1e-6)
|
||||
Tkvtot = Tktot;
|
||||
else
|
||||
Tkvtot = ((m_Omega_m - m_Omega_b) * Tkvc + m_Omega_b * Tkvb) / m_Omega_m;
|
||||
Tkvtot = cosmo_params_["f_c"] * Tkvc + cosmo_params_["f_b"]* Tkvb;
|
||||
|
||||
m_linbaryoninterp |= Tkb < 0.0 || Tkvb < 0.0;
|
||||
|
||||
|
|
@ -180,15 +183,12 @@ private:
|
|||
}
|
||||
|
||||
public:
|
||||
transfer_CAMB_file_plugin(config_file &cf)
|
||||
: TransferFunction_plugin(cf)
|
||||
transfer_CAMB_file_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
|
||||
: TransferFunction_plugin(cf, cosmo_params)
|
||||
{
|
||||
music::wlog << "The CAMB file plugin is not well tested! Proceed with checks of correctness of output before running a simulation!" << std::endl;
|
||||
|
||||
m_filename_Tk = pcf_->get_value<std::string>("cosmology", "transfer_file");
|
||||
m_Omega_m = cf.get_value<double>("cosmology", "Omega_m"); //MvD
|
||||
m_Omega_b = cf.get_value<double>("cosmology", "Omega_b"); //MvD
|
||||
m_zstart = cf.get_value<double>("setup", "zstart"); //MvD
|
||||
|
||||
read_table();
|
||||
|
||||
|
|
@ -381,7 +381,7 @@ public:
|
|||
|
||||
case vtotal0:
|
||||
case vtotal:
|
||||
return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot)); //MvD
|
||||
return pow(10.0, gsl_spline_eval(spline_vtot, lk, acc_vtot));
|
||||
|
||||
case vcdm0:
|
||||
case vcdm:
|
||||
|
|
|
|||
|
|
@ -28,12 +28,16 @@
|
|||
#include <general.hh>
|
||||
#include <config_file.hh>
|
||||
#include <transfer_function_plugin.hh>
|
||||
#include <ic_generator.hh>
|
||||
|
||||
#include <math/interpolate.hh>
|
||||
|
||||
class transfer_CLASS_plugin : public TransferFunction_plugin
|
||||
{
|
||||
|
||||
private:
|
||||
|
||||
using TransferFunction_plugin::cosmo_params_;
|
||||
|
||||
interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_;
|
||||
interpolated_function_1d<true, true, false> delta_c0_, delta_b0_, delta_n0_, delta_m0_, theta_c0_, theta_b0_, theta_n0_, theta_m0_;
|
||||
|
||||
|
|
@ -42,7 +46,8 @@ private:
|
|||
// gsl_spline *gsl_sp_Cplus_, *gsl_sp_Cminus_;
|
||||
// std::vector<double> tab_Cplus_, tab_Cminus_;
|
||||
|
||||
double Omega_m_, Omega_b_, N_ur_, zstart_, ztarget_, kmax_, kmin_, h_, astart_, atarget_, A_s_, n_s_, sigma8_, Tcmb_, tnorm_;
|
||||
//double Omega_m_, Omega_b_, N_ur_, zstart_, ztarget_, kmax_, kmin_, h_, astart_, atarget_, A_s_, n_s_, sigma8_, Tcmb_, tnorm_;
|
||||
double zstart_, ztarget_, astart_, atarget_, kmax_, kmin_, h_, tnorm_;
|
||||
|
||||
ClassParams pars_;
|
||||
std::unique_ptr<ClassEngine> the_ClassEngine_;
|
||||
|
|
@ -70,11 +75,11 @@ private:
|
|||
add_class_parameter("gauge", "synchronous");
|
||||
|
||||
//--- cosmological parameters, densities --------------------------
|
||||
add_class_parameter("h", h_);
|
||||
add_class_parameter("h", cosmo_params_.get("h"));
|
||||
|
||||
add_class_parameter("Omega_b", Omega_b_);
|
||||
add_class_parameter("Omega_cdm", Omega_m_ - Omega_b_);
|
||||
add_class_parameter("Omega_k", 0.0);
|
||||
add_class_parameter("Omega_b", cosmo_params_.get("Omega_b"));
|
||||
add_class_parameter("Omega_cdm", cosmo_params_.get("Omega_c"));
|
||||
add_class_parameter("Omega_k", cosmo_params_.get("Omega_k"));
|
||||
add_class_parameter("Omega_fld", 0.0);
|
||||
add_class_parameter("Omega_scf", 0.0);
|
||||
|
||||
|
|
@ -88,7 +93,7 @@ private:
|
|||
#if 1
|
||||
//default off
|
||||
// add_class_parameter("Omega_ur",0.0);
|
||||
add_class_parameter("N_eff", N_ur_);
|
||||
add_class_parameter("N_eff", cosmo_params_.get("Neff"));
|
||||
add_class_parameter("N_ncdm", 0);
|
||||
|
||||
#else
|
||||
|
|
@ -102,15 +107,15 @@ private:
|
|||
//--- cosmological parameters, primordial -------------------------
|
||||
add_class_parameter("P_k_ini type", "analytic_Pk");
|
||||
|
||||
if( A_s_ > 0.0 ){
|
||||
add_class_parameter("A_s", A_s_);
|
||||
if( cosmo_params_.get("A_s") > 0.0 ){
|
||||
add_class_parameter("A_s", cosmo_params_.get("A_s"));
|
||||
}else{
|
||||
add_class_parameter("sigma8", sigma8_);
|
||||
add_class_parameter("sigma8", cosmo_params_.get("sigma_8"));
|
||||
}
|
||||
add_class_parameter("n_s", n_s_);
|
||||
add_class_parameter("n_s", cosmo_params_.get("n_s"));
|
||||
add_class_parameter("alpha_s", 0.0);
|
||||
add_class_parameter("T_cmb", Tcmb_);
|
||||
add_class_parameter("YHe", 0.248);
|
||||
add_class_parameter("T_cmb", cosmo_params_.get("Tcmb"));
|
||||
add_class_parameter("YHe", cosmo_params_.get("YHe"));
|
||||
|
||||
// additional parameters
|
||||
add_class_parameter("reio_parametrization", "reio_none");
|
||||
|
|
@ -170,14 +175,15 @@ private:
|
|||
|
||||
the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm);
|
||||
|
||||
real_t fc = (Omega_m_ - Omega_b_) / Omega_m_;
|
||||
real_t fb = Omega_b_ / Omega_m_;
|
||||
double h = cosmo_params_.get("h");
|
||||
double fc = cosmo_params_.get("f_c");
|
||||
double fb = cosmo_params_.get("f_b");
|
||||
|
||||
for (size_t i = 0; i < k.size(); ++i)
|
||||
{
|
||||
// convert to 'CAMB' format, since we interpolate loglog and
|
||||
// don't want negative numbers...
|
||||
auto ik2 = 1.0 / (k[i] * k[i]) * h_ * h_;
|
||||
auto ik2 = 1.0 / (k[i] * k[i]) * h * h;
|
||||
dc[i] = -dc[i] * ik2;
|
||||
db[i] = -db[i] * ik2;
|
||||
dn[i] = -dn[i] * ik2;
|
||||
|
|
@ -190,33 +196,30 @@ private:
|
|||
}
|
||||
|
||||
public:
|
||||
explicit transfer_CLASS_plugin(config_file &cf)
|
||||
: TransferFunction_plugin(cf)
|
||||
explicit transfer_CLASS_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
|
||||
: TransferFunction_plugin(cf,cosmo_params)
|
||||
{
|
||||
this->tf_isnormalised_ = true;
|
||||
|
||||
ofs_class_input_.open("input_class_parameters.ini", std::ios::trunc);
|
||||
|
||||
h_ = pcf_->get_value<double>("cosmology", "H0") / 100.0;
|
||||
Omega_m_ = pcf_->get_value<double>("cosmology", "Omega_m");
|
||||
Omega_b_ = pcf_->get_value<double>("cosmology", "Omega_b");
|
||||
N_ur_ = pcf_->get_value_safe<double>("cosmology", "Neff", 3.046);
|
||||
// all cosmological parameters need to be passed through the_cosmo_calc
|
||||
|
||||
ztarget_ = pcf_->get_value_safe<double>("cosmology", "ztarget", 0.0);
|
||||
atarget_ = 1.0 / (1.0 + ztarget_);
|
||||
zstart_ = pcf_->get_value<double>("setup", "zstart");
|
||||
astart_ = 1.0 / (1.0 + zstart_);
|
||||
A_s_ = pcf_->get_value_safe<double>("cosmology", "A_s", -1.0);
|
||||
n_s_ = pcf_->get_value<double>("cosmology", "nspec");
|
||||
Tcmb_ = cf.get_value_safe<double>("cosmology", "Tcmb", 2.7255);
|
||||
|
||||
if (A_s_ > 0) {
|
||||
music::ilog << "CLASS: Using A_s=" << A_s_<< " to normalise the transfer function." << std::endl;
|
||||
h_ = cosmo_params_["h"];
|
||||
|
||||
if (cosmo_params_["A_s"] > 0.0) {
|
||||
music::ilog << "CLASS: Using A_s=" << cosmo_params_["A_s"] << " to normalise the transfer function." << std::endl;
|
||||
}else{
|
||||
sigma8_ = pcf_->get_value_safe<double>("cosmology", "sigma_8", -1.0);
|
||||
if( sigma8_ < 0 ){
|
||||
double sigma8 = cosmo_params_["sigma_8"];
|
||||
if( sigma8 < 0 ){
|
||||
throw std::runtime_error("Need to specify either A_s or sigma_8 for CLASS plugin...");
|
||||
}
|
||||
music::ilog << "CLASS: Using sigma8_ =" << sigma8_<< " to normalise the transfer function." << std::endl;
|
||||
music::ilog << "CLASS: Using sigma8_ =" << sigma8<< " to normalise the transfer function." << std::endl;
|
||||
}
|
||||
|
||||
// determine highest k we will need for the resolution selected
|
||||
|
|
@ -226,11 +229,11 @@ public:
|
|||
|
||||
// initialise CLASS and get the normalisation
|
||||
this->init_ClassEngine();
|
||||
A_s_ = the_ClassEngine_->get_A_s(); // this either the input one, or the one computed from sigma8
|
||||
double A_s_ = the_ClassEngine_->get_A_s(); // this either the input one, or the one computed from sigma8
|
||||
|
||||
// compute the normalisation to interface with MUSIC
|
||||
double k_p = pcf_->get_value_safe<double>("cosmology", "k_p", 0.05);
|
||||
tnorm_ = std::sqrt(2.0 * M_PI * M_PI * A_s_ * std::pow(1.0 / k_p * h_, n_s_ - 1) / std::pow(2.0 * M_PI, 3.0));
|
||||
tnorm_ = std::sqrt(2.0 * M_PI * M_PI * A_s_ * std::pow(1.0 / k_p * cosmo_params["h"], cosmo_params["n_s"] - 1) / std::pow(2.0 * M_PI, 3.0));
|
||||
|
||||
// compute the transfer function at z=0 using CLASS engine
|
||||
std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm;
|
||||
|
|
|
|||
|
|
@ -208,6 +208,8 @@ struct eisenstein_transfer
|
|||
*/
|
||||
class transfer_eisenstein_plugin : public TransferFunction_plugin
|
||||
{
|
||||
using TransferFunction_plugin::cosmo_params_;
|
||||
|
||||
protected:
|
||||
eisenstein_transfer etf_;
|
||||
|
||||
|
|
@ -218,13 +220,13 @@ public:
|
|||
\param Tcmb mean temperature of the CMB fluctuations (defaults to
|
||||
Tcmb = 2.726 if not specified)
|
||||
*/
|
||||
transfer_eisenstein_plugin(config_file &cf)
|
||||
: TransferFunction_plugin(cf)
|
||||
transfer_eisenstein_plugin(config_file &cf, const cosmology::parameters& cp)
|
||||
: TransferFunction_plugin(cf, cp)
|
||||
{
|
||||
double Tcmb = pcf_->get_value_safe<double>("cosmology", "Tcmb", 2.726);
|
||||
double H0 = pcf_->get_value<double>("cosmology", "H0");
|
||||
double Omega_m = pcf_->get_value<double>("cosmology", "Omega_m");
|
||||
double Omega_b = pcf_->get_value<double>("cosmology", "Omega_b");
|
||||
double Tcmb = cosmo_params_["Tcmb"];
|
||||
double H0 = cosmo_params_["H0"];
|
||||
double Omega_m = cosmo_params_["Omega_m"];
|
||||
double Omega_b = cosmo_params_["Omega_b"];
|
||||
|
||||
etf_.set_parameters(H0, Omega_m, Omega_b, Tcmb);
|
||||
|
||||
|
|
@ -252,6 +254,8 @@ public:
|
|||
#include <map>
|
||||
class transfer_eisenstein_wdm_plugin : public TransferFunction_plugin
|
||||
{
|
||||
using TransferFunction_plugin::cosmo_params_;
|
||||
|
||||
protected:
|
||||
real_t m_WDMalpha, m_h0;
|
||||
double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_;
|
||||
|
|
@ -268,13 +272,13 @@ protected:
|
|||
};
|
||||
|
||||
public:
|
||||
transfer_eisenstein_wdm_plugin(config_file &cf)
|
||||
: TransferFunction_plugin(cf)
|
||||
transfer_eisenstein_wdm_plugin(config_file &cf, const cosmology::parameters& cp)
|
||||
: TransferFunction_plugin(cf, cp)
|
||||
{
|
||||
double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726);
|
||||
omegam_ = pcf_->get_value<double>("cosmology", "Omega_m");
|
||||
omegab_ = pcf_->get_value<double>("cosmology", "Omega_b");
|
||||
H0_ = pcf_->get_value<double>("cosmology", "H0");
|
||||
double Tcmb = cosmo_params_["Tcmb"];
|
||||
H0_ = cosmo_params_["H0"];
|
||||
omegam_ = cosmo_params_["Omega_m"];
|
||||
omegab_ = cosmo_params_["Omega_b"];
|
||||
m_h0 = H0_ / 100.0;
|
||||
wdmm_ = pcf_->get_value<double>("cosmology", "WDMmass");
|
||||
|
||||
|
|
@ -345,20 +349,22 @@ public:
|
|||
// CDM Bino type WIMP small-scale damped spectrum from Green, Hofmann & Schwarz (2004)
|
||||
class transfer_eisenstein_cdmbino_plugin : public TransferFunction_plugin
|
||||
{
|
||||
using TransferFunction_plugin::cosmo_params_;
|
||||
|
||||
protected:
|
||||
|
||||
real_t m_h0;
|
||||
double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_;
|
||||
eisenstein_transfer etf_;
|
||||
|
||||
public:
|
||||
transfer_eisenstein_cdmbino_plugin(config_file &cf)
|
||||
: TransferFunction_plugin(cf)
|
||||
transfer_eisenstein_cdmbino_plugin(config_file &cf, const cosmology::parameters& cp)
|
||||
: TransferFunction_plugin(cf, cp)
|
||||
{
|
||||
double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726);
|
||||
|
||||
omegam_ = pcf_->get_value<double>("cosmology", "Omega_m");
|
||||
omegab_ = pcf_->get_value<double>("cosmology", "Omega_b");
|
||||
H0_ = pcf_->get_value<double>("cosmology", "H0");
|
||||
double Tcmb = cosmo_params_["Tcmb"];
|
||||
H0_ = cosmo_params_["H0"];
|
||||
omegam_ = cosmo_params_["Omega_m"];
|
||||
omegab_ = cosmo_params_["Omega_b"];
|
||||
m_h0 = H0_ / 100.0;
|
||||
|
||||
etf_.set_parameters(H0_, omegam_, omegab_, Tcmb);
|
||||
|
|
@ -397,53 +403,10 @@ public:
|
|||
}
|
||||
};
|
||||
|
||||
class transfer_eisenstein_cutoff_plugin : public TransferFunction_plugin
|
||||
{
|
||||
protected:
|
||||
real_t m_h0;
|
||||
double omegam_, H0_, omegab_, mcdm_, Tkd_, kfs_, kd_;
|
||||
double Rcut_;
|
||||
eisenstein_transfer etf_;
|
||||
|
||||
public:
|
||||
transfer_eisenstein_cutoff_plugin(config_file &cf)
|
||||
: TransferFunction_plugin(cf)
|
||||
{
|
||||
double Tcmb = pcf_->get_value_safe("cosmology", "Tcmb", 2.726);
|
||||
|
||||
omegam_ = pcf_->get_value<double>("cosmology", "Omega_m");
|
||||
omegab_ = pcf_->get_value<double>("cosmology", "Omega_b");
|
||||
H0_ = pcf_->get_value<double>("cosmology", "H0");
|
||||
m_h0 = H0_ / 100.0;
|
||||
|
||||
etf_.set_parameters(H0_, omegam_, omegab_, Tcmb);
|
||||
|
||||
Rcut_ = pcf_->get_value_safe<double>("cosmology", "Rcut", 1.0);
|
||||
}
|
||||
|
||||
inline double compute(double k, tf_type type) const
|
||||
{
|
||||
//std::cerr << k << " " << Rcut_ << " " << std::exp(-k*k*Rcut_*Rcut_) << std::endl;
|
||||
double cut = std::exp(-k*k*Rcut_*Rcut_);
|
||||
// std::cerr << k << " " << cut << ", " << std::endl;
|
||||
return etf_.at_k(k) * cut;
|
||||
}
|
||||
|
||||
inline double get_kmin(void) const
|
||||
{
|
||||
return 1e-4;
|
||||
}
|
||||
|
||||
inline double get_kmax(void) const
|
||||
{
|
||||
return 1.e4;
|
||||
}
|
||||
};
|
||||
|
||||
namespace
|
||||
{
|
||||
TransferFunction_plugin_creator_concrete<transfer_eisenstein_plugin> creator("eisenstein");
|
||||
TransferFunction_plugin_creator_concrete<transfer_eisenstein_wdm_plugin> creator2("eisenstein_wdm");
|
||||
TransferFunction_plugin_creator_concrete<transfer_eisenstein_cdmbino_plugin> creator3("eisenstein_cdmbino");
|
||||
// TransferFunction_plugin_creator_concrete<transfer_eisenstein_cutoff_plugin> creator4("eisenstein_cutoff");
|
||||
} // namespace
|
||||
|
|
|
|||
|
|
@ -39,7 +39,7 @@ void print_TransferFunction_plugins()
|
|||
music::ilog << std::endl;
|
||||
}
|
||||
|
||||
std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf)
|
||||
std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_file &cf, const cosmology::parameters& cosmo_param)
|
||||
{
|
||||
std::string tfname = cf.get_value<std::string>("cosmology", "transfer");
|
||||
|
||||
|
|
@ -57,5 +57,5 @@ std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_f
|
|||
music::ilog << std::setw(32) << std::left << "Transfer function plugin" << " : " << tfname << std::endl;
|
||||
}
|
||||
|
||||
return std::move(the_TransferFunction_plugin_creator->create(cf));
|
||||
return std::move(the_TransferFunction_plugin_creator->create(cf, cosmo_param));
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in a new issue