Homec4science
Diffusion ARRAKIHS (master) Browse

arrakihs/Workspace/CumNumOfSat/sto.satComputeNCumvsLvWithErrors.py
0d50a2737ed2master

Last Change
Raw File

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)
c4science · Help