Homec4science
Diffusion Gear (master) Browse

Gear/PyChem/scripts/pychem_plot_chemivol
b87f34a0ce7emaster

Last Change
Raw File

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