Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_plot_yieldsAbundances_withZ
b9a4fa236396master

Last Change
Raw File

pychem_plot_yieldsAbundances_withZ
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("--NSNII",
action="store",
dest="NSNII",
type="float",
default = 1,
help="number of SNII",
metavar=" FLOAT")
parser.add_option("--NSNIa",
action="store",
dest="NSNIa",
type="float",
default = 1,
help="number of SNIa",
metavar=" FLOAT")
parser.add_option("--NDYIN",
action="store",
dest="NDYIN",
type="float",
default = 1,
help="number of NDYIN",
metavar=" FLOAT")
parser.add_option("-Z","--Z",
action="store",
dest="Z",
type="float",
default = 0.02,
help="metallicity",
metavar=" FLOAT")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
datas = []
linest = pt.LineStyles()
if opt.elts !=None:
opt.elts = string.split(opt.elts,',')
# read files
for file in files:
chemistry.init_chimie(file)
nelts = chemistry.get_allnelts()
labels = chemistry.get_allelts_labels()
print labels
mmax = chemistry.get_Mmax()
mmin = chemistry.get_Mmin()
# mass range
ms = arange(mmin,mmax,1e-12).astype(float32)
metals = zeros(chemistry.get_nelts(),float)
metals[-1]= opt.Z
SolarAbun = chemistry.get_elts_SolarMassAbundances()
################
# get values
################
if opt.elts==None:
elts = labels
elts.remove('Ej')
elts.remove('Ejnp')
elts.remove('Metals')
else:
elts = opt.elts
################################################
# SNII (no Z dependence)
################################################
colors = pt.Colors(n=len(elts))
for elt in elts:
nelt = labels.index(elt)
nFe = labels.index('Fe')
print "(SNII)",elt,nelt
x = ms
y = zeros(len(x))
for i,m in enumerate(ms):
EjectedMass = chemistry.SNII_single_mass_ejection(m)
EjectedGasMass = EjectedMass[0]
EjectedEltMass = EjectedMass[nelt]
EjectedFeMass = EjectedMass[nFe]
if EjectedGasMass>0:
y[i] = log10( (EjectedEltMass/EjectedFeMass) / (SolarAbun[elt]/SolarAbun['Fe']) + 1.0e-100)
else:
y[i] = -100
# do some cleaning
c = y>-100
x = compress(c,x)
y = compress(c,y)
x = x / (1.989e+33/UnitMass_in_g) # to solar mass
xlabel = r"$\rm{Star\,\,Mass [M_{\odot}]}$"
ylabel = r"$\rm{[X/Fe]}$"
#label = "%s (%s)"%(labels[nelt],file)
label = "%s"%(labels[nelt])
data = pt.DataPoints(x,y,label=label,color=colors.get(),tpe='line',linestyle=linest.current())
datas.append(data)
linest.next()
# now, plot
for d in datas:
if d.tpe=='points' or d.tpe=='both':
pt.scatter(d.x,d.y,c=d.color,s=5,linewidths=0,marker='o',vmin=opt.zmin,vmax=opt.zmax)
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color,ls=d.linestyle)
pt.text(d.x[-1]+5,d.y[-1],d.label,color=d.color, horizontalalignment='left', verticalalignment='center')
if opt.legend:
pt.LegendFromDataPoints(datas,loc='lower right')
################################################
# SNIa (no Z dependence, no mass dependence)
################################################
colors = pt.Colors(n=len(elts))
for elt in elts:
nelt = labels.index(elt)
nFe = labels.index('Fe')
label = elt
print "(SNIa)",elt,nelt
EjectedMass = chemistry.SNIa_single_mass_ejection(1)
EjectedGasMass = EjectedMass[0]
EjectedEltMass = EjectedMass[nelt]
EjectedFeMass = EjectedMass[nFe]
if EjectedGasMass>0:
z = log10( (EjectedEltMass/EjectedFeMass) / (SolarAbun[elt]/SolarAbun['Fe']) + 1.0e-100)
x = (0,3)
y = (z,z)
color = colors.get()
pt.plot(x,y,color=color)
pt.text(x[0]-3,y[0],label,color=color, horizontalalignment='right', verticalalignment='center')
data = pt.DataPoints(x,y)
datas.append(data)
# set limits and draw axis
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.grid(False)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
c4science · Help