Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_SSP_discret_evolution2
f7e36798d4e7master

Last Change
Raw File

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