Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F105418383
pychem_SSP_discret_evolution2
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
Mon, Mar 17, 01:57
Size
27 KB
Mime Type
text/x-python
Expires
Wed, Mar 19, 01:57 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
24963967
Attached To
rGEAR Gear
pychem_SSP_discret_evolution2
View Options
#!/usr/bin/env python
from optparse import OptionParser
from PyChem import chemistry
from numpy import *
import sys
import string
import Ptools as pt
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser = pt.add_postscript_options(parser)
parser.add_option("--x",
action="store",
dest="x",
type="string",
default = 'Fe',
help="x value to plot",
metavar=" STRING")
parser.add_option("--y",
action="store",
dest="y",
type="string",
default = 'Mg',
help="y value to plot",
metavar=" STRING")
parser.add_option("--dt",
action="store",
dest="dt",
type="float",
default = 0.1,
help="dt",
metavar=" FLOAT")
parser.add_option("--mstar",
action="store",
dest="mstar",
type="float",
default = 1e5,
help="initial mass of the SSP in solar mass",
metavar=" FLOAT")
parser.add_option("--tstar",
action="store",
dest="tstar",
type="float",
default = 0,
help="formation time of the SSP",
metavar=" FLOAT")
parser.add_option("-o",
action="store",
dest="obs",
type="string",
default = 'Y',
help="observable to plot",
metavar=" STRING")
parser.add_option("--timeunit",
action="store",
dest="timeunit",
type="string",
default = None,
help="unit of time",
metavar=" STRING")
parser.add_option("--NumberOfTables",
action="store",
dest="NumberOfTables",
type="int",
default = 1,
help="NumberOfTables",
metavar=" INT")
parser.add_option("--DefaultTable",
action="store",
dest="DefaultTable",
type="int",
default = 0,
help="DefaultTable",
metavar=" INT")
parser.add_option("--discsn",
action="store_true",
dest="discsn",
default = False,
help="discretize sn",
metavar=" FLOAT")
parser.add_option("--seed",
action="store",
dest="seed",
type="int",
default = None,
help="initial random seed",
metavar=" INT")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
########################################################################
# M A I N
########################################################################
SOLAR_MASS = 1.989e33
UnitLength_in_cm = 3.085e+21
UnitMass_in_g = 1.989e+43
UnitVelocity_in_cm_per_s = 20725573.785998672
UnitTime_in_s = 148849920000000.0
feedback_energy = 10**51 /UnitMass_in_g/(UnitVelocity_in_cm_per_s)**2
# parse options
files,opt = parse_options()
if opt.seed!=None:
random.seed(opt.seed)
# units
opt.ux = 1.0
if opt.timeunit == None:
opt.ux = 1.0
elif opt.timeunit == 'Myr':
opt.ux = UnitTime_in_s/(3600*24*365)/1e6
elif opt.timeunit == 'Gyr':
opt.ux = UnitTime_in_s/(3600*24*365)/1e9
# to plot
opt.obs = string.split(opt.obs,',')
# init chimie
chemistry.init_chimie(files[0],opt.NumberOfTables,opt.DefaultTable)
elts = chemistry.get_elts_labels()
nelts = chemistry.get_nelts()
print elts
# some parameters
dt = opt.dt
tmax = opt.xmax / opt.ux
t_star = opt.tstar # formation time of the SSP
M0 = opt.mstar *SOLAR_MASS/UnitMass_in_g # gas mass of the SSP (in code unit)
Mstacont = M0 # star mass
Mstaj = M0
Mstajcont = M0
Mgascont = 0.015*M0 # gas mass
Mgasj = 0.015*M0 # gas mass
Mgasjcont = 0.015*M0 # gas mass
# stellar initial element abondance
metals = zeros(chemistry.get_nelts(),float)
#metals[0] = 0.001771 *1
#metals[1] = 0.00091245 *1
#metals[2] = 0.0 *1
#metals[3] = 0.02 *1 # metal abondance (=metallicity z)
# gas initial element abondance
metalgcont = zeros(chemistry.get_nelts(),float)
metalgj = zeros(chemistry.get_nelts(),float)
metalgjcont = zeros(chemistry.get_nelts(),float)
# minimum live time of a star
minlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmax())
maxlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmin())
print minlivetime
print maxlivetime
# times
times = arange(dt,tmax,dt).astype(float)
nt = len(times)
# output arrays
TIMES = zeros(nt,float) # time
NSNIas = zeros(nt,float) # number of SNIa
NSNIIs = zeros(nt,float) # number of SNII
NSNIaconts = zeros(nt,float) # number of SNIa (continuous value)
NSNIIconts = zeros(nt,float) # number of SNII (continuous value)
NDYINs = zeros(nt,float) # number of dying stars
NDYINconts = zeros(nt,float) # number of dying stars
Es = zeros(nt,float) # injected energy
Mecont = zeros(nt,float) # ejected mass (continuous)
Mej = zeros(nt,float) # ejected mass (computed from NSN)
Mejcont = zeros(nt,float) # ejected mass (computed from NSN,continuous)
Mstaconts = zeros(nt,float) # star mass (continuous)
Mstajs = zeros(nt,float) # star mass
Mstajconts = zeros(nt,float) # star mass
Mgasconts = zeros(nt,float) # gas mass
Mgasjs = zeros(nt,float) # gas mass
Mgasjconts = zeros(nt,float) # gas mass
# elements to follow
# ejected
dataecont = {}
dataej = {}
dataejcont = {}
for elt in elts:
dataecont[elt] = zeros(nt,float)
dataej[elt] = zeros(nt,float)
dataejcont[elt] = zeros(nt,float)
# gas
datagcont = {}
datagj = {}
datagjcont = {}
for elt in elts:
datagcont[elt] = zeros(nt,float)
datagj[elt] = zeros(nt,float)
datagjcont[elt] = zeros(nt,float)
###############################################
# main loop
###############################################
for ti,time in enumerate(times):
tpdt = time + dt
if (tpdt - t_star >= minlivetime) and (tpdt - t_star < maxlivetime):
m1 = chemistry.star_mass_from_age(metals[-1],tpdt - t_star)
if (time - t_star >= minlivetime):
m2 = chemistry.star_mass_from_age(metals[-1],time - t_star)
else:
m2 = chemistry.get_Mmax()
if (m1>m2):
m1=m2
raise "m1>m2 !!!"
# number of SNIa
NSNIa = chemistry.SNIa_rate(m1,m2)*M0
# number of SNII
NSNII = chemistry.SNII_rate(m1,m2)*M0
# number of dying stars
NDYIN = chemistry.DYIN_rate(m1,m2)*M0
#########################
# discrete stuff
NSNIIcont = NSNII
NSNIacont = NSNIa
NDYINcont = NDYIN
if opt.discsn:
fNSNIa = NSNIa-floor(NSNIa) # fraction
NSNIa = floor(NSNIa) # integer
if random.random() < fNSNIa:
NSNIa = NSNIa+1
if opt.discsn:
fNSNII = NSNII-floor(NSNII) # fraction
NSNII = floor(NSNII) # integer
if random.random() < fNSNII:
NSNII = NSNII+1
if opt.discsn:
fNDYIN = NDYIN-floor(NDYIN) # fraction
NDYIN = floor(NDYIN) # integer
if random.random() < fNDYIN:
NDYIN = NDYIN+1
#
#########################
##################################################
# mass and yields ejection
##################################################
# standard way (continuous explosion)
EjectedMass = chemistry.Total_mass_ejection(m1,m2,M0,metals)
TotalEjectedEltMasscont = EjectedMass[2:]
TotalEjectedGasMasscont = EjectedMass[0]
# new way (discrete explosion)
EjectedMassj = chemistry.Total_single_mass_ejection(0.5*(m1+m2),metals,float(NSNII),float(NSNIa),float(NDYIN))
TotalEjectedEltMassj = EjectedMassj[2:]
TotalEjectedGasMassj = EjectedMassj[0]
# new way (continuous explosion)
EjectedMassjcont = chemistry.Total_single_mass_ejection(0.5*(m1+m2),metals,NSNIIcont,NSNIacont,NDYINcont)
TotalEjectedEltMassjcont = EjectedMassjcont[2:]
TotalEjectedGasMassjcont = EjectedMassjcont[0]
# other check
#SNII_EjectedMass = chemistry.SNII_Total_single_mass_ejection( 0.5*(m1+m2) ,metals) * NSNII
#SNIa_EjectedMass = chemistry.SNIa_Total_single_mass_ejection( 0.5*(m1+m2) ,metals) * NSNIa
#TotalEjectedEltMassjcont = SNII_EjectedMass[2:] + SNIa_EjectedMass[2:]
#TotalEjectedGasMassjcont = SNII_EjectedMass[0] + SNIa_EjectedMass[0]
'''
# all dying stars (= SNII + lighter masses)
DYIN_TotalEjectedGasMassj = chemistry.DYIN_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIa # !!! need to add SNIa
DYIN_TotalEjectedGasMassjcont = chemistry.DYIN_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIacont # !!! need to add SNIa
'''
##################################################
# energy
##################################################
# energy injected
TotalInjectedEnergy = feedback_energy* (NSNIa + NSNII)
# mass fraction of ejected elements
emetalcont = TotalEjectedEltMasscont /TotalEjectedGasMasscont
emetalj = TotalEjectedEltMassj /TotalEjectedGasMassj
emetaljcont= TotalEjectedEltMassjcont/TotalEjectedGasMassjcont
# metal enrichment
metalgcont = ((metalgcont *Mgascont) + TotalEjectedEltMasscont )
metalgj = ((metalgj *Mgasj) + TotalEjectedEltMassj )
metalgjcont = ((metalgjcont*Mgasjcont) + TotalEjectedEltMassjcont )
# Total Mass
Mstacont = Mstacont - TotalEjectedGasMasscont
Mstaj = Mstaj - TotalEjectedGasMassj
Mstajcont = Mstajcont - TotalEjectedGasMassjcont
# Gas Mass
Mgascont = Mgascont + TotalEjectedGasMasscont
Mgasj = Mgasj + TotalEjectedGasMassj
Mgasjcont = Mgasjcont + TotalEjectedGasMassjcont
# metal enrichment
metalgcont = metalgcont /Mgascont
metalgj = metalgj /Mgasj
metalgjcont = metalgjcont/Mgasjcont
# record output
TIMES[ti] = time
NSNIas[ti] = NSNIa
NSNIIs[ti] = NSNII
NSNIaconts[ti] = NSNIacont
NSNIIconts[ti] = NSNIIcont
NDYINs[ti] = NDYIN
NDYINconts[ti] = NDYINcont
Es[ti] = TotalInjectedEnergy
Mecont[ti] = TotalEjectedGasMasscont
Mej[ti] = TotalEjectedGasMassj
Mejcont[ti]= TotalEjectedGasMassjcont
Mstaconts[ti] = Mstacont # stellar mass
Mstajs[ti] = Mstaj # stellar mass
Mstajconts[ti] = Mstajcont # stellar mass
Mgasconts[ti] = Mgascont
Mgasjs[ti] = Mgasj
Mgasjconts[ti] = Mgasjcont
for j,elt in enumerate(elts):
dataecont[elt][ti] = emetalcont[j]
datagcont[elt][ti] = metalgcont[j]
dataej[elt][ti] = emetalj[j]
datagj[elt][ti] = metalgj[j]
dataejcont[elt][ti] = emetaljcont[j]
datagjcont[elt][ti] = metalgjcont[j]
# remove time=0 values
c = (TIMES!=0)
TIMES = compress(c,TIMES) * opt.ux
NSNIas = compress(c,NSNIas)
NSNIIs = compress(c,NSNIIs)
NSNIaconts = compress(c,NSNIaconts)
NSNIIconts = compress(c,NSNIIconts)
NDYINs = compress(c,NDYINs)
NDYINconts = compress(c,NDYINconts)
Es = compress(c,Es)
Mecont = compress(c,Mecont)
Mej = compress(c,Mej)
Mejcont = compress(c,Mejcont)
Mstaconts = compress(c,Mstaconts)
Mstajs = compress(c,Mstajs)
Mstajconts= compress(c,Mstajconts)
Mgasconts = compress(c,Mgasconts)
Mgasjs = compress(c,Mgasjs)
Mgasjconts= compress(c,Mgasjconts)
for elt in elts:
dataecont[elt] = compress(c,dataecont[elt])
datagcont[elt] = compress(c,datagcont[elt])
dataej[elt] = compress(c,dataej[elt])
datagj[elt] = compress(c,datagj[elt])
dataejcont[elt] = compress(c,dataejcont[elt])
datagjcont[elt] = compress(c,datagjcont[elt])
# get Solar Abundances
SolarAbun = chemistry.get_elts_SolarMassAbundances()
Xsconts = {}
Xsjs = {}
Xsjconts= {}
Xgconts = {}
Xgjs = {}
Xgjconts= {}
for elt in elts:
Xsconts[elt] = log10(dataecont[elt] / SolarAbun[elt] + 1.0e-10)
Xgconts[elt] = log10(datagcont[elt] / SolarAbun[elt] + 1.0e-10)
Xsjs[elt] = log10(dataej[elt] / SolarAbun[elt] + 1.0e-10)
Xgjs[elt] = log10(datagj[elt] / SolarAbun[elt] + 1.0e-10)
Xsjconts[elt] = log10(dataejcont[elt] / SolarAbun[elt] + 1.0e-10)
Xgjconts[elt] = log10(datagjcont[elt] / SolarAbun[elt] + 1.0e-10)
Xscont = Xsconts[opt.x]
Yscont = Xsconts[opt.y]
Xgcont = Xgconts[opt.x]
Ygcont = Xgconts[opt.y]
Xsj = Xsjs[opt.x]
Ysj = Xsjs[opt.y]
Xgj = Xgjs[opt.x]
Ygj = Xgjs[opt.y]
Xsjcont = Xsjconts[opt.x]
Ysjcont = Xsjconts[opt.y]
Xgjcont = Xgjconts[opt.x]
Ygjcont = Xgjconts[opt.y]
#########################################################
# plot
#########################################################
pt.InitPlot([],opt)
pt.pcolors
fig = pt.gcf()
#fig.subplots_adjust(left=0.02)
#fig.subplots_adjust(right=0.99)
#fig.subplots_adjust(bottom=0.1)
#fig.subplots_adjust(top=0.99)
#fig.subplots_adjust(wspace=0.0)
fig.subplots_adjust(hspace=0.25)
np = len(opt.obs)
def PlotYields(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{[%s/H],[%s/H],[%s/%s]}$"%(opt.x,opt.y,opt.y,opt.x)
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,Xscont,tpe='line',label='Xs',color='b')
datas.append(data)
data = pt.DataPoints(TIMES,Yscont,tpe='line',label='Ys',color='g')
datas.append(data)
data = pt.DataPoints(TIMES,Yscont-Xscont,tpe='line',label='Ys-Xs',color='r',linestyle='--')
datas.append(data)
data = pt.DataPoints(TIMES,Xsj,tpe='line',label='Xs',color='b',linestyle='--')
datas.append(data)
data = pt.DataPoints(TIMES,Ysj,tpe='line',label='Ys',color='g',linestyle='--')
datas.append(data)
data = pt.DataPoints(TIMES,Ysj-Xsj,tpe='line',label='Ys-Xs',color='r',linestyle='-')
datas.append(data)
data = pt.DataPoints(TIMES,Xsjcont,tpe='line',label='Xs',color='b')
#datas.append(data)
data = pt.DataPoints(TIMES,Ysjcont,tpe='line',label='Ys',color='g')
#datas.append(data)
data = pt.DataPoints(TIMES,Ysjcont-Xsjcont,tpe='line',label='Ys-Xs',color='r',linestyle='--')
#datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color,ls=d.linestyle)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
ymin = -2
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotGasElts(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{gas\,[%s/H],[%s/H],[%s/%s]}$"%(opt.x,opt.y,opt.y,opt.x)
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,Xgcont,tpe='line',label='Xs',color='b')
datas.append(data)
data = pt.DataPoints(TIMES,Ygcont,tpe='line',label='Ys',color='g')
datas.append(data)
data = pt.DataPoints(TIMES,Ygcont-Xgcont,tpe='line',label='Ys-Xs',color='r',linestyle='--')
datas.append(data)
data = pt.DataPoints(TIMES,Xgj,tpe='line',label='Xs',color='b',linestyle='--')
datas.append(data)
data = pt.DataPoints(TIMES,Ygj,tpe='line',label='Ys',color='g',linestyle='--')
datas.append(data)
data = pt.DataPoints(TIMES,Ygj-Xgj,tpe='line',label='Ys-Xs',color='r',linestyle='-')
datas.append(data)
data = pt.DataPoints(TIMES,Xgjcont,tpe='line',label='Xs',color='b')
#datas.append(data)
data = pt.DataPoints(TIMES,Ygjcont,tpe='line',label='Ys',color='g')
#datas.append(data)
data = pt.DataPoints(TIMES,Ygjcont-Xgjcont,tpe='line',label='Ys-Xs',color='r',linestyle='--')
#datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color,ls=d.linestyle)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
ymin = -2
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotEnergy(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Supernova\,Energy}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,Es,tpe='line',label='Es',color='k')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
if opt.log==None:
log='y'
else:
log = opt.log + 'y'
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotEjectedMass(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Ejected\,Mass}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,Mecont,tpe='line',label='Es',color='r')
datas.append(data)
data = pt.DataPoints(TIMES,Mej,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,Mejcont,tpe='line',label='Es',color='g')
#datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
if opt.log==None:
log=opt.log
else:
log = opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotCumEjectedMass(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Cumulative\,Ejected\,Mass}$"
pt.subplot(np,1,i)
cumcont = add.accumulate(Mecont)
cumj = add.accumulate(Mej)
cumjcont = add.accumulate(Mejcont)
datas = []
data = pt.DataPoints(TIMES,cumcont,tpe='line',label='Es',color='r')
datas.append(data)
data = pt.DataPoints(TIMES,cumj,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,cumjcont,tpe='line',label='Es',color='g')
#datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
if opt.log==None:
log=opt.log
else:
log = opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotStellarMass(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Stellar\,Mass}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,Mstaconts,tpe='line',label='Es',color='r')
datas.append(data)
data = pt.DataPoints(TIMES,Mstajs,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,Mstajconts,tpe='line',label='Es',color='g')
#datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
if opt.log==None:
log = opt.log
else:
log = opt.log
opt.ymin = 0
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotNDYIN(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Number\,dying\,stars}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,NDYINs,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,NDYINconts,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotCumNDYIN(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Cum\,Number\,of\,dying\,stars}$"
pt.subplot(np,1,i)
cum = add.accumulate(NDYINs)
cumc= add.accumulate(NDYINconts)
datas = []
data = pt.DataPoints(TIMES,cum,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,cumc,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotNSNII(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Number\,of\,SNII}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,NSNIIs,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,NSNIIconts,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotCumNSNII(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Cum\,Number\,of\,SNII}$"
pt.subplot(np,1,i)
cum = add.accumulate(NSNIIs)
cumc= add.accumulate(NSNIIconts)
datas = []
data = pt.DataPoints(TIMES,cum,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,cumc,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotNSNIa(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Number\,of\,SNIa}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,NSNIas,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,NSNIaconts,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotCumNSNIa(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Cum\,Number\,of\,SNIa}$"
pt.subplot(np,1,i)
cum = add.accumulate(NSNIas)
cumc= add.accumulate(NSNIaconts)
datas = []
data = pt.DataPoints(TIMES,cum,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,cumc,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotNSNIIDYIN(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Number\,of\,SNII\,and\,dying}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,NSNIIs+NDYINs,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,NSNIIconts+NDYINconts,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotFevsTime(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Fe}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,Ygj,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,Ygjcont,tpe='line',label='Es',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log=opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotDiff(i):
xlabel = r"$\rm{[%s/H]$"%(opt.x)
ylabel = r"$\rm{[%s/%s]}$"%(opt.y,opt.x)
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(Xgcont,Ygcont-Xgcont ,tpe='line',label='Xs',color='r')
datas.append(data)
data = pt.DataPoints(Xgj,Ygj-Xgj ,tpe='line',label='Xs',color='k')
datas.append(data)
#data = pt.DataPoints(Xgjcont,Ygjcont-Xgjcont,tpe='line',label='Xs',color='g')
#datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log = opt.log
#xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
xmin = -4
xmax = 0
ymin = -1
ymax = 2
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotStarMass(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Mass}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,Mstas,tpe='line',label='Mstas',color='k')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
log = opt.log
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pdict = {'Y':PlotYields,'Yg':PlotGasElts,'ESN':PlotEnergy,'NSNII':PlotNSNII,'CumNSNII':PlotCumNSNII,'NSNIa':PlotNSNIa,'CumNSNIa':PlotCumNSNIa,'NDYIN':PlotNDYIN,'CumNDYIN':PlotCumNDYIN,
'D':PlotDiff,'MSTAR':PlotStarMass,'EjMass':PlotEjectedMass,'CumEjMass':PlotCumEjectedMass,'StellarMass':PlotStellarMass,
'NSNIIDYIN':PlotNSNIIDYIN,'Fe':PlotFevsTime}
ip = 1
for obs in opt.obs:
pdict[obs](ip)
ip = ip + 1
pt.EndPlot(files,opt)
Event Timeline
Log In to Comment