Page MenuHomec4science
Files F66492439

gaccretion
No OneTemporary
Actions

File Metadata

Created
Mon, Jun 10, 21:41

gaccretion
View Options

#!/usr/bin/env python
'''
Extract and plot energy and mass values contained in the
output Gadget file called by default "accretion.txt".
Yves Revaz
mer sep 6 09:05:35 CEST 2006
'''
from numarray import *
from Nbody import *
import SM
import string
import sys
import os
from Nbody.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 = add_postscript_options(parser)
parser = add_color_options(parser)
parser = add_limits_options(parser)
parser = add_log_options(parser)
parser = add_units_options(parser)
parser.add_option("-o",
action="store",
dest="obs",
type="string",
default = 'MTot',
help="observable name",
metavar=" NAME")
parser.add_option("-r",
action="store",
dest="reduc",
type="int",
default = 1,
help="reduction number",
metavar=" INT")
(options, args) = parser.parse_args()
if options.colors!=None:
exec("options.colors = array([%s])"%(options.colors))
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#############################
# graph
#############################
# get options
files,options = parse_options()
ps = options.ps
col = options.colors
xmin = options.xmin
xmax = options.xmax
ymin = options.ymin
ymax = options.ymax
log = options.log
obs = options.obs
reduc = options.reduc
localsystem = Set_SystemUnits_From_Options(options)
#######################################
# define output system of unit
#######################################
outputunits = UnitSystem('mks',[Unit_kpc, Unit_Ms, Unit_yr , Unit_K, Unit_mol, Unit_C])
#######################################
# open sm
#######################################
g = Graph_Init(ps)
Graph_SetDefaultsGraphSettings(g)
colors = Graph_SetColorsForFiles(files,col)
#######################################
# LOOP
#######################################
# read files
for file in files:
try:
vals = io.read_accretion(file)
except:
print "problem reading %s"%(file)
sys.exit()
x = vals['Time']
if obs == 'dMdt':
y = 0.
mvals = ['M1','M1a','M1b','M2','M3','M4','M5','M6']
for mval in mvals:
if vals.has_key(mval):
print mval
y = y + vals[mval]
n = reduc
c = fmod(arange(len(x)),n)
x = compress((c==0),x)
y = compress((c==0),y)
dt = x[1:]-x[:-1]
dM = y[1:]-y[:-1]
x = x[1:]
y = dM/dt
FACT = PhysCte(1,localsystem.UnitDic['kg']/localsystem.UnitDic['s'])
FACT = FACT.into(outputunits)
y = y * FACT
elif obs == 'MTot':
y = 0.
mvals = ['M1','M1a','M1b','M2','M3','M4','M5','M6']
for mval in mvals:
if vals.has_key(mval):
print mval
y = y + vals[mval]
elif vals.has_key(obs):
y = vals[obs]
else:
print "unknown observable %s"%(obs)
sys.exit()
# use log
if log != None:
x,y = Graph_UseLog(x,y,log)
if file == files[0]:
xmin,xmax,ymin,ymax = Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y)
g.box()
# plot points
g.ctype(colors[file])
g.connect(x,y)
# labels
g.ctype(0)
g.xlabel('T')
g.ylabel('%s'%obs)
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log T')
g.ylabel('log %s'%obs)
elif log == 'x':
g.xlabel('log T')
g.ylabel('%s'%obs)
elif log == 'y':
g.xlabel('T')
g.ylabel('log %s'%obs)
else:
g.xlabel('T')
g.ylabel('%s'%obs)
# -- end ---
Graph_End(g,ps)

Event Timeline

c4science · Help