Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80560112
cosmo_generate_power_spectum.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Sep 1, 09:39
Size
5 KB
Mime Type
text/x-python
Expires
Tue, Sep 3, 09:39 (2 d)
Engine
blob
Format
Raw Data
Handle
20359885
Attached To
rARRAKIHS ARRAKIHS
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
Log In to Comment