Homec4science
Diffusion ARRAKIHS (master) Browse

arrakihs/Workspace/CumNumOfSat/satComputeNCumvsLv.py
0d50a2737ed2master

Last Change
Raw File

satComputeNCumvsLv.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
####################################################################
# 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=1e10,
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("--LvvsMh_file",
action="store",
type=str,
dest="LvvsMh_file",
default = None,
help="Lv Mh file")
####################################################################
# main
####################################################################
opt = parser.parse_args()
myblue = (0,128/255.,255/255.)
myred = (255/255.,137/255.,137/255.)
# read Mh vs Lv
f = open(opt.LvvsMh_file,"rb")
Mhbins = pickle.load(f,encoding="latin1")
Lvmins = pickle.load(f,encoding="latin1")
Lvmaxs = pickle.load(f,encoding="latin1")
f.close()
# 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()
Nbins = 100
Nsas_CDM = np.zeros((opt.Ngal,Nbins-1))
Nsas_WDM = np.zeros((opt.Ngal,Nbins-1))
for j in range(opt.Ngal):
print(j)
#################################################
# generate N haloes for the CDM
#################################################
Nh = int(satlib.CDM_CumulativeMass(opt.Mmin))
Mhs = fct_CDM_Sample(np.random.random(Nh))
###################
# 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,10,Nbins)
Ns,Lv = np.histogram(logLv,bins)
# now, need to cumulate
# invert vector first and reinvert after
Nsa = np.add.accumulate(Ns[-np.arange(1,len(Ns)+1)])
Nsa = Nsa[-np.arange(1,len(Nsa)+1)]
Lv = Lv[1:]
# add
Nsas_CDM[j] = Nsa
#################################################
# generate N haloes for the WDM
#################################################
Nh = int(satlib.fctWDM_CumulativeMass(opt.Mmin))
Mhs = fct_WDM_Sample(np.random.random(Nh))
###################
# 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,10,Nbins)
Ns,Lv = np.histogram(logLv,bins)
# now, need to cumulate
# invert vector first and reinvert after
Nsa = np.add.accumulate(Ns[-np.arange(1,len(Ns)+1)])
Nsa = Nsa[-np.arange(1,len(Nsa)+1)]
Lv = Lv[1:]
# add
Nsas_WDM[j] = Nsa
Ns_CDM_mean = Nsas_CDM.mean(axis=0)
Ns_CDM_std = Nsas_CDM.std(axis=0)
Ns_CDM_meanp = Ns_CDM_mean + Ns_CDM_std
Ns_CDM_meanm = Ns_CDM_mean - Ns_CDM_std
Ns_WDM_mean = Nsas_WDM.mean(axis=0)
Ns_WDM_std = Nsas_WDM.std(axis=0)
Ns_WDM_meanp = Ns_WDM_mean + Ns_WDM_std
Ns_WDM_meanm = Ns_WDM_mean - Ns_WDM_std
pt.plot(10**Lv,Ns_CDM_mean,label=r"$\textrm{CDM}$",lw=5)
pt.plot(10**Lv,Ns_WDM_mean,label=r"$\textrm{WDM}$",lw=5)
#ax = pt.gca()
#ax.fill_between(10**Lv,Ns_CDM_meanp,Ns_CDM_meanm,facecolor='k',alpha=0.1,label=r"$\textrm{CDM}$")
#ax.fill_between(10**Lv,Ns_WDM_meanp,Ns_WDM_meanm,facecolor='k',alpha=0.1,label=r"$\textrm{WDM}$")
data = True
if 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='blue',alpha=0.1,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='blue',alpha=0.1)
###################
# 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)}$")
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)}$")
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)}$")
###########################
# 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)
pt.ylabel(ylabel)
pt.grid(False)
pt.legend()
pt.show()
c4science · Help