Homec4science
Diffusion Gtools (master) Browse

gtools/scripts/gplot_energy
8cfdf441e0b6master

Last Change
Raw File

gplot_energy
View Options

#!/usr/bin/env python
'''
Extract and plot energy and mass values contained in the
output Gadget file called by default "energy.txt".
Yves Revaz
ven jun 9 10:43:59 CEST 2006
'''
from numpy import *
from pNbody import *
import string
import sys
import os
from pNbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
from Gtools import io
import Ptools as pt
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser = pt.add_cmd_options(parser)
parser.add_option("--legend",
action="store_true",
dest="legend",
default = False,
help="add a legend")
(options, args) = parser.parse_args()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
palette = pt.GetPalette()
colors = pt.SetColorsForFiles(files,palette)
labels = []
#######################################
# LOOP
#######################################
# set ltype
ltype = [0,1,2,3,4,5,6]
i = 0
# read files
for file in files:
try:
vals = io.read_new_energy(file)
except "NotNewEgyFileError":
vals = io.read_energy(file,iobs=None)
x = vals['Time']
y = 0
xs = array([],float)
ys = array([],float)
EnergyInt = None
EnergyPot = None
EnergyKin = None
EnergyRadSph = None
EnergyRadSticky = None
EnergySfr = None
EnergyFeedbackWind = None
EnergyBubbles = None
EnergyAGNHeat = None
EnergyThermalFeedback = None
EnergyKineticFeedback = None
if vals.has_key('EnergyInt'):
EnergyInt = vals['EnergyInt']
y = y + EnergyInt
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyInt))
if vals.has_key('EnergyPot'):
EnergyPot = vals['EnergyPot']
y = y + EnergyPot
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyPot))
if vals.has_key('EnergyKin'):
EnergyKin = vals['EnergyKin']
y = y + EnergyKin
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyKin))
if vals.has_key('EnergyRadSph'):
EnergyRadSph = vals['EnergyRadSph']
y = y + EnergyRadSph
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyRadSph))
if vals.has_key('EnergyRadSticky'):
EnergyRadSticky = vals['EnergyRadSticky']
y = y + EnergyRadSticky
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyRadSticky))
if vals.has_key('EnergySfr'):
EnergySfr = vals['EnergySfr']
y = y + EnergySfr
xs = concatenate((xs,x))
ys = concatenate((ys,EnergySfr))
if vals.has_key('EnergyFeedbackWind'):
EnergyFeedbackWind = vals['EnergyFeedbackWind']
y = y + EnergyFeedbackWind
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyFeedbackWind))
if vals.has_key('EnergyBubbles'):
EnergyBubbles = vals['EnergyBubbles']
y = y + EnergyBubbles
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyBubbles))
if vals.has_key('EnergyAGNHeat'):
EnergyAGNHeat = vals['EnergyAGNHeat']
y = y + EnergyAGNHeat
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyAGNHeat))
if vals.has_key('EnergyThermalFeedback'):
EnergyThermalFeedback = vals['EnergyThermalFeedback']
y = y + EnergyThermalFeedback
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyThermalFeedback))
if vals.has_key('EnergyKineticFeedback'):
EnergyKineticFeedback = vals['EnergyKineticFeedback']
y = y + EnergyKineticFeedback
xs = concatenate((xs,x))
ys = concatenate((ys,EnergyKineticFeedback))
xmin,xmax,ymin,ymax = pt.SetLimits(opt.xmin,opt.xmax,opt.ymin,opt.ymax,x,y,opt.log)
# plot points
t = x
curves = []
labels = []
color = pt.SetColor(0,palette)
curves.append(pt.plot(t,y,c=color,label='EnergyTot'))
labels.append('EnergyTot')
if EnergyInt != None:
color = pt.SetColor(160,palette)
curves.append(pt.plot(t,EnergyInt,c=color,label='EnergyInt'))
labels.append('EnergyInt')
print 'EnergyInt'
if EnergyPot != None:
color = pt.SetColor(192,palette)
curves.append(pt.plot(t,EnergyPot,c=color,label='EnergyPot'))
labels.append('EnergyPot')
print 'EnergyPot'
if EnergyKin != None:
color = pt.SetColor(128,palette)
curves.append(pt.plot(t,EnergyKin,c=color,label='EnergyKin'))
labels.append('EnergyKin')
print 'EnergyKin'
if EnergyRadSph != None:
color = pt.SetColor(64,palette)
curves.append(pt.plot(t,EnergyRadSph,c=color,label='EnergyRadSph'))
labels.append('EnergyRadSph')
print 'EnergyRadSph'
if EnergyRadSticky != None:
color = pt.SetColor(96,palette)
curves.append(pt.plot(t,EnergyRadSticky,c=color,label='EnergyRadSticky'))
labels.append('EnergyRadSticky')
print 'EnergyRadSticky'
if EnergySfr != None:
color = pt.SetColor(132,palette)
curves.append(pt.plot(t,EnergySfr,c=color,label='EnergySfr'))
labels.append('EnergySfr')
print 'EnergySfr'
if EnergyFeedbackWind != None:
color = pt.SetColor(250,palette)
curves.append(pt.plot(t,EnergyFeedbackWind,c=color,label='EnergyFeedbackWind'))
labels.append('EnergyFeedbackWind')
print 'EnergyFeedbackWind'
if EnergyBubbles != None:
color = pt.SetColor(175,palette)
curves.append(pt.plot(t,EnergyBubbles,c=color,label='EnergyBubbles'))
labels.append('EnergyBubbles')
print 'EnergyBubbles'
if EnergyAGNHeat != None:
color = pt.SetColor(168,palette)
curves.append(pt.plot(t,EnergyAGNHeat,c=color,label='EnergyAGNHeat'))
labels.append('EnergyAGNHeat')
print 'EnergyAGNHeat'
if EnergyThermalFeedback != None:
color = pt.SetColor(175,palette)
curves.append(pt.plot(t,EnergyThermalFeedback,c=color,label='EnergyThermalFeedback'))
labels.append('EnergyThermalFeedback')
print 'EnergyThermalFeedback'
if EnergyKineticFeedback != None:
color = pt.SetColor(168,palette)
curves.append(pt.plot(t,EnergyKineticFeedback,c=color,label='EnergyKineticFeedback'))
labels.append('EnergyKineticFeedback')
print 'EnergyKineticFeedback'
i=i+1
# labels
obs = "Energy"
#pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
if opt.legend:
pt.legend( labels , loc=3)
if log == 'xy' or log == 'yx':
pt.xlabel(r'$\textrm{log Time}$',fontsize=pt.labelfont)
pt.ylabel(r'$\textrm{log %s}$'%(obs),fontsize=pt.labelfont)
elif log == 'x':
pt.xlabel(r'$\textrm{log Time}$',fontsize=pt.labelfont)
pt.ylabel(r'$\textrm{%s}$'%(obs),fontsize=pt.labelfont)
elif log == 'y':
pt.xlabel(r'$\textrm{Time}$',fontsize=pt.labelfont)
pt.ylabel(r'$\textrm{log %s}$'%(obs),fontsize=pt.labelfont)
else:
pt.xlabel(r'$\textrm{Time}$',fontsize=pt.labelfont)
pt.ylabel(r'$\textrm{%s}$'%(obs),fontsize=pt.labelfont)
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
#pt.figure(figsize=(8*2,6*2))
#pt.figure(dpi=10)
pt.pcolors
#fig = pt.gcf()
#fig.subplots_adjust(left=0.1)
#fig.subplots_adjust(right=1)
#fig.subplots_adjust(bottom=0.12)
#fig.subplots_adjust(top=0.95)
#fig.subplots_adjust(wspace=0.25)
#fig.subplots_adjust(hspace=0.02)
MakePlot(files,opt)
pt.EndPlot(files,opt)
c4science · Help