Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_get_IMF_info_and_ejectats
d89b118970f3master

Last Change
Raw File

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
c4science · Help