mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
improvements to WDM transfer function parameterisations, now support for Bode et al. 2001 and Viel et al. 2005 parameterisations, fixed missing factor of h in old Bode+2001 parameterisation
This commit is contained in:
Oliver Hahn
parent
41aca84102
commit
77183c6e20
1 changed files with 64 additions and 302 deletions
|
|
@ -199,25 +199,81 @@ public:
|
|||
};
|
||||
|
||||
|
||||
|
||||
#include <map>
|
||||
class transfer_eisenstein_wdm_plugin : public transfer_eisenstein_plugin
|
||||
{
|
||||
protected:
|
||||
using transfer_eisenstein_plugin::TFfit_onek;
|
||||
real_t m_WDMalpha;
|
||||
double omegam_, wdmm_, wdmgx_, H0_, omegab_;
|
||||
double omegam_, wdmm_, wdmgx_, wdmnu_, H0_, omegab_;
|
||||
std::string type_;
|
||||
std::map< std::string, int > typemap_;
|
||||
|
||||
enum wdmtyp { wdm_bode, wdm_viel, wdm_bode_wrong=99};
|
||||
|
||||
public:
|
||||
transfer_eisenstein_wdm_plugin( config_file &cf )
|
||||
: transfer_eisenstein_plugin( cf )
|
||||
{
|
||||
typemap_.insert( std::pair<std::string,int>( "BODE", wdm_bode ) );
|
||||
typemap_.insert( std::pair<std::string,int>( "VIEL", wdm_viel ) ); // add the other types
|
||||
typemap_.insert( std::pair<std::string,int>( "BODE_WRONG", wdm_bode_wrong ) ); // add the other types
|
||||
|
||||
|
||||
|
||||
omegam_ = cf.getValue<double>("cosmology","Omega_m");
|
||||
omegab_ = cf.getValue<double>("cosmology","Omega_b");
|
||||
wdmm_ = cf.getValue<double>("cosmology","WDMmass");
|
||||
wdmgx_ = cf.getValue<double>("cosmology","WDMg_x");
|
||||
H0_ = cf.getValue<double>("cosmology","H0");
|
||||
m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15)
|
||||
*pow(H0_/0.65,1.3)*pow(wdmm_,-1.15)
|
||||
*pow(1.5/wdmgx_,0.29);
|
||||
|
||||
|
||||
H0_ = cf.getValue<double>("cosmology","H0");
|
||||
type_ = cf.getValueSafe<std::string>("cosmology","WDMtftype","BODE");
|
||||
|
||||
//type_ = std::string( toupper( type_.c_str() ) );
|
||||
|
||||
if( typemap_.find( type_ ) == typemap_.end() )
|
||||
throw std::runtime_error("unknown transfer function fit for WDM");
|
||||
|
||||
m_WDMalpha = 1.0;
|
||||
|
||||
switch( typemap_[type_] )
|
||||
{
|
||||
//... parameterisation from Bode et al. (2001), ApJ, 556, 93
|
||||
case wdm_bode:
|
||||
wdmnu_ = cf.getValueSafe<double>("cosmology","WDMnu",1.0);
|
||||
wdmgx_ = cf.getValueSafe<double>("cosmology","WDMg_x",1.5);
|
||||
m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15)
|
||||
*pow(H0_*0.01/0.65,1.3)*pow(wdmm_,-1.15)
|
||||
*pow(1.5/wdmgx_,0.29);
|
||||
|
||||
break;
|
||||
|
||||
//... parameterisation from Viel et al. (2005), Phys Rev D, 71
|
||||
case wdm_viel:
|
||||
wdmnu_ = cf.getValueSafe<double>("cosmology","WDMnu",1.12);
|
||||
m_WDMalpha = 0.049 * pow( omegam_/0.25,0.11)
|
||||
*pow(H0_*0.01/0.7,1.22)*pow(wdmm_,-1.11);
|
||||
break;
|
||||
|
||||
|
||||
//.... below is for historical reasons due to the buggy parameterisation
|
||||
//.... in early versions of MUSIC, but apart from H instead of h, Bode et al.
|
||||
case wdm_bode_wrong:
|
||||
wdmnu_ = cf.getValueSafe<double>("cosmology","WDMnu",1.0);
|
||||
wdmgx_ = cf.getValueSafe<double>("cosmology","WDMg_x",1.5);
|
||||
m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15)
|
||||
*pow(H0_/0.65,1.3)*pow(wdmm_,-1.15)
|
||||
*pow(1.5/wdmgx_,0.29);
|
||||
break;
|
||||
|
||||
default:
|
||||
wdmnu_ = cf.getValueSafe<double>("cosmology","WDMnu",1.0);
|
||||
wdmgx_ = cf.getValueSafe<double>("cosmology","WDMg_x",1.5);
|
||||
m_WDMalpha = 0.05 * pow( omegam_/0.4,0.15)
|
||||
*pow(H0_*0.01/0.65,1.3)*pow(wdmm_,-1.15)
|
||||
*pow(1.5/wdmgx_,0.29);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
inline real_t compute( real_t k, tf_type type )
|
||||
|
|
@ -230,307 +286,13 @@ public:
|
|||
|
||||
double tf = fb*tfb+fc*tfcdm;
|
||||
|
||||
return tf*pow(1.0+(m_WDMalpha*k)*(m_WDMalpha*k),-5.0);
|
||||
return tf*pow(1.0+pow(m_WDMalpha*k,2.0*wdmnu_),-5.0/wdmnu_);
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
#if 0
|
||||
|
||||
class TransferFunction_EisensteinNeutrino : public TransferFunction
|
||||
{
|
||||
/* Fitting Formulae for CDM + Baryon + Massive Neutrino (MDM) cosmologies. */
|
||||
/* Daniel J. Eisenstein & Wayne Hu, Institute for Advanced Study */
|
||||
|
||||
/* There are two primary routines here, one to set the cosmology, the
|
||||
other to construct the transfer function for a single wavenumber k.
|
||||
You should call the former once (per cosmology) and the latter as
|
||||
many times as you want. */
|
||||
|
||||
/* TFmdm_set_cosm() -- User passes all the cosmological parameters as
|
||||
arguments; the routine sets up all of the scalar quantites needed
|
||||
computation of the fitting formula. The input parameters are:
|
||||
1) omega_matter -- Density of CDM, baryons, and massive neutrinos,
|
||||
in units of the critical density.
|
||||
2) omega_baryon -- Density of baryons, in units of critical.
|
||||
3) omega_hdm -- Density of massive neutrinos, in units of critical
|
||||
4) degen_hdm -- (Int) Number of degenerate massive neutrino species
|
||||
5) omega_lambda -- Cosmological constant
|
||||
6) hubble -- Hubble constant, in units of 100 km/s/Mpc
|
||||
7) redshift -- The redshift at which to evaluate */
|
||||
|
||||
/* TFmdm_onek_mpc() -- User passes a single wavenumber, in units of Mpc^-1.
|
||||
Routine returns the transfer function from the Eisenstein & Hu
|
||||
fitting formula, based on the cosmology currently held in the
|
||||
internal variables. The routine returns T_cb (the CDM+Baryon
|
||||
density-weighted transfer function), although T_cbn (the CDM+
|
||||
Baryon+Neutrino density-weighted transfer function) is stored
|
||||
in the global variable tf_cbnu. */
|
||||
|
||||
/* We also supply TFmdm_onek_hmpc(), which is identical to the previous
|
||||
routine, but takes the wavenumber in units of h Mpc^-1. */
|
||||
|
||||
/* We hold the internal scalar quantities in global variables, so that
|
||||
the user may access them in an external program, via "extern" declarations. */
|
||||
|
||||
/* Please note that all internal length scales are in Mpc, not h^-1 Mpc! */
|
||||
|
||||
/* -------------------------- Prototypes ----------------------------- */
|
||||
|
||||
/* int TFmdm_set_cosm(float omega_matter, float omega_baryon, float omega_hdm,
|
||||
int degen_hdm, float omega_lambda, float hubble, float redshift);
|
||||
float TFmdm_onek_mpc(float kk);
|
||||
float TFmdm_onek_hmpc(float kk);
|
||||
*/
|
||||
|
||||
/* ------------------------- Global Variables ------------------------ */
|
||||
|
||||
/* The following are set in TFmdm_set_cosm() */
|
||||
float alpha_gamma, /* sqrt(alpha_nu) */
|
||||
alpha_nu, /* The small-scale suppression */
|
||||
beta_c, /* The correction to the log in the small-scale */
|
||||
num_degen_hdm, /* Number of degenerate massive neutrino species */
|
||||
f_baryon, /* Baryon fraction */
|
||||
f_bnu, /* Baryon + Massive Neutrino fraction */
|
||||
f_cb, /* Baryon + CDM fraction */
|
||||
f_cdm, /* CDM fraction */
|
||||
f_hdm, /* Massive Neutrino fraction */
|
||||
growth_k0, /* D_1(z) -- the growth function as k->0 */
|
||||
growth_to_z0, /* D_1(z)/D_1(0) -- the growth relative to z=0 */
|
||||
hhubble, /* Need to pass Hubble constant to TFmdm_onek_hmpc() */
|
||||
k_equality, /* The comoving wave number of the horizon at equality*/
|
||||
obhh, /* Omega_baryon * hubble^2 */
|
||||
omega_curv, /* = 1 - omega_matter - omega_lambda */
|
||||
omega_lambda_z, /* Omega_lambda at the given redshift */
|
||||
omega_matter_z, /* Omega_matter at the given redshift */
|
||||
omhh, /* Omega_matter * hubble^2 */
|
||||
onhh, /* Omega_hdm * hubble^2 */
|
||||
p_c, /* The correction to the exponent before drag epoch */
|
||||
p_cb, /* The correction to the exponent after drag epoch */
|
||||
sound_horizon_fit, /* The sound horizon at the drag epoch */
|
||||
theta_cmb, /* The temperature of the CMB, in units of 2.7 K */
|
||||
y_drag, /* Ratio of z_equality to z_drag */
|
||||
z_drag, /* Redshift of the drag epoch */
|
||||
z_equality; /* Redshift of matter-radiation equality */
|
||||
|
||||
/* The following are set in TFmdm_onek_mpc() */
|
||||
float gamma_eff, /* Effective \Gamma */
|
||||
growth_cb, /* Growth factor for CDM+Baryon perturbations */
|
||||
growth_cbnu, /* Growth factor for CDM+Baryon+Neutrino pert. */
|
||||
max_fs_correction, /* Correction near maximal free streaming */
|
||||
qq, /* Wavenumber rescaled by \Gamma */
|
||||
qq_eff, /* Wavenumber rescaled by effective Gamma */
|
||||
qq_nu, /* Wavenumber compared to maximal free streaming */
|
||||
tf_master, /* Master TF */
|
||||
tf_sup, /* Suppressed TF */
|
||||
y_freestream; /* The epoch of free-streaming for a given scale */
|
||||
|
||||
/* Finally, TFmdm_onek_mpc() and TFmdm_onek_hmpc() give their answers as */
|
||||
float tf_cb, /* The transfer function for density-weighted
|
||||
CDM + Baryon perturbations. */
|
||||
tf_cbnu; /* The transfer function for density-weighted
|
||||
CDM + Baryon + Massive Neutrino perturbations. */
|
||||
|
||||
/* By default, these functions return tf_cb */
|
||||
|
||||
/* ------------------------- TFmdm_set_cosm() ------------------------ */
|
||||
int TFmdm_set_cosm(float omega_matter, float omega_baryon, float omega_hdm,
|
||||
int degen_hdm, float omega_lambda, float hubble, float redshift)
|
||||
/* This routine takes cosmological parameters and a redshift and sets up
|
||||
all the internal scalar quantities needed to compute the transfer function. */
|
||||
/* INPUT: omega_matter -- Density of CDM, baryons, and massive neutrinos,
|
||||
in units of the critical density. */
|
||||
/* omega_baryon -- Density of baryons, in units of critical. */
|
||||
/* omega_hdm -- Density of massive neutrinos, in units of critical */
|
||||
/* degen_hdm -- (Int) Number of degenerate massive neutrino species */
|
||||
/* omega_lambda -- Cosmological constant */
|
||||
/* hubble -- Hubble constant, in units of 100 km/s/Mpc */
|
||||
/* redshift -- The redshift at which to evaluate */
|
||||
/* OUTPUT: Returns 0 if all is well, 1 if a warning was issued. Otherwise,
|
||||
sets many global variables for use in TFmdm_onek_mpc() */
|
||||
{
|
||||
float z_drag_b1, z_drag_b2, omega_denom;
|
||||
int qwarn;
|
||||
qwarn = 0;
|
||||
|
||||
theta_cmb = 2.728/2.7; /* Assuming T_cmb = 2.728 K */
|
||||
|
||||
/* Look for strange input */
|
||||
if (omega_baryon<0.0) {
|
||||
fprintf(stderr,
|
||||
"TFmdm_set_cosm(): Negative omega_baryon set to trace amount.\n");
|
||||
qwarn = 1;
|
||||
}
|
||||
if (omega_hdm<0.0) {
|
||||
fprintf(stderr,
|
||||
"TFmdm_set_cosm(): Negative omega_hdm set to trace amount.\n");
|
||||
qwarn = 1;
|
||||
}
|
||||
if (hubble<=0.0) {
|
||||
fprintf(stderr,"TFmdm_set_cosm(): Negative Hubble constant illegal.\n");
|
||||
exit(1); /* Can't recover */
|
||||
} else if (hubble>2.0) {
|
||||
fprintf(stderr,"TFmdm_set_cosm(): Hubble constant should be in units of 100 km/s/Mpc.\n");
|
||||
qwarn = 1;
|
||||
}
|
||||
if (redshift<=-1.0) {
|
||||
fprintf(stderr,"TFmdm_set_cosm(): Redshift < -1 is illegal.\n");
|
||||
exit(1);
|
||||
} else if (redshift>99.0) {
|
||||
fprintf(stderr,
|
||||
"TFmdm_set_cosm(): Large redshift entered. TF may be inaccurate.\n");
|
||||
qwarn = 1;
|
||||
}
|
||||
if (degen_hdm<1) degen_hdm=1;
|
||||
num_degen_hdm = (float) degen_hdm;
|
||||
/* Have to save this for TFmdm_onek_mpc() */
|
||||
/* This routine would crash if baryons or neutrinos were zero,
|
||||
so don't allow that */
|
||||
if (omega_baryon<=0) omega_baryon=1e-5;
|
||||
if (omega_hdm<=0) omega_hdm=1e-5;
|
||||
|
||||
omega_curv = 1.0-omega_matter-omega_lambda;
|
||||
omhh = omega_matter*SQR(hubble);
|
||||
obhh = omega_baryon*SQR(hubble);
|
||||
onhh = omega_hdm*SQR(hubble);
|
||||
f_baryon = omega_baryon/omega_matter;
|
||||
f_hdm = omega_hdm/omega_matter;
|
||||
f_cdm = 1.0-f_baryon-f_hdm;
|
||||
f_cb = f_cdm+f_baryon;
|
||||
f_bnu = f_baryon+f_hdm;
|
||||
|
||||
/* Compute the equality scale. */
|
||||
z_equality = 25000.0*omhh/SQR(SQR(theta_cmb)); /* Actually 1+z_eq */
|
||||
k_equality = 0.0746*omhh/SQR(theta_cmb);
|
||||
|
||||
/* Compute the drag epoch and sound horizon */
|
||||
z_drag_b1 = 0.313*pow(omhh,-0.419)*(1+0.607*pow(omhh,0.674));
|
||||
z_drag_b2 = 0.238*pow(omhh,0.223);
|
||||
z_drag = 1291*pow(omhh,0.251)/(1.0+0.659*pow(omhh,0.828))*
|
||||
(1.0+z_drag_b1*pow(obhh,z_drag_b2));
|
||||
y_drag = z_equality/(1.0+z_drag);
|
||||
|
||||
sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1.0+10.0*pow(obhh,0.75));
|
||||
|
||||
/* Set up for the free-streaming & infall growth function */
|
||||
p_c = 0.25*(5.0-sqrt(1+24.0*f_cdm));
|
||||
p_cb = 0.25*(5.0-sqrt(1+24.0*f_cb));
|
||||
|
||||
omega_denom = omega_lambda+SQR(1.0+redshift)*(omega_curv+
|
||||
omega_matter*(1.0+redshift));
|
||||
omega_lambda_z = omega_lambda/omega_denom;
|
||||
omega_matter_z = omega_matter*SQR(1.0+redshift)*(1.0+redshift)/omega_denom;
|
||||
growth_k0 = z_equality/(1.0+redshift)*2.5*omega_matter_z/
|
||||
(pow(omega_matter_z,4.0/7.0)-omega_lambda_z+
|
||||
(1.0+omega_matter_z/2.0)*(1.0+omega_lambda_z/70.0));
|
||||
growth_to_z0 = z_equality*2.5*omega_matter/(pow(omega_matter,4.0/7.0)
|
||||
-omega_lambda + (1.0+omega_matter/2.0)*(1.0+omega_lambda/70.0));
|
||||
growth_to_z0 = growth_k0/growth_to_z0;
|
||||
|
||||
/* Compute small-scale suppression */
|
||||
alpha_nu = f_cdm/f_cb*(5.0-2.*(p_c+p_cb))/(5.-4.*p_cb)*
|
||||
pow(1+y_drag,p_cb-p_c)*
|
||||
(1+f_bnu*(-0.553+0.126*f_bnu*f_bnu))/
|
||||
(1-0.193*sqrt(f_hdm*num_degen_hdm)+0.169*f_hdm*pow(num_degen_hdm,0.2))*
|
||||
(1+(p_c-p_cb)/2*(1+1/(3.-4.*p_c)/(7.-4.*p_cb))/(1+y_drag));
|
||||
alpha_gamma = sqrt(alpha_nu);
|
||||
beta_c = 1/(1-0.949*f_bnu);
|
||||
/* Done setting scalar variables */
|
||||
hhubble = hubble; /* Need to pass Hubble constant to TFmdm_onek_hmpc() */
|
||||
return qwarn;
|
||||
}
|
||||
|
||||
/* ---------------------------- TFmdm_onek_mpc() ---------------------- */
|
||||
|
||||
float TFmdm_onek_mpc(float kk)
|
||||
/* Given a wavenumber in Mpc^-1, return the transfer function for the
|
||||
cosmology held in the global variables. */
|
||||
/* Input: kk -- Wavenumber in Mpc^-1 */
|
||||
/* Output: The following are set as global variables:
|
||||
growth_cb -- the transfer function for density-weighted
|
||||
CDM + Baryon perturbations.
|
||||
growth_cbnu -- the transfer function for density-weighted
|
||||
CDM + Baryon + Massive Neutrino perturbations. */
|
||||
/* The function returns growth_cb */
|
||||
{
|
||||
float tf_sup_L, tf_sup_C;
|
||||
float temp1, temp2;
|
||||
|
||||
qq = kk/omhh*SQR(theta_cmb);
|
||||
|
||||
/* Compute the scale-dependent growth functions */
|
||||
y_freestream = 17.2*f_hdm*(1+0.488*pow(f_hdm,-7.0/6.0))*
|
||||
SQR(num_degen_hdm*qq/f_hdm);
|
||||
temp1 = pow(growth_k0, 1.0-p_cb);
|
||||
temp2 = pow(growth_k0/(1+y_freestream),0.7);
|
||||
growth_cb = pow(1.0+temp2, p_cb/0.7)*temp1;
|
||||
growth_cbnu = pow(pow(f_cb,0.7/p_cb)+temp2, p_cb/0.7)*temp1;
|
||||
|
||||
/* Compute the master function */
|
||||
gamma_eff =omhh*(alpha_gamma+(1-alpha_gamma)/
|
||||
(1+SQR(SQR(kk*sound_horizon_fit*0.43))));
|
||||
qq_eff = qq*omhh/gamma_eff;
|
||||
|
||||
tf_sup_L = log(2.71828+1.84*beta_c*alpha_gamma*qq_eff);
|
||||
tf_sup_C = 14.4+325/(1+60.5*pow(qq_eff,1.11));
|
||||
tf_sup = tf_sup_L/(tf_sup_L+tf_sup_C*SQR(qq_eff));
|
||||
|
||||
qq_nu = 3.92*qq*sqrt(num_degen_hdm/f_hdm);
|
||||
max_fs_correction = 1+1.2*pow(f_hdm,0.64)*pow(num_degen_hdm,0.3+0.6*f_hdm)/
|
||||
(pow(qq_nu,-1.6)+pow(qq_nu,0.8));
|
||||
tf_master = tf_sup*max_fs_correction;
|
||||
|
||||
/* Now compute the CDM+HDM+baryon transfer functions */
|
||||
tf_cb = tf_master*growth_cb/growth_k0;
|
||||
tf_cbnu = tf_master*growth_cbnu/growth_k0;
|
||||
return tf_cb;
|
||||
}
|
||||
|
||||
/* ---------------------------- TFmdm_onek_hmpc() ---------------------- */
|
||||
|
||||
float TFmdm_onek_hmpc(float kk)
|
||||
/* Given a wavenumber in h Mpc^-1, return the transfer function for the
|
||||
cosmology held in the global variables. */
|
||||
/* Input: kk -- Wavenumber in h Mpc^-1 */
|
||||
/* Output: The following are set as global variables:
|
||||
growth_cb -- the transfer function for density-weighted
|
||||
CDM + Baryon perturbations.
|
||||
growth_cbnu -- the transfer function for density-weighted
|
||||
CDM + Baryon + Massive Neutrino perturbations. */
|
||||
/* The function returns growth_cb */
|
||||
{
|
||||
return TFmdm_onek_mpc(kk*hhubble);
|
||||
}
|
||||
|
||||
|
||||
|
||||
TransferFunction_EisensteinNeutrino( Cosmology aCosm, real_t Omega_HDM, int degen_HDM, real_t Tcmb = 2.726 )
|
||||
: TransferFunction(aCosm)//, m_h0( aCosm.H0*0.01 )
|
||||
{
|
||||
TFmdm_set_cosm( aCosm.Omega_m, aCosm.Omega_b, Omega_HDM, degen_HDM, aCosm.Omega_L, aCosm.H0, aCosm.astart );
|
||||
|
||||
}
|
||||
|
||||
//! Computes the transfer function for k in Mpc/h by calling TFfit_onek
|
||||
virtual inline real_t compute( real_t k ){
|
||||
return TFmdm_onek_hmpc( k );
|
||||
|
||||
|
||||
/*real_t tfb, tfcdm, fb, fc; //, tfull
|
||||
TFfit_onek( k*m_h0, &tfb, &tfcdm );
|
||||
|
||||
fb = m_Cosmology.Omega_b/(m_Cosmology.Omega_m);
|
||||
fc = (m_Cosmology.Omega_m-m_Cosmology.Omega_b)/(m_Cosmology.Omega_m) ;
|
||||
|
||||
return fb*tfb+fc*tfcdm;*/
|
||||
//return 1.0;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
namespace{
|
||||
transfer_function_plugin_creator_concrete< transfer_eisenstein_plugin > creator("eisenstein");
|
||||
|
|
|
|||
Loading…
Reference in a new issue