Page MenuHomec4science
Files F108935408

genergy
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 19, 04:04

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.Colors(n=len(files))
labels = []
datas = []
# read files
for file in files:
try:
vals = io.read_new_energy(file)
except:
vals = io.read_energy(file,iobs=None)
# check if this is a swift file
if vals.has_key("E_tot"):
n_vals = {}
n_vals["Time"] = vals["Time"]
n_vals["EnergyInt"] = vals["E_int"]
n_vals["EnergyPot"] = vals["E_pot"]
n_vals["EnergyKin"] = vals["E_kin"]
vals = n_vals
x = vals['Time']
if opt.obs == 'EnergyTot':
y = 0.
nrjvals = ['EnergyInt','EnergyPot','EnergyKin','EnergyRadSph','EnergyRadSticky','EnergySfr','EnergyFeedbackWind','EnergyBubbles','EnergyAGNHeat','EnergyThermalFeedback','EnergyKineticFeedback','EnergyDissipationForces','EnergyICDissipation']
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)
data = pt.DataPoints(x,y,color=colors.get(),label=file,tpe='line')
datas.append(data)
# now, plot
for d in datas:
pt.plot(d.x,d.y,color=d.color)
# set limits and draw axis
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
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