Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F110191357
pychem_get_IMF_info_and_ejectats
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
Fri, Apr 25, 02:47
Size
6 KB
Mime Type
text/x-python
Expires
Sun, Apr 27, 02:47 (2 d)
Engine
blob
Format
Raw Data
Handle
25792673
Attached To
rGEAR Gear
pychem_get_IMF_info_and_ejectats
View Options
#!/usr/bin/env python
import sys
from PyChem import chemistry
from numpy import *
from pNbody import units
from optparse import OptionParser
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("--t1",
action="store",
dest="t1",
type="float",
default = None,
help="t1 in Myrs",
metavar=" FLOAT")
parser.add_option("--t2",
action="store",
dest="t2",
type="float",
default = None,
help="t2 in Myrs",
metavar=" FLOAT")
parser.add_option("--M0",
action="store",
dest="M0",
type="float",
default = 1,
help="mass of SSP in Msol",
metavar=" FLOAT")
parser.add_option("--Mgas",
action="store",
dest="Mgas",
type="float",
default = 1,
help="mass of gas reservoir in Msol",
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)
toMsol= system_of_units.convertionFactorTo(out_units.UnitMass)
# parse options
files,opt = parse_options()
# init chimie
chemistry.init_chimie(files[0],opt.NumberOfTables,opt.DefaultTable)
# assume SSP with zero metallicity
metals = zeros(chemistry.get_nelts(),float)
tmin = chemistry.star_lifetime(metals[0],chemistry.get_Mmax())
tmax = chemistry.star_lifetime(metals[0],chemistry.get_Mmin())
if opt.t1!=None:
opt.t1=opt.t1/toMyrs
opt.t1=max(tmin,opt.t1)
m2 = chemistry.star_mass_from_age(metals[-1],opt.t1)
else:
m2 = chemistry.get_Mmax()
if opt.t2!=None:
opt.t2=opt.t2/toMyrs
opt.t2=min(tmax,opt.t2)
m1 = chemistry.star_mass_from_age(metals[-1],opt.t2)
else:
m1 = chemistry.get_Mmin()
livetime1 = chemistry.star_lifetime(metals[0],m1)
livetime2 = chemistry.star_lifetime(metals[0],m2)
labels = chemistry.get_elts_labels()
SolarAbun = chemistry.get_elts_SolarAbundances()
####################################################
# global info
####################################################
print
print "min star mass = %5.2f Msol livetime = %5.2f Myrs"%(m1*toMsol,livetime1*toMyrs)
print "max star mass = %5.2f Msol livetime = %5.2f Myrs"%(m2*toMsol,livetime2*toMyrs)
print
print "initial SSP mass = %g Msol"%(opt.M0)
print "gas reservoir mass = %g Msol"%(opt.Mgas)
print
##############################
# Ejection from all stars
##############################
EjectedMass = chemistry.Total_mass_ejection(m1,m2,1,metals)
TotalEjectedMass = EjectedMass[0]*opt.M0
TotalEjectedEltMass = EjectedMass[2:]*opt.M0
print "Total ejected mass = %9.3g Msol"%(TotalEjectedMass)
for i,label in enumerate(labels):
print " Ejected mass in %07s = %9.3g Msol"%(label,TotalEjectedEltMass[i])
print
for i,label in enumerate(labels):
XH = log10( (TotalEjectedEltMass[i]/(opt.Mgas+TotalEjectedMass)) / SolarAbun[label] )
print " [%07s/H] = %9.3g"%(label,XH)
print
##############################
# Ejection from SNII
##############################
EjectedMass = chemistry.SNII_mass_ejection(m1,m2)
TotalEjectedMass = EjectedMass[0]*opt.M0
TotalEjectedEltMass = EjectedMass[2:]*opt.M0
print "Total ejected mass (SNII) = %9.3g Msol"%(TotalEjectedMass)
for i,label in enumerate(labels):
print " Ejected mass in %07s = %9.3g Msol"%(label,TotalEjectedEltMass[i])
print
for i,label in enumerate(labels):
XH = log10( (TotalEjectedEltMass[i]/(opt.Mgas+TotalEjectedMass)) / SolarAbun[label] )
print " [%07s/H] = %9.3g"%(label,XH)
print
##############################
# Ejection from SNIa
##############################
EjectedMass = chemistry.SNIa_mass_ejection(m1,m2)
TotalEjectedMass = EjectedMass[0]*opt.M0
TotalEjectedEltMass = EjectedMass[2:]*opt.M0
print "Total ejected mass (SNIa) = %9.3g Msol"%(TotalEjectedMass)
for i,label in enumerate(labels):
print " Ejected mass in %07s = %9.3g Msol"%(label,TotalEjectedEltMass[i])
print
for i,label in enumerate(labels):
XH = log10( (TotalEjectedEltMass[i]/(opt.Mgas+TotalEjectedMass)) / SolarAbun[label] )
print " [%07s/H] = %9.3g"%(label,XH)
print
Event Timeline
Log In to Comment