Homec4science
Diffusion ARRAKIHS (master) Browse

arrakihs/Workspace/CumNumOfSat/satComputeNCumWithErrors.py
0d50a2737ed2master

Last Change
Raw File

satComputeNCumWithErrors.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("--M0",
action="store",
dest="M0",
metavar='FLOAT',
type=float,
default=1e12,
help='host halo mass')
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("--DMmass",
action="store",
dest="DMmass",
metavar='FLOAT',
type=float,
default=2.0,
help='dark matter particle 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("-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")
####################################################################
# main
####################################################################
def MakePlot(opt):
# halo bins
lnMh = np.linspace(np.log10(opt.Mmin),np.log10(opt.Mmax),opt.N)
Mh0 = 10**lnMh
# define the DM models
#CDM = satlib.CDM_Sawala_Model(opt.Mmin,opt.Mmax,opt.N)
#WDM = satlib.CDM_Forouhar_Model(opt.Mmin,opt.Mmax,opt.N)
CDM = satlib.GEN_Model(opt.Mmin,opt.Mmax,opt.N,MDM=None,fWDM=0,M0=opt.M0)
WDM = satlib.GEN_Model(opt.Mmin,opt.Mmax,opt.N,MDM=opt.DMmass,fWDM=1,M0=opt.M0)
# cumulative number of haloes
NM_CDM = CDM.CumulativeMass()
NM_WDM = WDM.CumulativeMass()
pt.plot(Mh0,NM_WDM,c='k',ls='--',label=r"$\rm{theory}$")
pt.plot(Mh0,NM_CDM,c='k',ls='--')
# number of haloes for the CDM
Nh_CDM = CDM.NumberOfHaloes()
# number of haloes for the WDM
Nh_WDM = WDM.NumberOfHaloes()
Nb = 100
Ncs_CDM = np.zeros((opt.Ngal,Nb-1))
Ncs_WDM = np.zeros((opt.Ngal,Nb-1))
Nc_CDM_realisations = np.zeros((opt.Nrealisations,Nb-1))
Nc_WDM_realisations = np.zeros((opt.Nrealisations,Nb-1))
for realisation in tqdm(range(opt.Nrealisations)):
for j in range(opt.Ngal):
# generate N haloes for the CDM
Mhs = CDM.fctSample(np.random.random(Nh_CDM))
lnm = np.log10(Mhs)
# do the histogram
logMh = np.linspace(np.log10(opt.Mmin),np.log10(opt.Mmax),Nb)
n, bins = np.histogram(lnm, logMh)
Ms = 10**bins
dM = Ms[1:]-Ms[:-1]
dNdM = n/dM
Nc_CDM = np.add.accumulate(np.flip(n,0))
Nc_CDM = np.flip(Nc_CDM,0)
Ncs_CDM[j] = Nc_CDM
# generate N haloes for the WDM
#Mhs = fct_WDM_Sample(np.random.random(Nh_WDM))
Mhs = WDM.fctSample(np.random.random(Nh_WDM))
lnm = np.log10(Mhs)
# do the histogram
logMh = np.linspace(np.log10(opt.Mmin),np.log10(opt.Mmax),Nb)
n, bins = np.histogram(lnm, logMh)
Ms = 10**bins
dM = Ms[1:]-Ms[:-1]
dNdM = n/dM
Nc_WDM = np.add.accumulate(np.flip(n,0))
Nc_WDM = np.flip(Nc_WDM,0)
Ncs_WDM[j] = Nc_WDM
#pt.plot(Ms[1:],Nc_CDM,'blue',alpha=0.5)
#pt.plot(Ms[1:],Nc_WDM,'red',alpha=0.5)
# compute the mean over the Ngal "galaxies"
Nc_CDM_mean = Ncs_CDM.mean(axis=0)
Nc_WDM_mean = Ncs_WDM.mean(axis=0)
Nc_CDM_realisations[realisation] = Nc_CDM_mean
Nc_WDM_realisations[realisation] = Nc_WDM_mean
#pt.plot(Ms[1:],Nc_CDM_mean,'blue',alpha=0.5)
#pt.plot(Ms[1:],Nc_WDM_mean,'red',alpha=0.5)
Nc_CDM_realisations_mean = Nc_CDM_realisations.mean(axis=0)
Nc_CDM_realisations_std = Nc_CDM_realisations.std(axis=0)
Nc_CDM_meanp = Nc_CDM_realisations_mean + 3*Nc_CDM_realisations_std
Nc_CDM_meanm = Nc_CDM_realisations_mean - 3*Nc_CDM_realisations_std
Nc_WDM_realisations_mean = Nc_WDM_realisations.mean(axis=0)
Nc_WDM_realisations_std = Nc_WDM_realisations.std(axis=0)
Nc_WDM_meanp = Nc_WDM_realisations_mean + 3*Nc_WDM_realisations_std
Nc_WDM_meanm = Nc_WDM_realisations_mean - 3*Nc_WDM_realisations_std
ax = pt.gca()
ax.fill_between(Ms[1:],Nc_CDM_meanp,Nc_CDM_meanm,facecolor='blue',alpha=0.25,label=r"$\textrm{CDM}$")
ax.fill_between(Ms[1:],Nc_WDM_meanp,Nc_WDM_meanm,facecolor='red',alpha=0.25,label=r"$\textrm{WDM}$")
###########################
# finalize
###########################
xmin = opt.Mmin
xmax = opt.Mmax
ymin = None
ymax = None
xlabel = r"$\rm{Halo\,\,Mass\,\,(M_{\rm halo})}\,\,[M_{\odot}]$"
ylabel = r"$N_{\rm halo}(>M_{\rm halo})$"
pt.SetAxis(xmin,xmax,ymin,ymax,log="xy")
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.grid(False)
pt.title(r"$\rm{%d\,\,galaxies\,\,observed}$"%opt.Ngal,fontsize=pt.labelfont)
pt.legend()
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