Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_mass_ejection
73ca2cd44116master

Last Change
Raw File

pychem_mass_ejection
View Options

#!/usr/bin/env python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
import string
from scipy import optimize
from PyChem import chemistry
UnitLength_in_cm = 3.085e+21
UnitMass_in_g = 1.989e+43
UnitVelocity_in_cm_per_s = 20725573.785998672
UnitTime_in_s = 148849920000000.0
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_reduc_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_select_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_display_options(parser)
parser = pt.add_info_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser.add_option("--x",
action="store",
dest="x",
type="string",
default = 'r',
help="x value to plot",
metavar=" STRING")
parser.add_option("--y",
action="store",
dest="y",
type="string",
default = 'T',
help="y value to plot",
metavar=" STRING")
parser.add_option("--z",
action="store",
dest="z",
type="string",
default = None,
help="z value to plot",
metavar=" STRING")
parser.add_option("--legend",
action="store_true",
dest="legend",
default = False,
help="add a legend")
parser.add_option("--elts",
action="store",
type="string",
dest="elts",
default = None,
help="element name")
parser.add_option("--m1",
action="store",
type="float",
dest="m1",
default = None,
help=" min mass in Msol")
parser.add_option("--m2",
action="store",
type="float",
dest="m2",
default = None,
help=" max mass in Msol")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
########################################################################
# MAIN
########################################################################
files,opt = parse_options()
file = files[0]
chemistry.init_chimie(file)
mmax = chemistry.get_Mmax()
mmin = chemistry.get_Mmin()
if opt.m1 != None:
m1 = opt.m1
m1 = m1 * (1.989e+33/UnitMass_in_g)
else:
m1 = mmin
if opt.m2 != None:
m2 = opt.m2
m2 = m2 * (1.989e+33/UnitMass_in_g)
else:
m2 = mmax
M0 = 1.
metals = zeros(chemistry.get_nelts(),float)
print "mmin = %f"%(mmin/ (1.989e+33/UnitMass_in_g))
print "mmax = %f"%(mmax/ (1.989e+33/UnitMass_in_g))
# Contribution of SNII
EjectedMass_SNII = chemistry.SNII_mass_ejection(m1,m2)[0]
print "Ejected mass fraction due to SNII = %f"%EjectedMass_SNII
# Contribution of SNIa
NSNIa = chemistry.SNIa_rate(m1,m2)*M0;
Mco = chemistry.get_Mco() # in code unit
EjectedMass_SNIa = Mco * NSNIa;
print "Ejected mass fraction due to SNIa = %f"%EjectedMass_SNIa
# Total from chemistry
EjectedMass = chemistry.Total_mass_ejection(m1,m2,M0,metals)
print "Total Ejected mass fraction = %f"%EjectedMass[0]
print "Total Ejected mass fraction = %f"%(EjectedMass_SNIa+EjectedMass_SNII)
c4science · Help