satComputeNCumvsLvWithErrors_forOliver.py
0d50a2737ed2master

satComputeNCumvsLvWithErrors_forOliver.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("--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')


			parser.add_argument("--filter",
			action="store",
			dest="filter",
			metavar='STR',
			type=str,
			default=None,
			help='filter name')

			parser.add_argument("--distance",
			action="store",
			dest="distance",
			type=float,
			default=35,
			help="distance in Mpc",
			metavar=" FLOAT")

			####################################################################
			# 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


			# 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)



			# number of haloes for the CDM
			Nh_CDM = CDM.NumberOfHaloes()


			Nbins = 100
			Nsats_CDM = np.zeros((opt.Ngal,Nbins-1))

			Ns_CDM_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 = CDM.fctSample(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





			Ns_CDM_mean = Nsats_CDM.mean(axis=0)

			Ns_CDM_realisations[realisation] = Ns_CDM_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



			if opt.filter is not None:
			# convert Lv assumed to be Msol in Mag

			Mag = np.zeros(len(Lv))
			for i in range(len(Lv)):
			Mag[i] = satlib.Msol2Mag(10**Lv[i],filter=opt.filter)

			# to apparent magnitude
			opt.distance = opt.distance * 1e6 # Mpc to pc
			Mag = Mag + 5*np.log10(opt.distance) - 5


			ax = pt.gca()
			ax.fill_between(Mag,Ns_CDM_meanp,Ns_CDM_meanm,facecolor='blue',alpha=0.5,label=r"$\textrm{CDM}$")

			xmin = 30
			xmax = 15
			ymin = 1
			ymax = 500

			pt.SetAxis(xmin,xmax,ymin,ymax,log="y")
			xlabel = r"$\rm{%s\,\,mag}$"%(opt.filter)

			# 1e5 Msol
			pt.plot([23.2,23.2],[0.1,1000],c='k',ls=':')


			else:
			ax = pt.gca()
			ax.fill_between(10**Lv,Ns_CDM_meanp,Ns_CDM_meanm,facecolor='blue',alpha=0.5,label=r"$\textrm{CDM}$")

			xmin = 5e2
			xmax = 5e9
			ymin = 1
			ymax = 500

			pt.SetAxis(xmin,xmax,ymin,ymax,log="xy")
			xlabel = r"$\rm{V-band\,\,Luminosity}$"

			pt.plot([1e5,1e5],[0.1,1000],c='k',ls=':')


			'''

			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)

			'''

			# Oliver data
			Lv = np.array([5.91561634e+06, 2.14783047e+07, 8.55066713e+07, 1.12719746e+06,
			3.73250158e+07, 1.95884467e+08, 2.83139200e+06, 3.73250158e+07,
			1.78648757e+06, 1.02801630e+08, 2.83139200e+05, 4.09260660e+08,
			5.91561634e+06])

			logLv = np.log10(Lv)

			bins = np.array([5,5.5,6,6.5,7,7.5,8])
			bins = np.linspace(4,9,20)

			h,bins = np.histogram(logLv, bins=bins)

			h = np.flip(h)
			h = np.add.accumulate(h)
			h = np.flip(h)

			#pt.scatter(10**bins[1:],h,marker='s',color="k",s=50,label=r'$\rm{M83}$')
			pt.plot(10**bins[1:],h,color="k",label=r'$\rm{M83}$')



			###########################
			# finalize
			###########################





			ylabel = r"$\rm{Cumulative\,\,number\,\,of\,\,galaxies}\,\,(D<300\,\rm{kpc})$"


			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=(82,62))
			# 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)

arrakihs/Workspace/CumNumOfSat/satComputeNCumvsLvWithErrors_forOliver.py0d50a2737ed2master

satComputeNCumvsLvWithErrors_forOliver.pyView Options

arrakihs/Workspace/CumNumOfSat/satComputeNCumvsLvWithErrors_forOliver.py
0d50a2737ed2master

satComputeNCumvsLvWithErrors_forOliver.py
View Options