Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F83480350
sto.satComputeNCumvsLvWithErrors.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
Tue, Sep 17, 09:14
Size
10 KB
Mime Type
text/x-python
Expires
Thu, Sep 19, 09:14 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20843763
Attached To
rARRAKIHS ARRAKIHS
sto.satComputeNCumvsLvWithErrors.py
View Options
#!/usr/bin/env python3
import
sys
,
os
,
string
import
argparse
import
satlib
import
numpy
as
np
import
Ptools
as
pt
import
pickle
from
scipy.interpolate
import
splrep
,
splev
from
tqdm
import
tqdm
####################################################################
# option parser
####################################################################
description
=
""
epilog
=
""""""
parser
=
argparse
.
ArgumentParser
(
description
=
description
,
epilog
=
epilog
,
formatter_class
=
argparse
.
RawDescriptionHelpFormatter
)
parser
.
add_argument
(
"--Mmin"
,
action
=
"store"
,
dest
=
"Mmin"
,
metavar
=
'FLOAT'
,
type
=
float
,
default
=
5e7
,
help
=
'the min halo mass'
)
parser
.
add_argument
(
"--Mmax"
,
action
=
"store"
,
dest
=
"Mmax"
,
metavar
=
'FLOAT'
,
type
=
float
,
default
=
1e11
,
help
=
'the max halo mass'
)
parser
.
add_argument
(
"-N"
,
action
=
"store"
,
dest
=
"N"
,
metavar
=
'INT'
,
type
=
int
,
default
=
1000
,
help
=
'number of bins'
)
parser
.
add_argument
(
"--Ngal"
,
action
=
"store"
,
dest
=
"Ngal"
,
metavar
=
'INT'
,
type
=
int
,
default
=
100
,
help
=
'number of galaxies'
)
parser
.
add_argument
(
"--Nrealisations"
,
action
=
"store"
,
dest
=
"Nrealisations"
,
metavar
=
'INT'
,
type
=
int
,
default
=
100
,
help
=
'number of realizations'
)
parser
.
add_argument
(
"--LvvsMh_model"
,
action
=
"store"
,
type
=
str
,
dest
=
"LvvsMh_model"
,
default
=
'rj2018'
,
help
=
"Lv Mh model"
)
parser
.
add_argument
(
"-p"
,
action
=
"store"
,
dest
=
"ps"
,
metavar
=
'FILE NAME'
,
type
=
str
,
default
=
None
,
help
=
'output file name'
)
parser
.
add_argument
(
"--dpi"
,
action
=
"store"
,
dest
=
"dpi"
,
type
=
float
,
default
=
300
,
help
=
"DPI of the saved file"
,
metavar
=
" FLOAT"
)
parser
.
add_argument
(
"--data"
,
action
=
"store_true"
,
dest
=
"data"
,
default
=
False
,
help
=
'add data'
)
####################################################################
# main
####################################################################
def
MakePlot
(
opt
):
# get Mh vs Lv relation
Mhbins
,
Lvmins
,
Lvmaxs
=
satlib
.
GetStellarMassHaloMassRelation
(
model
=
opt
.
LvvsMh_model
)
# halo bins
lnMh
=
np
.
linspace
(
np
.
log10
(
opt
.
Mmin
),
np
.
log10
(
opt
.
Mmax
),
opt
.
N
)
Mh0
=
10
**
lnMh
# cumulative number of haloes
NM_CDM
=
satlib
.
CDM_CumulativeMass
(
Mh0
)
NM_WDM
=
satlib
.
fctWDM_CumulativeMass
(
Mh0
)
# define the CDM sampling function
fct_CDM_Sample
=
np
.
vectorize
(
lambda
x
:
10
**
np
.
interp
(
x
,
NM_CDM
[::
-
1
]
/
NM_CDM
.
max
(),
np
.
log10
(
Mh0
[::
-
1
]))
)
maxCDM
=
(
NM_CDM
/
NM_CDM
.
max
())
.
min
()
# define the WDM sampling function
fct_WDM_Sample
=
np
.
vectorize
(
lambda
x
:
10
**
np
.
interp
(
x
,
NM_WDM
[::
-
1
]
/
NM_WDM
.
max
(),
np
.
log10
(
Mh0
[::
-
1
]))
)
maxWDM
=
(
NM_WDM
/
NM_WDM
.
max
())
.
min
()
Nh_CDM
=
int
(
satlib
.
CDM_CumulativeMass
(
opt
.
Mmin
))
Nh_WDM
=
int
(
satlib
.
fctWDM_CumulativeMass
(
opt
.
Mmin
))
Nbins
=
100
Nsats_CDM
=
np
.
zeros
((
opt
.
Ngal
,
Nbins
-
1
))
Nsats_WDM
=
np
.
zeros
((
opt
.
Ngal
,
Nbins
-
1
))
Ns_CDM_realisations
=
np
.
zeros
((
opt
.
Nrealisations
,
Nbins
-
1
))
Ns_WDM_realisations
=
np
.
zeros
((
opt
.
Nrealisations
,
Nbins
-
1
))
for
realisation
in
tqdm
(
range
(
opt
.
Nrealisations
)):
for
j
in
range
(
opt
.
Ngal
):
#print(realisation,j)
#################################################
# generate N haloes for the CDM
#################################################
Mhs
=
fct_CDM_Sample
(
np
.
random
.
random
(
Nh_CDM
))
###################
# loop over haloes
logLv
=
np
.
zeros
(
len
(
Mhs
))
for
i
,
Mh
in
enumerate
(
Mhs
):
if
Mh
<
1e8
:
continue
# find the nearest values in Mbins
# and set a Luminosity
idx
=
np
.
argmin
(
np
.
fabs
(
Mh
-
Mhbins
)
)
logLvmins
=
np
.
log10
(
Lvmins
[
idx
])
logLvmaxs
=
np
.
log10
(
Lvmaxs
[
idx
])
logLv
[
i
]
=
np
.
random
.
uniform
(
logLvmins
,
logLvmaxs
)
######################
# cumulative
bins
=
np
.
linspace
(
2
,
12
,
Nbins
)
Ns
,
Lv
=
np
.
histogram
(
logLv
,
bins
)
# now, need to cumulate
# invert vector first and reinvert after
Nsat
=
np
.
add
.
accumulate
(
Ns
[
-
np
.
arange
(
1
,
len
(
Ns
)
+
1
)])
Nsat
=
Nsat
[
-
np
.
arange
(
1
,
len
(
Nsat
)
+
1
)]
Lv
=
Lv
[
1
:]
# add
Nsats_CDM
[
j
]
=
Nsat
#################################################
# generate N haloes for the WDM
#################################################
Mhs
=
fct_WDM_Sample
(
np
.
random
.
random
(
Nh_WDM
))
###################
# loop over haloes
logLv
=
np
.
zeros
(
len
(
Mhs
))
for
i
,
Mh
in
enumerate
(
Mhs
):
if
Mh
<
1e8
:
continue
# find the nearest values in Mbins
# and set a Luminosity
idx
=
np
.
argmin
(
np
.
fabs
(
Mh
-
Mhbins
)
)
logLvmins
=
np
.
log10
(
Lvmins
[
idx
])
logLvmaxs
=
np
.
log10
(
Lvmaxs
[
idx
])
logLv
[
i
]
=
np
.
random
.
uniform
(
logLvmins
,
logLvmaxs
)
######################
# cumulative
bins
=
np
.
linspace
(
2
,
12
,
Nbins
)
Ns
,
Lv
=
np
.
histogram
(
logLv
,
bins
)
# now, need to cumulate
# invert vector first and reinvert after
Nsat
=
np
.
add
.
accumulate
(
Ns
[
-
np
.
arange
(
1
,
len
(
Ns
)
+
1
)])
Nsat
=
Nsat
[
-
np
.
arange
(
1
,
len
(
Nsat
)
+
1
)]
Lv
=
Lv
[
1
:]
# add
Nsats_WDM
[
j
]
=
Nsat
# plot individual galaxies
#pt.plot(10**Lv,Nsats_CDM[j],'blue',alpha=0.5)
#pt.plot(10**Lv,Nsats_WDM[j],'red' ,alpha=0.5)
Ns_CDM_mean
=
Nsats_CDM
.
mean
(
axis
=
0
)
Ns_WDM_mean
=
Nsats_WDM
.
mean
(
axis
=
0
)
Ns_CDM_realisations
[
realisation
]
=
Ns_CDM_mean
Ns_WDM_realisations
[
realisation
]
=
Ns_WDM_mean
# plot means
#pt.plot(10**Lv,Ns_CDM_mean,lw=1,alpha=0.5, c='blue')
#pt.plot(10**Lv,Ns_WDM_mean,lw=1,alpha=0.5, c='red')
Ns_CDM_realisations_mean
=
Ns_CDM_realisations
.
mean
(
axis
=
0
)
Ns_CDM_realisations_std
=
Ns_CDM_realisations
.
std
(
axis
=
0
)
Ns_CDM_meanp
=
Ns_CDM_realisations_mean
+
3
*
Ns_CDM_realisations_std
Ns_CDM_meanm
=
Ns_CDM_realisations_mean
-
3
*
Ns_CDM_realisations_std
Ns_WDM_realisations_mean
=
Ns_WDM_realisations
.
mean
(
axis
=
0
)
Ns_WDM_realisations_std
=
Ns_WDM_realisations
.
std
(
axis
=
0
)
Ns_WDM_meanp
=
Ns_WDM_realisations_mean
+
3
*
Ns_WDM_realisations_std
Ns_WDM_meanm
=
Ns_WDM_realisations_mean
-
3
*
Ns_WDM_realisations_std
ax
=
pt
.
gca
()
ax
.
fill_between
(
10
**
Lv
,
Ns_CDM_meanp
,
Ns_CDM_meanm
,
facecolor
=
'blue'
,
alpha
=
0.5
,
label
=
r"$\textrm{CDM}$"
)
ax
.
fill_between
(
10
**
Lv
,
Ns_WDM_meanp
,
Ns_WDM_meanm
,
facecolor
=
'red'
,
alpha
=
0.5
,
label
=
r"$\textrm{WDM}$"
)
if
opt
.
data
:
################################################################################
# add data from Nadler, Drlica etc.
ax
=
pt
.
gca
()
s
=
0.0
k
=
1
def
read_data
(
f
):
data
=
np
.
loadtxt
(
f
,
delimiter
=
','
)
Mv
=
data
[:,
0
]
N
=
data
[:,
1
]
Lv
=
pow
(
10
,
(
4.74
-
Mv
)
/
2.5
)
return
Lv
,
N
###################
# Nadler
Lvp
,
Np
=
read_data
(
"data/Nadler_rawt1.txt"
)
Lvm
,
Nm
=
read_data
(
"data/Nadler_rawb1.txt"
)
ap
=
splrep
(
Lvp
,
Np
,
s
=
s
,
k
=
k
)
am
=
splrep
(
Lvm
,
Nm
,
s
=
s
,
k
=
k
)
Np
=
splev
(
Lvp
,
ap
)
Nm
=
splev
(
Lvp
,
am
)
#pt.plot(Lvp,Np,'k-')
#pt.plot(Lvm,Nm,'k-')
ax
.
fill_between
(
Lvp
,
Np
,
Nm
,
facecolor
=
'green'
,
alpha
=
0.05
,
label
=
r"$\textrm{Nadler\,\,et\,\,al.\,\,2018a}$"
)
Lvp
,
Np
=
read_data
(
"data/Nadler_rawt2.txt"
)
Lvm
,
Nm
=
read_data
(
"data/Nadler_rawb2.txt"
)
ap
=
splrep
(
Lvp
,
Np
,
s
=
s
,
k
=
k
)
am
=
splrep
(
Lvm
,
Nm
,
s
=
s
,
k
=
k
)
Np
=
splev
(
Lvp
,
ap
)
Nm
=
splev
(
Lvp
,
am
)
#pt.plot(Lvp,Np,'k-')
#pt.plot(Lvm,Nm,'k-')
ax
.
fill_between
(
Lvp
,
Np
,
Nm
,
facecolor
=
'green'
,
alpha
=
0.05
)
###################
# Drilca
Lv
,
N
=
read_data
(
"data/Drlica_raw_weighted.txt"
)
pt
.
plot
(
Lv
,
N
,
'k-'
,
label
=
r"$\textrm{DES+PS1,\,\,weighted\,\,(Drlica-Wagner\,\,et\,\,al.\,\,2020)}$"
,
alpha
=
0.5
)
Lv
,
N
=
read_data
(
"data/Drlica_raw_detected.txt"
)
pt
.
plot
(
Lv
,
N
,
'k--'
,
label
=
r"$\textrm{DES+PS1,\,\,detected\,\,(Drlica-Wagner\,\,et\,\,al.\,\,2020)}$"
,
alpha
=
0.5
)
Lv
,
N
=
read_data
(
"data/Drlica_raw_all.txt"
)
pt
.
plot
(
Lv
,
N
,
'r-'
,
label
=
r"$\textrm{All\,\,known\,\,satellites\,\,(Drlica-Wagner\,\,et\,\,al.\,\,2020)}$"
,
alpha
=
0.5
)
###########################
# finalize
###########################
xmin
=
5e2
xmax
=
5e9
ymin
=
1
ymax
=
500
xlabel
=
r"$\rm{V-band\,\,Luminosity}$"
ylabel
=
r"$\rm{Cumulative\,\,number\,\,of\,\,galaxies}\,\,(D<300\,\rm{kpc})$"
pt
.
SetAxis
(
xmin
,
xmax
,
ymin
,
ymax
,
log
=
"xy"
)
pt
.
xlabel
(
xlabel
,
fontsize
=
pt
.
labelfont
)
pt
.
ylabel
(
ylabel
,
fontsize
=
pt
.
labelfont
)
pt
.
grid
(
False
)
pt
.
legend
()
pt
.
title
(
r"$\rm{
%d
\,\,galaxies\,\,observed\,\,+\,\,model=
%s
}$"
%
(
opt
.
Ngal
,
opt
.
LvvsMh_model
),
fontsize
=
pt
.
labelfont
)
pt
.
plot
([
1e5
,
1e5
],[
0.1
,
1000
],
c
=
'k'
,
ls
=
':'
)
if
__name__
==
'__main__'
:
opt
=
parser
.
parse_args
()
files
=
[]
pt
.
InitPlot
(
files
,
opt
)
pt
.
labelfont
=
20
# pt.figure(figsize=(8*2,6*2))
# pt.figure(dpi=10)
#fig = pt.gcf()
# fig.subplots_adjust(left=0.1)
# fig.subplots_adjust(right=1)
# fig.subplots_adjust(bottom=0.12)
# fig.subplots_adjust(top=0.95)
# fig.subplots_adjust(wspace=0.25)
# fig.subplots_adjust(hspace=0.02)
MakePlot
(
opt
)
pt
.
EndPlot
(
files
,
opt
)
Event Timeline
Log In to Comment