Page MenuHomec4science
Files F93067870

gplot_r-egrav
No OneTemporary
Actions

File Metadata

Created
Mon, Nov 25, 23:14

gplot_r-egrav
View Options

#!/usr/bin/env python
'''
Plot density Energy of the model as a function of radius
(isotropique density)
Yves Revaz
Thu Feb 23 15:00:11 CET 2006
'''
from numarray import *
from Nbody import *
import SM
import string
import sys
import os
from libjeans import *
from Nbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
def massr(nb,rmax,nbins):
rs = arange(0,rmax,rmax/nbins)
rs = rs[1:]
mr = array([],Float)
r = nb.rxyz()
for rp in rs:
c = (r<rp)
m = sum(compress(c,nb.mass))
mr = concatenate((mr,m))
return rs,mr
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_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("--ev",
action="store_true",
dest="ev",
default = 0,
help="give energy in eV")
parser.add_option("--rmax",
action="store",
dest="rmax",
type="float",
default=10.,
help="max radius")
parser.add_option("--nbins",
action="store",
dest="nbins",
type="int",
default=100,
help="number of bins")
(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
ev = options.ev
ftype = options.ftype
rmax = options.rmax
nbins = options.nbins
if cgs:
units = get_units_options(options)
#######################################
# open sm
#######################################
g = Graph_Init(ps)
Graph_SetDefaultsGraphSettings(g)
colors = Graph_SetColorsForFiles(files,col)
#######################################
# LOOP
#######################################
# read files
for file in files:
nbody = Nbody(file,ftype=ftype)
nbody.hdcenter()
if reduc!= None:
nbody = nbody.reduc(reduc)
if cgs:
nbody.pos = Length2cgs(nbody.pos,units)
x,y = massr(nbody,rmax,nbins)
# convert into cgs
if cgs:
#x = Length2cgs(x,units)
y = Mass2cgs(y,units)*GRAVITY/x
if ev:
y = y*PROTONMASS*1e-7/1.602e-19
# use log
if log != None:
x,y = Graph_UseLog(x,y,log)
if file == files[0]:
Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y)
g.box()
# plot points
g.ctype(colors[file])
g.connect(x,y)
g.ctype(0)
g.xlabel('R')
if not ev:
if log == 'xy' or log == 'yx':
g.xlabel('log R')
g.ylabel('log Egrav Spec')
elif log == 'x':
g.xlabel('log R')
g.ylabel('Egrav Spec')
elif log == 'y':
g.xlabel('R')
g.ylabel('log Egrav Spec')
else:
g.xlabel('R')
g.ylabel('Egrav Spec')
else:
if log == 'xy' or log == 'yx':
g.xlabel('log R')
g.ylabel('log Egrav [eV/part]')
elif log == 'x':
g.xlabel('log R')
g.ylabel('Egrav [eV/part]')
elif log == 'y':
g.xlabel('R')
g.ylabel('log Egrav [eV/part]')
else:
g.xlabel('R')
g.ylabel('Egrav [eV/part]')
# -- end ---
Graph_End(g,ps)

Event Timeline

c4science · Help