Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_SSP_discret_evolution3
e0b388c18b0amaster

Last Change
Raw File

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
print "list of elements"
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())
print
print "max and min stellar mass"
print chemistry.get_Mmax(),chemistry.get_Mmax()*1e10
print chemistry.get_Mmin(),chemistry.get_Mmin()*1e10
print
print "min and max livetimes"
print minlivetime, minlivetime*UnitTime_in_s/(31536000)/1e6
print maxlivetime, maxlivetime*UnitTime_in_s/(31536000)/1e6
print
print "max and min stellar mass by inverting age"
print chemistry.star_mass_from_age(0,minlivetime)
print chemistry.star_mass_from_age(0,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
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 NSNII > 0:
# print NSNII,0.5*(m1+m2) *1e10
#
#########################$
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)
'''
datas = []
if opt.plot_discrete:
data = pt.DataPoints(TIMES,Mstajs,tpe='line',label='Es',color='k')
datas.append(data)
if opt.plot_continuous:
data = pt.DataPoints(TIMES,Mstaconts,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 = 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)
c4science · Help