Page MenuHomec4science
Files F80560112

cosmo_generate_power_spectum.py
No OneTemporary
Actions

File Metadata

Created
Sun, Sep 1, 09:39

cosmo_generate_power_spectum.py
View Options

#!/usr/bin/env python
import sys,os,string
import argparse
import numpy as np
####################################################################
# option parser
####################################################################
description=""
epilog =""""""
parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('-f',
action="store",
type=str,
dest="file",
default = None,
nargs=1,
help="the reference power spectrum file")
parser.add_argument("--redshift",
action="store",
dest="redshift",
metavar='FLOAT',
type=float,
default=0.0,
help='redshift')
parser.add_argument("--DMmass",
action="store",
dest="DMmass",
metavar='FLOAT',
type=float,
default=2.0,
help='dark matter particle mass')
parser.add_argument("--fWDM",
action="store",
dest="fWDM",
metavar='FLOAT',
type=float,
default=1.0,
help='warm dark matter fraction')
parser.add_argument("-o","--outputfile",
action="store",
type=str,
dest="outputfile",
default = None,
help="output file name")
parser.add_argument("--Om",
action="store",
dest="Om",
metavar='FLOAT',
type=float,
default=0.315,
help='Omega matter')
parser.add_argument("--Ob",
action="store",
dest="Ob",
metavar='FLOAT',
type=float,
default=0.049,
help='Omega baryon')
parser.add_argument("--h0",
action="store",
dest="h0",
metavar='FLOAT',
type=float,
default=0.674,
help='Hubble parameter')
parser.add_argument("--As",
action="store",
dest="As",
metavar='FLOAT',
type=float,
default=2.07e-9,
help='primordial comoving curvature power spectrum amplitude')
parser.add_argument("--ns",
action="store",
dest="ns",
metavar='FLOAT',
type=float,
default=0.965,
help='scalar spectral index')
def power_spectrum(cosmopars):
"""
Calculate power spectrum with CLASS
"""
from classy import Class
Om, Ob, As, h0, ns, mWDM, fWDM = cosmopars
if (fWDM>0):
m_nu = 4.43*mWDM**(4./3)*((Om-Ob)*h0**2/0.1225)**(-1./3)
CLASSparams = {
'h': h0,
'T_cmb': 2.726,
'Omega_b': Ob,
'Omega_cdm': (1-fWDM)*(Om-Ob),
'Omega_ncdm': fWDM*(Om-Ob),
'N_ur': 3.04,
'N_ncdm': 1,
'm_ncdm': 1000*m_nu,
'T_ncdm': 0.715985,
'n_s': ns,
'A_s': As,
'P_k_max_h/Mpc': 200,
'k_per_decade_for_pk': 5,
'output': 'mPk',
'z_pk': 0.0,
'ncdm_fluid_approximation': 3,
}
else:
CLASSparams = {
'h': h0,
'T_cmb': 2.726,
'Omega_b': Ob,
'Omega_cdm': (Om-Ob),
'N_ur': 3.04,
'n_s': ns,
'A_s': As,
'P_k_max_h/Mpc': 200,
'k_per_decade_for_pk': 5,
'output': 'mPk',
'z_pk': 0.0,
'ncdm_fluid_approximation': 1,
}
CLASScosmo = Class()
CLASScosmo.set(CLASSparams)
CLASScosmo.compute()
s8 = CLASScosmo.sigma8()
k_bin = np.logspace(-3, 2, 1000) #in [h/Mpc]
z_bin = np.array([0])
pk_bin = []
for i in range(len(k_bin)):
pk_bin += [CLASScosmo.pk_lin(k_bin[i]*h0,0)]
pk_bin = np.array(pk_bin)*h0**3 # [Mpc/h]^3
CLASScosmo.struct_cleanup()
CLASScosmo.empty()
return k_bin, pk_bin
def Transfer(k,pk,DMmass,OmegaWDM,h):
"""
initially from Bode, Ostriker & Turok (2001)
but values from Viel 2005 (nu updated)
"""
OmegaWDM = OmegaWDM/0.25
h = h/0.7
nu = 1.12
alpha = 0.049*DMmass**(-1.11) * (OmegaWDM)**(0.11) * (h)**(1.22)
T = (1+(alpha*k)**(2*nu))**(-5/nu)
pk = pk*T**2
return pk
####################################################################
# main
####################################################################
if __name__ == '__main__':
opt = parser.parse_args()
# use an existing power spectrum
if opt.file is not None:
k, pk = np.loadtxt(opt.file[0], unpack=True, usecols = (0,1))
Pk = Transfer(k,pk,opt.DMmass,opt.Om,opt.h0)
# use Classy:
else:
# Cosmological parameters
Om = opt.Om
Ob = opt.Ob
h0 = opt.h0
As = opt.As
ns = opt.ns
# redshift (not used)
z0 = opt.redshift
MDM = opt.DMmass # dark matter mass in keV
fWDM = opt.fWDM # warm dark matter fraction
cosmopars = Om, Ob, As, h0, ns, MDM, fWDM
# calculate power spectrum
k, Pk = power_spectrum(cosmopars)
if opt.outputfile is not None:
body = np.transpose([k, Pk])
np.savetxt(opt.outputfile,body,delimiter='\t')

Event Timeline

c4science · Help