Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_plot_ejectats_vs_time
d89b118970f3master

Last Change
Raw File

pychem_plot_ejectats_vs_time
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
from pNbody import units
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")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
########################################################################
# M A I N
########################################################################
SOLAR_MASS = 1.989e33
params={}
params['UnitLength_in_cm']=3.085e+21
params['UnitMass_in_g']=1.989e+43
params['UnitVelocity_in_cm_per_s']=20725573.785998672
system_of_units = units.Set_SystemUnits_From_Params(params)
out_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_Myr,units.Unit_K])
toMyrs= system_of_units.convertionFactorTo(out_units.UnitTime)
# parse options
files,opt = parse_options()
# init chimie
chemistry.init_chimie(files[0],opt.NumberOfTables,opt.DefaultTable)
# some inits
t_star=0
M0 = 1.
metals = zeros(chemistry.get_nelts(),float)
# minimum live time of a star
maxSNIIlivetime = chemistry.star_lifetime(metals[0],chemistry.get_SNII_Mmin())
# minimum live time of a star
minlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmax())
maxlivetime = chemistry.star_lifetime(metals[-1],chemistry.get_Mmin())
#tmax = maxSNIIlivetime
tmax = maxSNIIlivetime
dt = tmax/1000.
times = arange(dt,tmax,dt).astype(float)
n = len(times)
mFes = zeros(n)
mMgs = zeros(n)
mZs = zeros(n)
mEjs = zeros(n)
Ms = zeros(n)
###############################################
# 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 !!!"
m2 = chemistry.get_Mmax()
# standard way (continuous explosion)
EjectedMass = chemistry.Total_mass_ejection(m1,m2,M0,metals)
TotalEjectedEltMasscont = EjectedMass[2:]
TotalEjectedGasMasscont = EjectedMass[0]
#print TotalEjectedEltMasscont[0],TotalEjectedEltMasscont[-1],m1,m2
mEjs[ti]=TotalEjectedGasMasscont
mFes[ti]=TotalEjectedEltMasscont[0]
mMgs[ti]=TotalEjectedEltMasscont[1]
mZs[ti] =TotalEjectedEltMasscont[-1]
Ms[ti]=m2
SolarAbun = chemistry.get_elts_SolarAbundances()
MgFe = log10((mMgs/mFes)/(SolarAbun['Mg']/SolarAbun['Fe']))
pt.plot(times*toMyrs,mFes,label='Fe')
pt.plot(times*toMyrs,mMgs,label='Mg')
pt.plot(times*toMyrs,mZs,label='Metals')
pt.plot(times*toMyrs,mEjs,label='Ejected Mass')
#pt.plot(times*toMyrs,MgFe)
pt.xlabel('Time [Myrs]')
pt.ylabel('Mass fraction')
pt.legend(loc="upper left")
pt.show()
print system_of_units.convertionFactorTo(out_units.UnitTime)
print "Fe",mFes[-1]
print "Mg",mMgs[-1]
print "Z",mZs[-1]
c4science · Help