Page MenuHomec4science
Files F91157678

genergy
No OneTemporary
Actions

File Metadata

Created
Fri, Nov 8, 11:47

genergy
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:06:11 CEST 2006
'''
import Ptools as pt
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
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("-o",
action="store",
dest="obs",
type="string",
default = 'EnergyTot',
help="observable name",
metavar=" NAME")
parser.add_option("--relative",
action="store_true",
dest="rel",
default = 0,
help="plot relative value")
(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 = []
# read files
for file in files:
try:
vals = io.read_new_energy(file)
except:
vals = io.read_energy(file,iobs=None)
x = vals['Time']
if opt.obs == 'EnergyTot':
y = 0.
nrjvals = ['EnergyInt','EnergyPot','EnergyKin','EnergyRadSph','EnergyRadSticky','EnergySfr','EnergyFeedbackWind','EnergyBubbles','EnergyAGNHeat','EnergyThermalFeedback','EnergyKineticFeedback','EnergyDissipationForces']
for nrjval in nrjvals:
if vals.has_key(nrjval):
print nrjval
y = y + vals[nrjval]
elif opt.obs == 'Virial':
y = -2*vals['EnergyKin'] / vals['EnergyPot']
elif vals.has_key(opt.obs):
y = vals[opt.obs]
else:
print "unknown observable %s"%(opt.obs)
sys.exit()
if opt.rel:
y = 100*(fabs(y-y[0]))/y[0]
y = fabs(y)
xmin,xmax,ymin,ymax = pt.SetLimits(opt.xmin,opt.xmax,opt.ymin,opt.ymax,x,y,opt.log)
# plot points
pt.plot(x,y)
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
# labels
pt.xlabel(r'$\textrm{Time}$',fontsize=pt.labelfont)
pt.ylabel(r'$\textrm{%s}$'%(opt.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)

Event Timeline

c4science · Help