Page MenuHomec4science
Files F121607956

gplot_histovals
No OneTemporary
Actions

File Metadata

Created
Sat, Jul 12, 10:12

gplot_histovals
View Options

#!/home/revaz/local/bin/python
'''
Plot histogram of a certain value of the model
Yves Revaz
ven oct 6 10:58:38 CEST 2006
'''
from numarray import *
from Nbody import *
import SM
import string
import sys
import os
import numarray.mlab as mlab
from libjeans import *
from Nbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
from Gtools import vanderwaals as vw
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_units_options(parser)
parser = add_log_options(parser)
parser = add_reduc_options(parser)
parser = add_select_options(parser)
parser.add_option("-t",
action="store",
dest="ftype",
type="string",
default = None,
help="type of the file",
metavar=" TYPE")
parser.add_option("--cgs",
action="store_true",
dest="cgs",
default = 0,
help="invert into cgs units")
parser.add_option("--nbins",
action="store",
dest="nbins",
type="int",
default=25,
help="number of bins")
parser.add_option("--cmd",
action="store",
dest="cmd",
type="string",
default = 'v = nb.T()',
help="command to extract value ex. 'v = nb.T()' ",
metavar=" STRING")
(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
reduc = options.reduc
cgs = options.cgs
ftype = options.ftype
select = options.select
nbins = options.nbins
cmd = options.cmd
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:
nb = Nbody(file,ftype=ftype)
# select
if select!=None:
print "select %s"%(select)
nb = nb.select(select)
if reduc!= None:
print "reducing %s"%(reduc)
nb = nb.reduc(reduc)
# get values
nb.localsystem_of_units = localsystem
exec(cmd)
if xmin==None:
xmin = min(v)
if xmax==None:
xmax = max(v)
x = arange(xmin,xmax,(xmax-xmin)/nbins)
y = histogram(v.astype(Float),x)[:-1]
x = x[:-1]
# print some info
print "vmin = %g vmax =%g"%(v.min(),v.max())
print "mean = %g std =%g"%(v.mean(),mlab.std(v))
print "median = %g "%(mlab.median(nb.rho))
print
print "xmin = %g xmax =%g"%(xmin,xmax)
print "ymin = %g ymax =%g"%(min(y),max(y))
print "number of bins = %d"%(len(y))
print
print "total number of particles = %d\n"%(sum(y))
# use log
if log != None:
x,y = Graph_UseLog(x,y,log)
g.ticksize(-1,0,-1,-1)
if file == files[0]:
xmin,xmax,ymin,ymax = Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y)
Graph_DrawBox(g,xmin,xmax,ymin,ymax,log)
g.ctype(colors[file])
# plot points
g.histogram(x,y)
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log %s'%cmd[4:])
g.ylabel('log N')
elif log == 'x':
g.xlabel('log %s'%cmd[4:])
g.ylabel('N')
elif log == 'y':
g.xlabel('%s'%cmd[4:])
g.ylabel('log N')
else:
g.xlabel('%s'%cmd[4:])
g.ylabel('N')
# -- end ---
Graph_End(g,ps)

Event Timeline

c4science · Help