/*
transfer_function.hh - This file is part of MUSIC -
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
#ifndef __TRANSFERFUNCTION_HH
#define __TRANSFERFUNCTION_HH
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "Numerics.hh"
#include "general.hh"
#include
//#define XI_SAMPLE
#include "config_file.hh"
enum tf_type{
total, cdm, baryon
};
//! Abstract base class for transfer functions
/*!
This class implements a purely virtual interface that can be
used to derive instances implementing various transfer functions.
*/
class transfer_function_plugin{
public:
Cosmology cosmo_; //!< cosmological parameter, read from config_file
config_file *pcf_; //!< pointer to config_file from which to read parameters
bool tf_distinct_; //!< bool if transfer function is distinct for baryons and DM
public:
//! constructor
transfer_function_plugin( config_file& cf )
: pcf_( &cf )
{
real_t zstart;
zstart = pcf_->getValue( "setup", "zstart" );
cosmo_.astart = 1.0/(1.0+zstart);
cosmo_.Omega_b = pcf_->getValue( "cosmology", "Omega_b" );
cosmo_.Omega_m = pcf_->getValue( "cosmology", "Omega_m" );
cosmo_.Omega_L = pcf_->getValue( "cosmology", "Omega_L" );
cosmo_.H0 = pcf_->getValue( "cosmology", "H0" );
cosmo_.sigma8 = pcf_->getValue( "cosmology", "sigma_8" );
cosmo_.nspect = pcf_->getValue( "cosmology", "nspec" );
}
//! destructor
virtual ~transfer_function_plugin(){ };
//! compute value of transfer function at waven umber
virtual double compute( double k, tf_type type) = 0;
//! return maximum wave number allowed
virtual double get_kmax( void ) = 0;
//! return minimum wave number allowed
virtual double get_kmin( void ) = 0;
//! return if transfer function is distinct for baryons and DM
bool tf_is_distinct( void )
{ return tf_distinct_; }
};
//! Implements abstract factory design pattern for transfer function plug-ins
struct transfer_function_plugin_creator
{
//! create an instance of a transfer function plug-in
virtual transfer_function_plugin * create( config_file& cf ) const = 0;
//! destroy an instance of a plug-in
virtual ~transfer_function_plugin_creator() { }
};
//! Write names of registered transfer function plug-ins to stdout
std::map< std::string, transfer_function_plugin_creator *>& get_transfer_function_plugin_map();
void print_transfer_function_plugins( void );
//! Concrete factory pattern for transfer function plug-ins
template< class Derived >
struct transfer_function_plugin_creator_concrete : public transfer_function_plugin_creator
{
//! register the plug-in by its name
transfer_function_plugin_creator_concrete( const std::string& plugin_name )
{
get_transfer_function_plugin_map()[ plugin_name ] = this;
}
//! create an instance of the plug-in
transfer_function_plugin * create( config_file& cf ) const
{
return new Derived( cf );
}
};
typedef transfer_function_plugin transfer_function;
transfer_function_plugin *select_transfer_function_plugin( config_file& cf );
/**********************************************************************/
/**********************************************************************/
/**********************************************************************/
//! k-space transfer function
class TransferFunction_k
{
public:
static transfer_function *ptf_;
static real_t nspec_;
double dplus_, pnorm_, sqrtpnorm_;
static tf_type type_;
TransferFunction_k( tf_type type, transfer_function *tf, real_t nspec, real_t pnorm, real_t dplus )
: dplus_(dplus), pnorm_(pnorm)
{
ptf_ = tf;
nspec_ = nspec;
sqrtpnorm_ = sqrt( pnorm_ );
type_ = type;
std::ofstream ofs("input_powerspec.txt");
double kmin=-3, kmax=3, dk=(kmax-kmin)/100.;
for( int i=0; i<100; ++i )
{
double k = pow(10.0,kmin+i*dk);
ofs << std::setw(16) << k
<< std::setw(16) << pow(dplus_*sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,type_),2)
<< std::endl;
}
}
inline real_t compute( real_t k ) const
{
return dplus_*sqrtpnorm_*pow(k,0.5*nspec_)*ptf_->compute(k,type_);
}
};
/**********************************************************************/
/**********************************************************************/
/**********************************************************************/
#define NZERO_Q
typedef std::complex complex;
class TransferFunction_real
{
public:
double Tr0_;
real_t Tmin_, Tmax_, Tscale_;
real_t rneg_, rneg2_, kny_;
static transfer_function *ptf_;
static real_t nspec_;
protected:
real_t krgood( real_t mu, real_t q, real_t dlnr, real_t kr )
{
double krnew = kr;
complex cdgamma, zm, zp;
double arg, iarg, xneg, xpos, y;
gsl_sf_result g_a, g_p;
xpos = 0.5*(mu+1.0+q);
xneg = 0.5*(mu+1.0-q);
y = M_PI/(2.0*dlnr);
zp=complex(xpos,y);
zm=complex(xneg,y);
gsl_sf_lngamma_complex_e (zp.real(), zp.imag(), &g_a, &g_p);
zp=std::polar(exp(g_a.val),g_p.val);
real_t zpa = g_p.val;
gsl_sf_lngamma_complex_e (zm.real(), zm.imag(), &g_a, &g_p);
zm=std::polar(exp(g_a.val),g_p.val);
real_t zma = g_p.val;
arg=log(2.0/kr)/dlnr+(zpa+zma)/M_PI;
iarg=(real_t)((int)(arg + 0.5));
if( arg!=iarg )
krnew=kr*exp((arg-iarg)*dlnr);
return krnew;
}
void transform( real_t pnorm, real_t dplus, unsigned N, real_t q, std::vector& rr, std::vector& TT )
{
const double mu = 0.5;
double qmin = 1.0e-6, qmax = 1.0e+6;
q = 0.0;
N = 16384;
#ifdef NZERO_Q
q=0.4;
#endif
double kmin = qmin, kmax=qmax;
double rmin = qmin, rmax = qmax;
double k0 = exp(0.5*(log(kmax)+log(kmin)));
double r0 = exp(0.5*(log(rmax)+log(rmin)));
double L = log(rmax)-log(rmin);
double k0r0 = k0*r0;
double dlnk = L/N, dlnr = L/N;
double sqrtpnorm = sqrt(pnorm);
double dir = 1.0;
double fftnorm = 1.0/N;
fftw_complex *in, *out;
in = new fftw_complex[N];
out = new fftw_complex[N];
//... perform anti-ringing correction from Hamilton (2000)
k0r0 = krgood( mu, q, dlnr, k0r0 );
std::ofstream ofsk("transfer_k.txt");
double sum_in = 0.0;
for( unsigned i=0; icompute( k, type_ );
#ifdef FFTW3
#ifdef XI_SAMPLE
in[i][0] = dplus*dplus*sqrtpnorm*sqrtpnorm*T*T*pow(k,nspec_)*pow(k,1.5-q)*exp(-k*lambda0);//*exp(k*lambda1);
#else
in[i][0] = dplus*sqrtpnorm*T*pow(k,0.5*nspec_)*pow(k,1.5-q);//*exp(-k*lambda0);//*exp(k*lambda1);
#endif
in[i][1] = 0.0;
sum_in += in[i][0];
ofsk << std::setw(16) << k < (int)N/2 )
ii -= N;
#ifndef NZERO_Q
double y=ii*M_PI/L;
complex zp((mu+1.0)*0.5,y);
gsl_sf_result g_a, g_p;
gsl_sf_lngamma_complex_e(zp.real(), zp.imag(), &g_a, &g_p);
double arg = 2.0*(log(2.0/k0r0)*y+g_p.val);
complex cu = complex(out[i].re,out[i].im)*std::polar(1.0,arg);
out[i].re = cu.real()*fftnorm;
out[i].im = cu.imag()*fftnorm;
#else
complex x(dir*q, (double)ii*2.0*M_PI/L);
gsl_sf_result g_a, g_p;
complex g1, g2, garg, U, phase;
complex twotox = pow(complex(2.0,0.0),x);
/////////////////////////////////////////////////////////
//.. evaluate complex Gamma functions
garg = 0.5*(mu+1.0+x);
gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p);
g1 = std::polar(exp(g_a.val),g_p.val);
garg = 0.5*(mu+1.0-x);
gsl_sf_lngamma_complex_e (garg.real(), garg.imag(), &g_a, &g_p);
g2 = std::polar(exp(g_a.val),g_p.val);
/////////////////////////////////////////////////////////
//.. compute U
if( (fabs(g2.real()) < 1e-19 && fabs(g2.imag()) < 1e-19) )
{
//std::cerr << "Warning : encountered possible singularity in TransferFunction_real::transform!\n";
g1 = 1.0; g2 = 1.0;
}
U = twotox * g1 / g2;
phase = pow(complex(k0r0,0.0),complex(0.0,2.0*M_PI*(double)ii/L));
#ifdef FFTW3
complex cu = complex(out[i][0],out[i][1])*U*phase*fftnorm;
out[i][0] = cu.real();
out[i][1] = cu.imag();
#else
complex cu = complex(out[i].re,out[i].im)*U*phase*fftnorm;
out[i].re = cu.real();
out[i].im = cu.imag();
#endif
/*if( (out[i].re != out[i].re)||(out[i].im != out[i].im) )
{ std::cerr << "NaN @ i=" << i << ", U= " << U << ", phase = " << phase << ", g1 = " << g1 << ", g2 = " << g2 << std::endl;
std::cerr << "mu+1+q = " << mu+1.0+q << std::endl;
//break;
}*/
#endif
}
/*out[N/2].im = 0.0;
out[N/2+1].im = 0.0;
out[N/2+1].re = out[N/2].re;
out[N/2].im = 0.0;*/
#ifdef FFTW3
fftw_execute(ip);
#else
fftw_one(ip, out, in);
#endif
rr.assign(N,0.0);
TT.assign(N,0.0);
r0 = k0r0/k0;
for( unsigned i=0; i m_xtable,m_ytable,m_dytable;
double m_xmin, m_xmax, m_dx, m_rdx;
static tf_type type_;
public:
TransferFunction_real( double boxlength, int nfull, tf_type type, transfer_function *tf,
real_t nspec, real_t pnorm, real_t dplus, real_t rmin, real_t rmax, real_t knymax, unsigned nr )
{
real_t q = 0.8;
ptf_ = tf;
nspec_ = nspec;
type_ = type;
kny_ = knymax;
/*****************************************************************/
//... compute the FFTlog transform of the k^n T(k) kernel
std::vector r,T;
transform( pnorm, dplus, nr, q, r, T );
gsl_set_error_handler_off ();
/*****************************************************************/
//... compute T(r=0) by 3D k-space integration
{
const double REL_PRECISION=1.e-5;
gsl_integration_workspace * ws = gsl_integration_workspace_alloc (1000);
gsl_function ff;
double kmin = 2.0*M_PI/boxlength;
double kmax = nfull*M_PI/boxlength;
//... integrate 0..kmax
double a[6];
a[3] = 0.0;
a[4] = kmax;
ff.function = &call_x;
ff.params = reinterpret_cast (a);
double res, err, res2, err2;
gsl_integration_qags( &ff, a[3], a[4], 0.0, REL_PRECISION, 1000, ws, &res, &err );
if( err/res > REL_PRECISION )
std::cerr << " - Warning: no convergence in \'TransferFunction_real\', rel. error=" << err/res << std::endl;
//... integrate 0..kmin
a[3] = 0.0;
a[4] = kmin;
gsl_integration_qags( &ff, a[3], a[4], 0.0, REL_PRECISION, 1000, ws, &res2, &err2 );
if( err2/res2 > 10*REL_PRECISION )
std::cerr << " - Warning: no convergence in \'TransferFunction_real\', rel. error=" << err2/res2 << std::endl;
gsl_integration_workspace_free ( ws );
//.. get kmin..kmax
res -= res2;
res *= 8.0*dplus*sqrt(pnorm);
Tr0_ = res;
}
#if 0
{
double sum = 0.0;
unsigned long long count = 0;
int nf2 = nfull/2;
double dk = 2.0*M_PI/boxlength;
#pragma omp parallel for reduction(+:sum,count)
for( int i=0; i<=nf2; ++i )
for( int j=0; j<=nf2; ++j )
for( int k=0; k<=nf2; ++k )
{
int ii(i),jj(j),kk(k);
double rr = (double)ii*(double)ii+(double)jj*(double)jj+(double)kk*(double)kk;
double kp = sqrt(rr)*dk;
if( i==0||j==0||k==0||i==nf2||j==nf2||k==nf2 )
{
sum += pow(kp,0.5*nspec_)*ptf_->compute( kp, type_ );
count++;
}
else
{
sum += 2.0*pow(kp,0.5*nspec_)*ptf_->compute( kp, type_ );
count+=2;
}
}
Tr0_ = pow(2.0*M_PI/boxlength*nfull,3.0)*dplus*sqrt(pnorm)*sum/(double)count;
}
#endif
/*****************************************************************/
//... store as table for spline interpolation
gsl_interp_accel *accp;
gsl_spline *splinep;
std::vector xsp, ysp;
for( unsigned i=0; i rmin && r[i] < rmax )
{
xsp.push_back( log10(r[i]) );
ysp.push_back( T[i]*r[i]*r[i] );
}
}
accp = gsl_interp_accel_alloc ();
//... spline interpolation is only marginally slower here
splinep = gsl_spline_alloc (gsl_interp_akima, xsp.size() );
//... set up everything for spline interpolation
gsl_spline_init (splinep, &xsp[0], &ysp[0], xsp.size() );
//.. build lookup table using spline interpolation
m_xmin = log10(rmin);
m_xmax = log10(rmax);
m_dx = (m_xmax-m_xmin)/nr;
m_rdx = 1.0/m_dx;
for(unsigned i=0; icompute( k, type_ );
//double kmax = M_PI/100.0*64.0;
#ifdef XI_SAMPLE
return 4.0*M_PI*a[0]*a[0]*T*T*pow(k,nspec_)*k*k;
#else
return 4.0*M_PI*a[0]*T*pow(k,0.5*nspec_)*k*k;// *sin(M_PI*k/kmax)/(M_PI*k/kmax);
#endif
}
static double call_x( double kx, void *arg )
{
gsl_integration_workspace * wx = gsl_integration_workspace_alloc (1000);
double *a = (double*)arg;
//double ky = a[1], kz = a[2];
double kmin = a[3], kmax = a[4];
a[0] = kx;
gsl_function FX;
FX.function = &call_y;
FX.params = reinterpret_cast (a);
double resx, errx;
gsl_integration_qags( &FX, kmin, kmax, 1e-5, 1e-5, 1000, wx, &resx, &errx );
gsl_integration_workspace_free (wx);
return resx;
}
static double call_y( double ky, void *arg )
{
gsl_integration_workspace * wy = gsl_integration_workspace_alloc (1000);
double *a = (double*)arg;
double kmin = a[3], kmax = a[4];
a[1] = ky;
gsl_function FY;
FY.function = &call_z;
FY.params = reinterpret_cast (a);
double resy, erry;
gsl_integration_qags( &FY, kmin, kmax, 1e-5, 1e-5, 1000, wy, &resy, &erry );
gsl_integration_workspace_free (wy);
return resy;
}
static double call_z( double kz, void *arg )
{
double *a = (double*)arg;
//double kmin = a[3], kmax = a[4];
double kx = a[0], ky = a[1];
double kk = sqrt(kx*kx+ky*ky+kz*kz);
double T = ptf_->compute( kk, type_ );
return pow(kk,0.5*nspec_)*T;
}
~TransferFunction_real()
{ }
inline real_t get_grad( real_t r2 ) const
{
double r = 0.5*fast_log10(r2);
double ii = (r-m_xmin)*m_rdx;
int i = (int)ii;
i=std::max(0,i);
i=std::min(i, (int)m_xtable.size()-2);
return (real_t)m_dytable[i];
}
//... fast version
inline real_t compute_real( real_t r2 ) const
{
const double EPS = 1e-8;
const double Reps2 = EPS*EPS;
if( r2 m_tab_k;
std::vector m_tab_Tk;
void read_table( void ){
#ifdef WITH_MPI
if( MPI::COMM_WORLD.Get_rank() == 0 ){
#endif
std::cerr << " - reading tabulated transfer function from file \'"
<< m_filename << "\'\n";
std::string line;
std::ifstream ifs( m_filename.c_str() );
m_tab_k.clear();
m_tab_Tk.clear();
while( !ifs.eof() ){
getline(ifs,line);
std::stringstream ss(line);
real_t k, Tk;
ss >> k;
ss >> Tk;
m_tab_k.push_back( k );
m_tab_Tk.push_back( Tk );
}
#ifdef WITH_MPI
}
unsigned n=m_tab_k.size();
MPI::COMM_WORLD.Bcast( &n, 1, MPI_UNSIGNED, 0 );
if( MPI::COMM_WORLD.Get_rank() > 0 ){
m_tab_k.assign(n,0);
m_tab_Tk.assign(n,0);
}
#ifdef SINGLE_PRECISION
MPI::COMM_WORLD.Bcast( &m_tab_k[0], n, MPI_FLOAT, 0 );
MPI::COMM_WORLD.Bcast( &m_tab_Tk[0], n, MPI_FLOAT, 0 );
#else
MPI::COMM_WORLD.Bcast( &m_tab_k[0], n, MPI_DOUBLE, 0 );
MPI::COMM_WORLD.Bcast( &m_tab_Tk[0], n, MPI_DOUBLE, 0 );
#endif
#endif
}
public:
TransferFunction_tab( Cosmology aCosm, std::string filename )
: TransferFunction( aCosm ), m_filename( filename )
{
read_table();
}
virtual inline real_t compute( real_t k ){
return linint( k, m_tab_k, m_tab_Tk );
}
inline real_t get_kmin( void ){
return m_tab_k[1];
}
inline real_t get_kmax( void ){
return m_tab_k[m_tab_k.size()-2];
}
};
class TransferFunction_CAMB : public TransferFunction{
public:
enum TFtype{
TF_baryon, TF_cdm, TF_total
};
private:
//Cosmology m_Cosmology;
std::string m_filename_Pk, m_filename_Tk;
std::vector m_tab_k, m_tab_Tk;
Spline_interp *m_psinterp;
gsl_interp_accel *acc;
gsl_spline *spline;
void read_table( TFtype iwhich ){
#ifdef WITH_MPI
if( MPI::COMM_WORLD.Get_rank() == 0 ){
#endif
std::cerr
<< " - reading tabulated transfer function data from file \n"
<< " \'" << m_filename_Tk << "\'\n";
std::string line;
std::ifstream ifs( m_filename_Tk.c_str() );
m_tab_k.clear();
m_tab_Tk.clear();
while( !ifs.eof() ){
getline(ifs,line);
if(ifs.eof()) break;
std::stringstream ss(line);
real_t k, Tkc, Tkb, Tkg, Tkr, Tknu, Tktot;
ss >> k;
ss >> Tkc;
ss >> Tkb;
ss >> Tkg;
ss >> Tkr;
ss >> Tknu;
ss >> Tktot;
m_tab_k.push_back( log10(k) );
switch( iwhich ){
case TF_total:
m_tab_Tk.push_back( log10(Tktot) );
break;
case TF_baryon:
m_tab_Tk.push_back( log10(Tkb) );
break;
case TF_cdm:
m_tab_Tk.push_back( log10(Tkc) );
break;
}
//m_tab_Tk_cdm.push_back( Tkc );
//m_tab_Tk_baryon.push_back( Tkb );
//m_tab_Tk_tot.push_back( Tktot );
}
ifs.close();
#ifdef WITH_MPI
}
unsigned n=m_tab_k.size();
MPI::COMM_WORLD.Bcast( &n, 1, MPI_UNSIGNED, 0 );
if( MPI::COMM_WORLD.Get_rank() > 0 ){
m_tab_k.assign(n,0);
m_tab_Tk.assign(n,0);
}
#ifdef SINGLE_PRECISION
MPI::COMM_WORLD.Bcast( &m_tab_k[0], n, MPI_FLOAT, 0 );
MPI::COMM_WORLD.Bcast( &m_tab_Tk[0], n, MPI_FLOAT, 0 );
#else
MPI::COMM_WORLD.Bcast( &m_tab_k[0], n, MPI_DOUBLE, 0 );
MPI::COMM_WORLD.Bcast( &m_tab_Tk[0], n, MPI_DOUBLE, 0 );
#endif
#endif
}
public:
TransferFunction_CAMB( Cosmology aCosm, std::string filename_Tk, TFtype iwhich )
: TransferFunction( aCosm ), m_filename_Tk( filename_Tk ), m_psinterp( NULL )
{
read_table( iwhich );
//m_psinterp = new Spline_interp( m_tab_k, m_tab_Tk );
acc = gsl_interp_accel_alloc();
//spline = gsl_spline_alloc( gsl_interp_cspline, m_tab_k.size() );
spline = gsl_spline_alloc( gsl_interp_akima, m_tab_k.size() );
gsl_spline_init (spline, &m_tab_k[0], &m_tab_Tk[0], m_tab_k.size() );
}
virtual ~TransferFunction_CAMB()
{
//delete m_psinterp;
gsl_spline_free (spline);
gsl_interp_accel_free (acc);
}
virtual inline real_t compute( real_t k ){
//real_t Tkc = linint( k, m_tab_k, m_tab_Tk_cdm );
//real_t Tkb = linint( k, m_tab_k, m_tab_Tk_baryon );
real_t lk = fast_log10(k);
if(lk < m_tab_k[1] )
return 1.0;
return pow(10.0, linint( lk, m_tab_k, m_tab_Tk ));
//return pow(10.0, m_psinterp->interp(log10(k)) );
//if( k0 */
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
#endif