Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88362182
satlib.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
Fri, Oct 18, 09:50
Size
14 KB
Mime Type
text/x-python
Expires
Sun, Oct 20, 09:50 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
21759915
Attached To
rARRAKIHS ARRAKIHS
satlib.py
View Options
import
numpy
as
np
from
scipy
import
special
from
scipy.integrate
import
quad
a
=
1.859e9
# to fit Sawala 2017
b
=
-
1.915
# to fit Sawala 2017
def
CDM_DistributionFunction
(
Mh
):
'''
dN/dM :
'''
dNdM
=
a
*
(
Mh
)
**
b
return
dNdM
def
CDM_CumulativeMass
(
Mh
):
NM
=
-
a
*
(
Mh
)
**
(
b
+
1
)
/
(
b
+
1
)
return
NM
def
WDMCDMRatio
(
Mhs
):
"""
Ratio of the halo mass function of the WDM to CDM (Forouhar 2023)
"""
lnMhs
=
np
.
log10
(
Mhs
)
f
=
np
.
array
([
0
,
0.004
,
0.008
,
0.025
,
0.092
,
0.253
,
0.396
,
0.485
,
0.620
,
0.746
,
0.827
,
0.869
,
0.873
,
0.936
,
1
])
lnMh
=
np
.
array
([
7
,
7.176
,
7.549
,
7.934
,
8.367
,
8.900
,
9.244
,
9.457
,
9.742
,
10.00
,
10.35
,
10.55
,
11.10
,
11.38
,
11.63
])
fs
=
np
.
interp
(
lnMhs
,
lnMh
,
f
)
return
fs
# cumulative mass fct from CDM_CumulativeMass
fctCDM_CumulativeMass
=
np
.
vectorize
(
lambda
x
:
quad
(
lambda
M
:
CDM_DistributionFunction
(
M
)
,
x
,
1e12
)[
0
])
# cumulative mass fct from CDM_CumulativeMass
fctWDM_CumulativeMass
=
np
.
vectorize
(
lambda
x
:
quad
(
lambda
M
:
CDM_DistributionFunction
(
M
)
*
WDMCDMRatio
(
M
)
,
x
,
1e12
)[
0
])
def
GetStellarMassHaloMassRelation
(
model
=
'rj2018'
):
if
model
==
"rj2018"
:
import
pickle
# read Mh vs Lv
f
=
open
(
"MsvsMh_2018.pkl"
,
"rb"
)
Mhbins
=
pickle
.
load
(
f
,
encoding
=
"latin1"
)
Lvmins
=
pickle
.
load
(
f
,
encoding
=
"latin1"
)
Lvmaxs
=
pickle
.
load
(
f
,
encoding
=
"latin1"
)
f
.
close
()
return
Mhbins
,
Lvmins
,
Lvmaxs
elif
model
==
"model1"
:
Mhbins
=
np
.
linspace
(
7
,
11
,
20
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_low_1.txt"
)
Mh1
=
data
[:,
0
]
Msl
=
data
[:,
1
]
Lvmins
=
np
.
interp
(
Mhbins
,
Mh1
,
Msl
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_high_1.txt"
)
Mh1
=
data
[:,
0
]
Msu
=
data
[:,
1
]
Lvmaxs
=
np
.
interp
(
Mhbins
,
Mh1
,
Msu
)
return
10
**
Mhbins
,
10
**
Lvmins
,
10
**
Lvmaxs
elif
model
==
"model2"
:
Mhbins
=
np
.
linspace
(
7
,
11
,
20
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_low_2.txt"
)
Mh1
=
data
[:,
0
]
Msl
=
data
[:,
1
]
Lvmins
=
np
.
interp
(
Mhbins
,
Mh1
,
Msl
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_high_2.txt"
)
Mh1
=
data
[:,
0
]
Msu
=
data
[:,
1
]
Lvmaxs
=
np
.
interp
(
Mhbins
,
Mh1
,
Msu
)
return
10
**
Mhbins
,
10
**
Lvmins
,
10
**
Lvmaxs
elif
model
==
"model3"
:
Mhbins
=
np
.
linspace
(
7
,
11
,
20
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_low_3.txt"
)
Mh1
=
data
[:,
0
]
Msl
=
data
[:,
1
]
Lvmins
=
np
.
interp
(
Mhbins
,
Mh1
,
Msl
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_high_3.txt"
)
Mh1
=
data
[:,
0
]
Msu
=
data
[:,
1
]
Lvmaxs
=
np
.
interp
(
Mhbins
,
Mh1
,
Msu
)
return
10
**
Mhbins
,
10
**
Lvmins
,
10
**
Lvmaxs
elif
model
==
"model4"
:
Mhbins
=
np
.
linspace
(
7
,
11
,
20
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_low_4.txt"
)
Mh1
=
data
[:,
0
]
Msl
=
data
[:,
1
]
Lvmins
=
np
.
interp
(
Mhbins
,
Mh1
,
Msl
)
data
=
np
.
loadtxt
(
"sales2022/MsMh_high_4.txt"
)
Mh1
=
data
[:,
0
]
Msu
=
data
[:,
1
]
Lvmaxs
=
np
.
interp
(
Mhbins
,
Mh1
,
Msu
)
return
10
**
Mhbins
,
10
**
Lvmins
,
10
**
Lvmaxs
class
DM_Model
():
def
__init__
(
self
,
Mmin
,
Mmax
,
N
):
self
.
Mmin
=
Mmin
# minimal halo mass
self
.
Mmax
=
Mmax
# maximal halo mass
self
.
N
=
N
# number of mass bins
# halo mass bins
self
.
Mh
=
10
**
np
.
linspace
(
np
.
log10
(
self
.
Mmin
),
np
.
log10
(
self
.
Mmax
),
self
.
N
)
def
DistributionFunction
(
self
,
Mh
):
'''
return the dark matter distribution function
i.e., dN/dM
'''
pass
def
CumulativeMass
(
self
,
Mh
):
'''
return the cumulative halo mass
'''
pass
def
NumberOfHaloes
(
self
):
'''
return the number of haloes
'''
return
int
(
self
.
CumulativeMass
(
self
.
Mmin
))
##############################################
# CDM from Sawala 2017 fit
##############################################
class
CDM_Sawala_Model
(
DM_Model
):
def
__init__
(
self
,
Mmin
,
Mmax
,
N
):
# do the initialisation
super
()
.
__init__
(
Mmin
,
Mmax
,
N
)
self
.
a
=
1.859e9
# to fit Sawala 2017
self
.
b
=
-
1.915
# to fit Sawala 2017
# cumulative number of haloes
NM
=
self
.
CumulativeMass
()
# compute the sampling function
self
.
fctSample
=
np
.
vectorize
(
lambda
x
:
10
**
np
.
interp
(
x
,
NM
[::
-
1
]
/
NM
.
max
(),
np
.
log10
(
self
.
Mh
[::
-
1
]))
)
def
DistributionFunction
(
self
,
Mh
=
None
):
'''
return the dark matter distribution function
i.e., dN/dM
'''
if
Mh
is
None
:
Mh
=
self
.
Mh
dNdM
=
self
.
a
*
(
Mh
)
**
self
.
b
return
dNdM
def
CumulativeMass
(
self
,
Mh
=
None
):
'''
return the cumulative halo mass
'''
if
Mh
is
None
:
Mh
=
self
.
Mh
NM
=
-
self
.
a
*
(
Mh
)
**
(
self
.
b
+
1
)
/
(
self
.
b
+
1
)
return
NM
##############################################
# WDM from Sawala 2017 + Forouhar 2023
##############################################
class
CDM_Forouhar_Model
(
DM_Model
):
def
__init__
(
self
,
Mmin
,
Mmax
,
N
):
# do the initialisation
super
()
.
__init__
(
Mmin
,
Mmax
,
N
)
self
.
CDM
=
CDM_Sawala_Model
(
Mmin
,
Mmax
,
N
)
# cumulative mass fct from CDM_CumulativeMass
self
.
fctCumulativeMass
=
np
.
vectorize
(
lambda
x
:
quad
(
lambda
M
:
self
.
DistributionFunction
(
M
),
x
,
1e12
)[
0
])
# cumulative number of haloes
NM
=
self
.
CumulativeMass
()
# compute the sampling function
self
.
fctSample
=
np
.
vectorize
(
lambda
x
:
10
**
np
.
interp
(
x
,
NM
[::
-
1
]
/
NM
.
max
(),
np
.
log10
(
self
.
Mh
[::
-
1
]))
)
def
DistributionFunction
(
self
,
Mh
=
None
):
'''
return the dark matter distribution function
i.e., dN/dM
'''
if
Mh
is
None
:
Mh
=
self
.
Mh
dNdM
=
self
.
CDM
.
DistributionFunction
(
Mh
)
*
self
.
WDMCDMRatio
(
Mh
)
return
dNdM
def
CumulativeMass
(
self
,
Mh
=
None
):
'''
return the cumulative halo mass
'''
if
Mh
is
None
:
Mh
=
self
.
Mh
return
self
.
fctCumulativeMass
(
Mh
)
def
WDMCDMRatio
(
self
,
Mh
=
None
):
"""
Ratio of the halo mass function of the WDM to CDM (Forouhar 2023)
"""
if
Mh
is
None
:
Mh
=
self
.
Mh
lnMhs
=
np
.
log10
(
Mh
)
f
=
np
.
array
([
0
,
0.004
,
0.008
,
0.025
,
0.092
,
0.253
,
0.396
,
0.485
,
0.620
,
0.746
,
0.827
,
0.869
,
0.873
,
0.936
,
1
])
lnMh
=
np
.
array
([
7
,
7.176
,
7.549
,
7.934
,
8.367
,
8.900
,
9.244
,
9.457
,
9.742
,
10.00
,
10.35
,
10.55
,
11.10
,
11.38
,
11.63
])
fs
=
np
.
interp
(
lnMhs
,
lnMh
,
f
)
return
fs
##############################################
# WDM from Bode, Ostriker & Turok (2001)
# Viel 2005 (nu updated)
##############################################
# Note : here we avoid to use classy to get
# the power spectrum
from
scipy.integrate
import
cumtrapz
,
trapz
import
genmassfct
as
gm
class
GEN_Model
(
DM_Model
):
def
__init__
(
self
,
Mmin
,
Mmax
,
N
,
MDM
=
2
,
fWDM
=
0
,
M0
=
1e12
,
Om
=
0.315
,
h0
=
0.674
):
# do the initialisation
super
()
.
__init__
(
Mmin
,
Mmax
,
N
)
# unperturbed power spectrum filename
filename
=
'./CDM_PS.dat'
# temporary file
self
.
filename
=
'/tmp/powerspectrum_pk.dat'
k
,
pk
=
np
.
loadtxt
(
filename
,
unpack
=
True
,
usecols
=
(
0
,
1
))
if
MDM
is
not
None
:
Pk
=
self
.
Transfer
(
k
,
pk
,
MDM
,
Om
,
h0
)
else
:
Pk
=
pk
# save (needed for )
body
=
np
.
transpose
([
k
,
Pk
])
np
.
savetxt
(
self
.
filename
,
body
,
delimiter
=
'
\t
'
)
# initialise parameters
par
=
gmf
.
par
()
par
.
window
.
window
=
"sharpk"
par
.
code
.
rmin
=
0.008
par
.
mf
.
q
=
1.0
par
.
mf
.
p
=
0.3
par
.
mf
.
c
=
2.5
# redshift
z0
=
0
# calculate stellite function
par
.
file
.
psfct
=
self
.
filename
self
.
mbin
,
dNsatdlnM
=
gmf
.
dNsatdlnm
(
M0
,
z0
,
par
)
self
.
Nsat
=
Nsat_integral
(
self
.
mbin
,
dNsatdlnM
)
# cumulative number of haloes
NM
=
self
.
CumulativeMass
()
# compute the sampling function
self
.
fctSample
=
np
.
vectorize
(
lambda
x
:
10
**
np
.
interp
(
x
,
NM
[::
-
1
]
/
NM
.
max
(),
np
.
log10
(
self
.
Mh
[::
-
1
]))
)
def
CumulativeMass
(
self
,
Mh
=
None
):
'''
return the cumulative halo mass
'''
if
Mh
is
None
:
Mh
=
self
.
Mh
return
np
.
interp
(
Mh
,
self
.
mbin
,
self
.
Nsat
)
def
Transfer
(
self
,
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
##############################################
# Generic model (Schneider 2014)
##############################################
# Note : deprecated
from
scipy.integrate
import
cumtrapz
,
trapz
from
classy
import
Class
import
genmassfct
as
gmf
def
ps_class
(
cosmopars
):
"""
Calculate power spectrum with 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'
:
10
,
'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'
:
10
,
'output'
:
'mPk'
,
'z_pk'
:
0.0
,
'ncdm_fluid_approximation'
:
1
,
}
CLASScosmo
=
Class
()
CLASScosmo
.
set
(
CLASSparams
)
CLASScosmo
.
compute
()
s8
=
CLASScosmo
.
sigma8
()
print
(
"s8 = "
,
s8
)
k_bin
=
np
.
logspace
(
np
.
log10
(
1e-3
),
np
.
log10
(
200
),
500
)
#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
Nsat_integral
(
mbin
,
dNsatdlnm
):
"""
Cumulative mass function Nsat(>M)
"""
Nsat
=
trapz
(
dNsatdlnm
/
mbin
,
mbin
)
-
cumtrapz
(
dNsatdlnm
/
mbin
,
mbin
,
initial
=
dNsatdlnm
[
0
]
/
mbin
[
0
])
return
Nsat
class
oldGEN_Model
(
DM_Model
):
def
__init__
(
self
,
Mmin
,
Mmax
,
N
,
MDM
=
2
,
fWDM
=
0
,
M0
=
1e12
):
# do the initialisation
super
()
.
__init__
(
Mmin
,
Mmax
,
N
)
# temporary file
self
.
filename
=
'/tmp/powerspectrum_pk.dat'
# Cosmologycal parameters
Om
=
0.315
Ob
=
0.048
h0
=
0.681
As
=
2.07e-9
ns
=
0.963
# redshift
z0
=
0
self
.
MDM
=
MDM
# dark matter mass in keV
self
.
fWDM
=
fWDM
# warm dark matter fraction
self
.
M0
=
M0
# Host halo mass (Msun/h)
cosmopars
=
Om
,
Ob
,
As
,
h0
,
ns
,
self
.
MDM
,
fWDM
# calculate power spectrum
kbin
,
PSbin
=
ps_class
(
cosmopars
)
body
=
np
.
transpose
([
kbin
,
PSbin
])
np
.
savetxt
(
self
.
filename
,
body
,
delimiter
=
'
\t
'
)
# initialise parameters
par
=
gmf
.
par
()
par
.
window
.
window
=
"sharpk"
par
.
code
.
rmin
=
0.008
par
.
mf
.
q
=
1.0
par
.
mf
.
p
=
0.3
par
.
mf
.
c
=
2.5
# calculate stellite function
par
.
file
.
psfct
=
self
.
filename
self
.
mbin
,
dNsatdlnM
=
gmf
.
dNsatdlnm
(
M0
,
z0
,
par
)
self
.
Nsat
=
Nsat_integral
(
self
.
mbin
,
dNsatdlnM
)
# cumulative number of haloes
NM
=
self
.
CumulativeMass
()
# compute the sampling function
self
.
fctSample
=
np
.
vectorize
(
lambda
x
:
10
**
np
.
interp
(
x
,
NM
[::
-
1
]
/
NM
.
max
(),
np
.
log10
(
self
.
Mh
[::
-
1
]))
)
def
CumulativeMass
(
self
,
Mh
=
None
):
'''
return the cumulative halo mass
'''
if
Mh
is
None
:
Mh
=
self
.
Mh
return
np
.
interp
(
Mh
,
self
.
mbin
,
self
.
Nsat
)
def
Msol2Mag
(
mass
,
Age
=
12
,
MH
=-
2
,
filter
=
"F475X"
):
'''
mass : in Msol
Age : in Gyr
MH : in [M/H]
'''
import
stars_class
mass
=
np
.
array
([
mass
])
Age
=
np
.
array
([
Age
])
MH
=
np
.
array
([
MH
])
###################################
# compute Magnitudes
###################################
# get the number of stars in each mass bin
Nstars
=
stars_class
.
Stars_fun
(
mass
,
None
,
None
,
'normed_3slope'
)
##############################
# F475X magnitude
if
filter
==
"F475X"
:
M
=
stars_class
.
HST475X_fun
(
None
,
Age
,
MH
)
# convert to flux (ignore the zero point)
F
=
10
**
(
-
M
/
2.5
)
# sum the contribution of the mass bins
F
=
np
.
sum
(
F
*
Nstars
,
axis
=
0
)
# compute the absolute magnitude in each pixel (as before we ignore the zero point)
M
=
-
2.5
*
np
.
log10
(
F
)
return
M
[
0
]
##############################
# VIS Euclid magnitude
if
filter
==
"VISeuclid"
:
M
=
stars_class
.
VISeuclid_fun
(
None
,
Age
,
MH
)
# convert to flux (ignore the zero point)
F
=
10
**
(
-
M
/
2.5
)
# sum the contribution of the mass bins
F
=
np
.
sum
(
F
*
Nstars
,
axis
=
0
)
# compute the absolute magnitude in each pixel (as before we ignore the zero point)
M
=
-
2.5
*
np
.
log10
(
F
)
return
M
[
0
]
##############################
# Y Euclid magnitude
if
filter
==
"Yeuclid"
:
M
=
stars_class
.
Yeuclid_fun
(
None
,
Age
,
MH
)
# convert to flux (ignore the zero point)
F
=
10
**
(
-
M
/
2.5
)
# sum the contribution of the mass bins
F
=
np
.
sum
(
F
*
Nstars
,
axis
=
0
)
# compute the absolute magnitude in each pixel (as before we ignore the zero point)
M
=
-
2.5
*
np
.
log10
(
F
)
return
M
[
0
]
##############################
# J Euclid magnitude
if
filter
==
"Jeuclid"
:
M
=
stars_class
.
Jeuclid_fun
(
None
,
Age
,
MH
)
# convert to flux (ignore the zero point)
F
=
10
**
(
-
M
/
2.5
)
# sum the contribution of the mass bins
F
=
np
.
sum
(
F
*
Nstars
,
axis
=
0
)
# compute the absolute magnitude in each pixel (as before we ignore the zero point)
M
=
-
2.5
*
np
.
log10
(
F
)
return
M
[
0
]
##############################
# Vasdekhis
if
filter
==
"V"
:
from
pNbody.SSP
import
libvazdekis
LObj
=
libvazdekis
.
VazdekisLuminosities
(
"/home/revaz/.pNbody/opt/SSP/vazdekis_kb_mu1.3.txt"
,
band
=
"V"
)
LObj
.
ExtrapolateMatrix
(
order
=
1
,
s
=
0
)
LObj
.
CreateInterpolator
()
LObj
.
Extrapolate2DMatrix
()
# luminosity
Lv
=
LObj
.
Luminosities
(
MH
,
Age
)
*
mass
# Magnitude
Mag_Vega
=
4.81
Mv
=
Mag_Vega
-
2.5
*
np
.
log10
(
Lv
)
return
Mv
Event Timeline
Log In to Comment