Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F115194259
pychem_plot_chemivol
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, May 30, 10:25
Size
9 KB
Mime Type
text/x-python
Expires
Sun, Jun 1, 10:25 (2 d)
Engine
blob
Format
Raw Data
Handle
26439238
Attached To
rGEAR Gear
pychem_plot_chemivol
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
import pickle
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)
parse = pt.add_units_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")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
########################################################################
# M A I N
########################################################################
def MakePlot(files,opt):
file = files[0]
f = open(file)
optf = pickle.load(f)
gear_units = pickle.load(f)
chimiefile = pickle.load(f)
IN = pickle.load(f)
f.close()
Times = IN["Times"]
TotStarMass = IN["TotStarMass"]
TotGasMass = IN["TotGasMass"]
TotMetals = IN["TotMetals"]
Sfr = IN["Sfr"]
Sfr_th = IN["Sfr_th"]
TotSNII = IN["TotSNII"]
TotSNIa = IN["TotSNIa"]
TotDYIN = IN["TotDYIN"]
NStars = IN["NStars"]
##############
# init chimie
##############
chemistry.init_chimie(chimiefile,optf.NumberOfTables,optf.DefaultTable)
elts = chemistry.get_elts_labels()
nelts = chemistry.get_nelts()
##############
# units
##############
# units sfr (Msol/yr)
sfr_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_yr,units.Unit_K])
# units time (Myr)
time_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_Myr,units.Unit_K])
# some conversion factors
UnitsConvEnergy_cgs2gear = 1.0/gear_units.convertionFactorTo(units.cgs.UnitEnergy)
UnitsConvSfr_sfr2gear = 1.0/( gear_units.convertionFactorTo(sfr_units.UnitMass) / gear_units.convertionFactorTo(sfr_units.UnitTime) )
UnitsConvSfr_gear2sfr = 1.0/UnitsConvSfr_sfr2gear
UnitsConvMass_Msol2gear = 1.0/gear_units.convertionFactorTo(sfr_units.UnitMass)
UnitsConvMass_gear2Msol = 1/UnitsConvMass_Msol2gear
UnitsConvTime_Myrs2gear = 1.0/gear_units.convertionFactorTo(time_units.UnitTime)
UnitsConvTime_gear2Myrs = 1/UnitsConvTime_Myrs2gear
SolarAbun = chemistry.get_elts_SolarMassAbundances()
TIMES = Times*UnitsConvTime_gear2Myrs
SFR = Sfr * UnitsConvSfr_gear2sfr
SFRTH = Sfr_th * UnitsConvSfr_gear2sfr
TOTSTARMASS = TotStarMass * UnitsConvMass_gear2Msol
TOTGASMASS = TotGasMass * UnitsConvMass_gear2Msol
elt = "Fe"
FeH = log10(TotMetals[:,elts.index(elt)] / SolarAbun[elt] + 1.0e-10)
elt = "Ba"
BaH = log10(TotMetals[:,elts.index(elt)] / SolarAbun[elt] + 1.0e-10)
elt = "Mg"
MgH = log10(TotMetals[:,elts.index(elt)] / SolarAbun[elt] + 1.0e-10)
np = 6
i = 0
#############################
# Stellar particle vs time
#############################
i+=1
xlabel = r'$\rm{Time} [\rm{Myrs}]$'
ylabel = r"$\rm{Nstars}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,NStars,tpe='line',label='sfr',color='k')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
#############################
# SFR vs time
#############################
i+=1
xlabel = r'$\rm{Time} [\rm{Myrs}]$'
ylabel = r"$\rm{SFR} [\rm{M_{\odot}/yr}]$"
pt.subplot(np,1,i)
datas = []
#data = pt.DataPoints(TIMES,SFR,tpe='line',label='sfr',color='k')
#datas.append(data)
data = pt.DataPoints(TIMES,SFRTH,tpe='line',label='sfr (imposed)',color='r')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
#############################
# SNs+DYIN vs time
#############################
i+=1
xlabel = r'$\rm{Time} [\rm{Myrs}]$'
ylabel = r"$\rm{Number} $"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,add.accumulate(TotSNII),tpe='line',label='SNII',color='r')
datas.append(data)
data = pt.DataPoints(TIMES,add.accumulate(TotSNIa),tpe='line',label='SNIa',color='b')
datas.append(data)
data = pt.DataPoints(TIMES,add.accumulate(TotDYIN),tpe='line',label='AGB',color='k')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color,label=d.label)
if opt.log==None:
log = "y"
else:
log = opt.log + 'y'
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.legend()
#############################
# M vs time
#############################
i+=1
xlabel = r'$\rm{Time} [\rm{Myrs}]$'
ylabel = r"$\rm{Mass} [\rm{M_{\odot}}]$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,TOTGASMASS,tpe='line',label='gas mass',color='b')
datas.append(data)
data = pt.DataPoints(TIMES,TOTSTARMASS,tpe='line',label='stars mass',color='r')
datas.append(data)
if opt.log==None:
log = "y"
else:
log = opt.log + 'y'
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color,label=d.label)
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.legend()
#############################
# [X/H] vs time
#############################
i+=1
xlabel = r'$\rm{Time} [\rm{Myrs}]$'
ylabel = r"$\rm{[X]}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(TIMES,FeH,tpe='line',label='Fe/H',color='r')
datas.append(data)
data = pt.DataPoints(TIMES,BaH,tpe='line',label='Ba/H',color='g')
datas.append(data)
data = pt.DataPoints(TIMES,MgH,tpe='line',label='Mg/H',color='c')
datas.append(data)
data = pt.DataPoints(TIMES,BaH-FeH,tpe='line',label='Ba/Fe',color='k')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color,label=d.label)
ymin = -5
ymax = 0
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,ymin,ymax,datas,opt.log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.legend()
#############################
# [X/H] vs [Fe/H]
#############################
i+=1
xlabel = r"$\rm{[Fe/H]}$"
ylabel = r"$\rm{[X/Fe]}$"
pt.subplot(np,1,i)
datas = []
data = pt.DataPoints(FeH,BaH-FeH,tpe='line',label='Ba',color='g')
datas.append(data)
data = pt.DataPoints(FeH,MgH-FeH,tpe='line',label='Mg',color='c')
datas.append(data)
for d in datas:
if d.tpe=='line' or d.tpe=='both':
pt.plot(d.x,d.y,color=d.color,label=d.label)
log = None
xmin = -5
xmax = 1
ymin = -3
ymax = 1
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(xmin,xmax,ymin,ymax,datas,log)
# plot
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.legend(loc="upper right")
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
Event Timeline
Log In to Comment