Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F114297714
pychem_SSP_discret_evolution3
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
Sun, May 25, 00:04
Size
36 KB
Mime Type
text/x-python
Expires
Tue, May 27, 00:04 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
26379547
Attached To
rGEAR Gear
pychem_SSP_discret_evolution3
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("--X",
action="store",
dest="X",
type="string",
default = None,
help="x value to plot",
metavar=" STRING")
parser.add_option("--Y",
action="store",
dest="Y",
type="string",
default = None,
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("--tmin",
action="store",
dest="tmin",
type="float",
default = 1,
help="tmin",
metavar=" FLOAT")
parser.add_option("--tmax",
action="store",
dest="tmax",
type="float",
default = 14000,
help="tmax",
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("--mgasfact",
action="store",
dest="mgasfact",
type="float",
default = 1.5,
help="ratio between mgas and mstars",
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("--seed",
action="store",
dest="seed",
type="int",
default = None,
help="initial random seed",
metavar=" INT")
parser.add_option("-Z","--Z",
action="store",
dest="Z",
type="float",
default = 0.02,
help="metallicity",
metavar=" FLOAT")
parser.add_option("--plot_continuous_only",
action="store_true",
dest="plot_continuous_only",
default = False,
help="plot continuous values only",
metavar=" FLOAT")
parser.add_option("--plot_discrete_only",
action="store_true",
dest="plot_discrete_only",
default = False,
help="plot discrete values only",
metavar=" FLOAT")
parser.add_option("--optimal_sampling",
action="store_true",
dest="optimal_sampling",
default = False,
help="use optimal sampling for the discrete case",
metavar=" FLOAT")
(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,',')
if opt.X!=None:
opt.X = string.split(opt.X,',')
if opt.Y!=None:
opt.Y = string.split(opt.Y,',')
# 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.tmax / 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 = opt.mgasfact*M0 # gas mass
Mgasj = opt.mgasfact*M0 # gas mass
Mgasjcont = opt.mgasfact*M0 # gas mass
# stellar initial element abondance
metals = zeros(chemistry.get_nelts(),float)
metals[-1]= opt.Z
#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())
# 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
Escont = zeros(nt,float) # injected energy
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 = {}
Meltconts = {} # ejected elt mass
Meltjs = {} # ejected elt mass
Meltjconts = {} # ejected elt mass
for elt in elts:
dataecont[elt] = zeros(nt,float)
dataej[elt] = zeros(nt,float)
dataejcont[elt] = zeros(nt,float)
Meltconts[elt] = zeros(nt,float) # ejected elt mass
Meltjs[elt] = zeros(nt,float) # ejected elt mass
Meltjconts[elt] = zeros(nt,float) # ejected elt mass
# gas
datagcont = {}
datagj = {}
datagjcont = {}
for elt in elts:
datagcont[elt] = zeros(nt,float)
datagj[elt] = zeros(nt,float)
datagjcont[elt] = zeros(nt,float)
# optimal
Mecl = opt.mstar # in Msol
m_max_star = chemistry.optimal_sampling_init_norm(Mecl)
current_mass = m_max_star # in Msol
M_diced = 0
##############################################################################################
#
# 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 !!!"
############################
# optimal sampling stuffs
if opt.optimal_sampling:
NSNII_opt = 0
NSNIa_opt = 0
NDYIN_opt = 0
RSNIa=0
while ( current_mass*SOLAR_MASS/UnitMass_in_g > m1 ):
if chemistry.optimal_sampling_stop_loop(current_mass):
print "optimal_sampling_stop_loop !!!"
sys.exit()
#print "current mass = %g m1=%g m2=%g"%(current_mass,m1*1e10,m2*1e10)
if ( (chemistry.get_SNII_Mmin() <= current_mass*SOLAR_MASS/UnitMass_in_g) and (current_mass*SOLAR_MASS/UnitMass_in_g <= chemistry.get_SNII_Mmax()) ):
#print m1*1e10,m2*1e10,"one SNII"
NSNII_opt = NSNII_opt+1
if ( (chemistry.get_DYIN_Mmin() <= current_mass*SOLAR_MASS/UnitMass_in_g) and (current_mass*SOLAR_MASS/UnitMass_in_g <= chemistry.get_DYIN_Mmax()) ):
#print m1*1e10,m2*1e10,"one DYIN"
NDYIN_opt = NDYIN_opt+1
RSNIa = RSNIa + chemistry.SNIa_single_rate(current_mass*SOLAR_MASS/UnitMass_in_g)
current_mass = chemistry.optimal_sampling_get_next_mass(current_mass)
# because the function SNIa_single_rate returns a continuous value,
# we need to discretize it
if RSNIa>0:
fNSNIa = RSNIa-floor(RSNIa)
NSNIa = floor(RSNIa)
if random.random() < fNSNIa:
NSNIa = NSNIa+1
NSNIa_opt = NSNIa
# 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
# compute discrete values (SNIa)
fNSNIa = NSNIa-floor(NSNIa) # fraction
NSNIa = floor(NSNIa) # integer
if random.random() < fNSNIa:
NSNIa = NSNIa+1
# compute discrete values (SNII)
fNSNII = NSNII-floor(NSNII) # fraction
NSNII = floor(NSNII) # integer
if random.random() < fNSNII:
NSNII = NSNII+1
# compute discrete values (DYIN)
fNDYIN = NDYIN-floor(NDYIN) # fraction
NDYIN = floor(NDYIN) # integer
if random.random() < fNDYIN:
NDYIN = NDYIN+1
#
#########################$
if opt.optimal_sampling:
NSNII = NSNII_opt
NSNIa = NSNIa_opt
NDYIN = NDYIN_opt
##################################################
# mass and yields ejection
##################################################
# standard way (continuous explosion) (_cont)
EjectedMass = chemistry.Total_mass_ejection_Z(m1,m2,M0,metals)
TotalEjectedEltMasscont = EjectedMass[2:]
TotalEjectedGasMasscont = EjectedMass[0]
# new way (discrete explosion) (_j)
EjectedMassj = chemistry.Total_single_mass_ejection_Z(0.5*(m1+m2),metals,float(NSNII),float(NSNIa),float(NDYIN))
TotalEjectedEltMassj = EjectedMassj[2:]
TotalEjectedGasMassj = EjectedMassj[0]
# new way (continuous explosion) (_jcont)
EjectedMassjcont = chemistry.Total_single_mass_ejection_Z(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
TotalInjectedEnergycont = feedback_energy* (NSNIacont + NSNIIcont)
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
Escont[ti] = TotalInjectedEnergycont
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]
Meltconts[elt][ti] = TotalEjectedEltMasscont[j]
Meltjs[elt][ti] = TotalEjectedEltMassj[j]
Meltjconts[elt][ti] = TotalEjectedEltMassjcont[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)
Escont = compress(c,Escont)
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])
Meltconts[elt] = compress(c,Meltconts[elt])
Meltjs[elt] = compress(c,Meltjs[elt])
Meltjconts[elt] = compress(c,Meltjconts[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
#
##################################################################################################################
# some plotting options
opt.plot_continuous=True
opt.plot_discrete=True
opt.plot_semicontnuous=False
if opt.plot_continuous_only:
opt.plot_continuous=True
opt.plot_discrete=False
opt.plot_semicontnuous=False
if opt.plot_discrete_only:
opt.plot_continuous=False
opt.plot_discrete=True
opt.plot_semicontnuous=False
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 = []
if opt.plot_continuous:
data = pt.DataPoints(TIMES,Xscont,tpe='line',label='Xs',color='b',linestyle='-')
datas.append(data)
data = pt.DataPoints(TIMES,Yscont,tpe='line',label='Ys',color='g',linestyle='-')
datas.append(data)
data = pt.DataPoints(TIMES,Yscont-Xscont,tpe='line',label='Ys-Xs',color='r',linestyle='-',linewidth=2)
datas.append(data)
if opt.plot_discrete:
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='--',linewidth=2)
datas.append(data)
if opt.plot_semicontnuous:
data = pt.DataPoints(TIMES,Xsjcont,tpe='line',label='Xs',color='b',linestyle=':')
datas.append(data)
data = pt.DataPoints(TIMES,Ysjcont,tpe='line',label='Ys',color='g',linestyle=':')
datas.append(data)
data = pt.DataPoints(TIMES,Ysjcont-Xsjcont,tpe='line',label='Ys-Xs',color='r',linestyle=':',linewidth=2)
datas.append(data)
for d in datas:
pt.plot(d.x,d.y,color=d.color,ls=d.linestyle,lw=d.linewidth)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
ymin = max(-2,ymin)
# 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,Escont,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,Es,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)
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 PlotCumEnergy(i):
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Supernova\,Energy}$"
pt.subplot(np,1,i)
CumEscont = add.accumulate(Escont)
CumEs = add.accumulate(Es)
datas = []
data = pt.DataPoints(TIMES,CumEscont,tpe='line',label='Es',color='k')
datas.append(data)
data = pt.DataPoints(TIMES,CumEs,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)
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 = []
if opt.plot_discrete:
data = pt.DataPoints(TIMES,NSNIIs,tpe='line',label='Es',color='k')
datas.append(data)
if opt.plot_continuous:
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 = []
if opt.plot_discrete:
data = pt.DataPoints(TIMES,cum,tpe='line',label='Es',color='k')
datas.append(data)
if opt.plot_continuous:
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 = []
if opt.plot_discrete:
data = pt.DataPoints(TIMES,NSNIas,tpe='line',label='Es',color='k')
datas.append(data)
if opt.plot_continuous:
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 PlotMassFevsTime(i):
'''
!!! this is not very well written. Need to be improved
'''
xlabel = r'$\rm{Time} [\rm{%s}]$'%(opt.timeunit)
ylabel = r"$\rm{Fe}$"
pt.subplot(np,1,i)
Xgj = add.accumulate(Meltconts['Fe'])
Ygj = add.accumulate(Meltconts['Mg'])
#Xgj = log10(Meltconts['Mg']+1e-20)
#Ygj = log10(Meltconts['Fe']+1e-20)
pt.plot(TIMES,Xgj)
pt.plot(TIMES,Ygj)
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=None
#xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
xmin = 0
xmax = 14000
ymin = -20
ymax = -8
# 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 = []
if opt.plot_continuous:
data = pt.DataPoints(Xgcont,Ygcont-Xgcont ,tpe='line',label='Xs',color='r')
datas.append(data)
if opt.plot_discrete:
data = pt.DataPoints(Xgj,Ygj-Xgj ,tpe='line',label='Xs',color='k')
datas.append(data)
if opt.plot_semicontnuous:
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 = -2
ymax = 2
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
def PlotDXvsDY(i):
# X = X/H
X1 = opt.X[0]
X2 = opt.X[1]
print X1,X2
Y1 = opt.Y[0]
Y2 = opt.Y[1]
print Y1,Y2
if X2=='H':
xcont = Xgconts[X1]
xj = Xgjs[X1]
xjcont = Xgjconts[X1]
xlabel = r"$\rm{[%s/%s]}$"%(X1,"H")
else:
xcont = Xgconts[X1] - Xgconts[X2]
xj = Xgjs[X1] - Xgjs[X2]
xjcont = Xgjconts[X1] - Xgjconts[X2]
xlabel = r"$\rm{[%s/%s]}$"%(X1,X2)
#Xgcont,Ygcont-Xgcont
#Xgconts[opt.x]
if Y2=='H':
ycont = Xgconts[Y1]
yj = Xgjs[Y1]
yjcont = Xgjconts[Y1]
ylabel = r"$\rm{[%s/%s]}$"%(Y1,"H")
else:
ycont = Xgconts[Y1] - Xgconts[Y2]
yj = Xgjs[Y1] - Xgjs[Y2]
yjcont = Xgjconts[Y1] - Xgjconts[Y2]
ylabel = r"$\rm{[%s/%s]}$"%(Y1,Y2)
pt.subplot(np,1,i)
datas = []
if opt.plot_continuous:
data = pt.DataPoints(xcont,ycont, tpe='line',label='Xs',color='r')
datas.append(data)
if opt.plot_discrete:
data = pt.DataPoints(xj,yj, tpe='line',label='Xs',color='k')
datas.append(data)
if opt.plot_semicontnuous:
data = pt.DataPoints(xjcont,yjcont,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 = None
xmin = opt.xmin
xmax = opt.xmax
ymin = opt.ymin
ymax = opt.ymax
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(xmin,xmax,ymin,ymax,datas,log)
print xlabel
print ylabel
# 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,'CumESN':PlotCumEnergy,'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,'DXvsDY':PlotDXvsDY,'MFe':PlotMassFevsTime}
ip = 1
for obs in opt.obs:
pdict[obs](ip)
ip = ip + 1
pt.EndPlot(files,opt)
Event Timeline
Log In to Comment