Page MenuHomec4science

pychem_SSP_evolution
No OneTemporary

File Metadata

Created
Sat, Apr 19, 12:51

pychem_SSP_evolution

#!/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 = 0,
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()
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
Mgas = 10*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
metalg = 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)
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
Mgass = zeros(nt,float) # gas mass
# elements to follow
# ejected
datae = {}
for elt in elts:
datae[elt] = zeros(nt,float)
# gas
datag = {}
for elt in elts:
datag[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
#########################
# SN trick
NSNIIcont = NSNII
NSNIacont = NSNIa
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
#
#########################
EjectedMass = chemistry.Total_mass_ejection(m1,m2,M0,metals)
TotalEjectedEltMass = EjectedMass[2:]
TotalEjectedGasMasscont = EjectedMass[0] # continuous value !!! inconsitent with the next lines, as
SNII_TotalEjectedGasMassj = chemistry.SNII_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNII # !!! need to add SNIa
SNII_TotalEjectedGasMassjcont = chemistry.SNII_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIIcont # !!! need to add SNIa
SNIa_TotalEjectedGasMassj = chemistry.SNIa_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIa # !!! need to add SNIa
SNIa_TotalEjectedGasMassjcont = chemistry.SNIa_Total_single_mass_ejection( 0.5*(m1+m2) ,metals)[0] * NSNIacont # !!! need to add SNIa
TotalEjectedGasMassj = SNII_TotalEjectedGasMassj + SNIa_TotalEjectedGasMassj
TotalEjectedGasMassjcont = SNII_TotalEjectedGasMassjcont+ SNIa_TotalEjectedGasMassjcont
# energy injected
TotalInjectedEnergy = feedback_energy* (NSNIa + NSNII)
# mass fraction of ejected elements
emetal = TotalEjectedEltMass/TotalEjectedGasMasscont
# metal enrichment
metalg = ((metalg*Mgas) + TotalEjectedEltMass )
# Total Mass
Mstacont = Mstacont - TotalEjectedGasMasscont
Mstaj = Mstaj - TotalEjectedGasMassj
Mstajcont = Mstajcont - TotalEjectedGasMassjcont
# Gas Mass
Mgas = Mgas + TotalEjectedGasMasscont
# metal enrichment
metalg = metalg/Mgas
# record output
TIMES[ti] = time
NSNIas[ti] = NSNIa
NSNIIs[ti] = NSNII
NSNIaconts[ti] = NSNIacont
NSNIIconts[ti] = NSNIIcont
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
Mgass[ti] = Mgas
for j,elt in enumerate(elts):
datae[elt][ti] = emetal[j]
datag[elt][ti] = metalg[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)
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)
Mgass = compress(c,Mgass)
for elt in elts:
datae[elt] = compress(c,datae[elt])
datag[elt] = compress(c,datag[elt])
# get Solar Abundances
SolarAbun = chemistry.get_elts_SolarAbundances()
Xss = {}
Xgs = {}
for elt in elts:
Xss[elt] = log10(datae[elt] / SolarAbun[elt] + 1.0e-10)
Xgs[elt] = log10(datag[elt] / SolarAbun[elt] + 1.0e-10)
Xs = Xss[opt.x]
Ys = Xss[opt.y]
Xg = Xgs[opt.x]
Yg = Xgs[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,Xs,tpe='line',label='Xs',color='b')
datas.append(data)
data = pt.DataPoints(TIMES,Ys,tpe='line',label='Ys',color='g')
datas.append(data)
data = pt.DataPoints(TIMES,Ys-Xs,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:
opt.log='y'
else:
opt.log = opt.log + 'y'
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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:
opt.log='y'
else:
opt.log = opt.log + 'y'
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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:
opt.log='y'
else:
opt.log = opt.log + 'y'
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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\,SN}$"
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)
opt.log=None
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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\,SN}$"
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)
opt.log=None
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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\,SN}$"
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)
opt.log=None
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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\,SN}$"
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)
opt.log=None
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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(Xg,Yg-Xg,tpe='line',label='Xs',color='b')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
opt.log = None
#xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
xmin = -4
xmax = 0
ymin = -1
ymax = 2
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.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)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pdict = {'Y':PlotYields,'ESN':PlotEnergy,'NSNII':PlotNSNII,'CumNSNII':PlotCumNSNII,'NSNIa':PlotNSNIa,'CumNSNIa':PlotCumNSNIa,'D':PlotDiff,'MSTAR':PlotStarMass,'EjMass':PlotEjectedMass,'StellarMass':PlotStellarMass}
ip = 1
for obs in opt.obs:
pdict[obs](ip)
ip = ip + 1
pt.EndPlot(files,opt)
#######################################
#pt.subplot(3,1,1)
#pt.plot(TIMES,Xs)
#pt.plot(TIMES,Ys)
#pt.plot(TIMES,Ys-Xs,'--')
##pt.semilogx()
#pt.axis([dt,tmax,-3,2])
#pt.xlabel('Time')
#pt.ylabel('%s,%s,%s-%s'%(opt.x,opt.y,opt.y,opt.x))
#######################################
#pt.subplot(3,1,2)
##pt.plot(TIMES,Yg-Xg)
##pt.axis([dt,tmax,-4,2])
#pt.plot(TIMES,datae['Mg']/(datae['Fe']+1e-20))
##pt.plot(TIMES,datae['Mg'])
#pt.semilogy()
##pt.axis([dt,tmax,1e-])
#pt.xlabel('Time')
#pt.ylabel('%s-%s'%(opt.y,opt.x))
#######################################
#pt.subplot(3,1,3)
#pt.plot(Xg,Yg-Xg)
#pt.axis([-4,0,-4,2])
#pt.xlabel('%s'%opt.x)
#pt.ylabel('%s-%s'%(opt.y,opt.x))
#'''
#######################################
#pt.subplot(6,1,1)
#pt.plot(TIMES,NSNIIs)
#pt.plot(TIMES,NSNIas)
#pt.semilogy()
#pt.semilogx()
#pt.axis([dt,tmax,1e-5,1e3])
#pt.xlabel('Time')
#pt.ylabel('Number of SN')
#######################################
#pt.subplot(6,1,2)
#pt.plot(TIMES,Es)
#pt.semilogy()
#pt.semilogx()
#pt.axis([dt,tmax,1e-12,1e-5])
#pt.xlabel('Time')
#pt.ylabel('Energy')
#######################################
#pt.subplot(6,1,3)
#pt.plot(TIMES,Xs)
#pt.plot(TIMES,Ys)
#pt.plot(TIMES,Ys-Xs,'--')
##pt.semilogx()
#pt.axis([dt,tmax,-3,2])
#pt.xlabel('Time')
#pt.ylabel('%s,%s,%s-%s'%(opt.x,opt.y,opt.y,opt.x))
#######################################
#pt.subplot(6,1,4)
#pt.plot(TIMES,Mstas)
##pt.plot(TIMES,Mgass)
##pt.plot(TIMES,Mgass+Mstas)
#pt.semilogx()
##pt.axis([dt,tmax,0,2.2*M0])
#pt.xlabel('Time')
#pt.ylabel('Total mass')
#######################################
#pt.subplot(6,1,5)
#pt.plot(TIMES,Xg)
#pt.plot(TIMES,Yg)
#pt.plot(TIMES,Yg-Xg)
#pt.semilogx()
#pt.axis([dt,tmax,-10,2])
#pt.xlabel('Time')
#pt.ylabel('%s,%s'%(opt.x,opt.y))
#######################################
#pt.subplot(6,1,6)
#pt.plot(Xg,Yg-Xg)
#pt.axis([-4,0,-4,2])
#pt.xlabel('%s'%opt.x)
#pt.ylabel('%s-%s'%(opt.y,opt.x))
#'''
#pt.show()

Event Timeline