Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73029408
gplot-sph_r-val
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jul 18, 04:02
Size
7 KB
Mime Type
text/x-python
Expires
Sat, Jul 20, 04:02 (2 d)
Engine
blob
Format
Raw Data
Handle
19133894
Attached To
rGTOOLS Gtools
gplot-sph_r-val
View Options
#!/usr/bin/env python
'''
Plot value computed in a spherical grid, as a function of the radius
Yves Revaz
Wed Jan 9 09:29:46 CET 2008
'''
from pNbody import *
from pNbody import myNumeric
from pNbody import libdisk
from pNbody import libgrid
import SM
import string
import sys
import os
from optparse import OptionParser
from Gtools import *
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_center_options(parser)
parser = add_select_options(parser)
parser = add_cmd_options(parser)
parser = add_display_options(parser)
parser.add_option("-t",
action="store",
dest="ftype",
type="string",
default = None,
help="type of the file",
metavar=" TYPE")
parser.add_option("--rmax",
action="store",
dest="rmax",
type="float",
default = 50.,
help="max radius of bins",
metavar=" FLOAT")
parser.add_option("--nr",
action="store",
dest="nr",
type="int",
default = 32,
help="number of bins in r",
metavar=" INT")
parser.add_option("--eps",
action="store",
dest="eps",
type="float",
default = 0.1,
help="smoothing length",
metavar=" FLOAT")
parser.add_option("--fct",
action="store",
dest="fct",
type="string",
default = None,
help="transformation function ex : 'r**3' should be of the form M(r)",
metavar=" STR")
parser.add_option("--fctm",
action="store",
dest="fctm",
type="string",
default = None,
help="inverse transformation function ex : 'r**(1/3.)'",
metavar=" STR")
parser.add_option("--mode",
action="store",
dest="mode",
type="string",
default = None,
help="mode",
metavar=" STR")
parser.add_option("--notree",
action="store_true",
dest="notree",
default = False,
help="do not use tree to compute potential or forces")
(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
logxy = options.log
reduc = options.reduc
ftype = options.ftype
center = options.center
select = options.select
cmd = options.cmd
display = options.display
rmax = options.rmax
nr = options.nr
eps = options.eps
fct = options.fct
fctm = options.fctm
mode = options.mode
notree = options.notree
#localsystem = Set_SystemUnits_From_Options(options)
#######################################
# define output system of unit
#######################################
#outputunits = units.UnitSystem('mks',[units.Unit_km, units.Unit_Ms, units.Unit_s , units.Unit_K, units.Unit_mol, units.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)
# set model units
#nb.localsystem_of_units = localsystem
# center the model
if center=='hdcenter':
print "centering %s"%(center)
nb.hdcenter()
elif center=='histocenter':
print "centering %s"%(center)
nb.histocenter()
# select
if select!=None and select!="None":
print "select %s"%(select)
nb = nb.select(select)
# reduc
if reduc!= None:
print "reducing %s"%(reduc)
nb = nb.reduc(reduc)
# cmd
if (cmd!= None) and (cmd!= 'None'):
print nb.nbody
print "exec : %s"%cmd
exec(cmd)
# display
if (display!= None) and (display!= 'None'):
nb.display(obs=None,view=display,marker='cross')
################
# get values
################
if fct!=None:
print "f =lambda r:(%s)"%fct
exec("f =lambda r:(%s)"%fct)
else:
f = None
if fctm!=None:
print "fm =lambda r:(%s)"%fctm
exec("fm =lambda r:(%s)"%fctm)
else:
fm = None
r = libgrid.get_r_Of_Spherical_1d_Grid(nr,rmax,f=f,fm=fm)
if mode == None:
val = libgrid.get_NumberMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'nr':
val = libgrid.get_NumberMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'mr':
val = libgrid.get_MassMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'pot':
val = libgrid.get_PotentialMap_On_Spherical_1d_Grid(nb,nr,rmax,eps,f=f,fm=fm)
elif mode == 'srf':
val = libgrid.get_SurfaceMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'vol':
val = libgrid.get_VolumeMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'den':
val = libgrid.get_DensityMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'lden':
val = libgrid.get_LinearDensityMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'accm':
val = libgrid.get_AccumulatedMassMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
elif mode == 'dpot':
pot = libgrid.get_PotentialMap_On_Spherical_1d_Grid(nb,nr,rmax,eps,f=f,fm=fm)
val = libgrid.get_First_Derivative(pot,r)
elif mode == 'rgrid':
r = libgrid.get_r_Of_Spherical_1d_Grid(nr,rmax,f=None,fm=None)
val = libgrid.get_r_Of_Spherical_1d_Grid(nr,rmax,f=f,fm=fm)
elif mode == 'dr':
r1 = libgrid.get_r_Of_Spherical_1d_Grid(nr,rmax,offr=0,f=f,fm=fm)
r2 = libgrid.get_r_Of_Spherical_1d_Grid(nr,rmax,offr=1,f=f,fm=fm)
val = r2-r1
elif mode == 'dpotden':
pot = libgrid.get_PotentialMap_On_Spherical_1d_Grid(nb,nr,rmax,eps,f=f,fm=fm)
dpot = libgrid.get_First_Derivative(pot,r)
den = libgrid.get_DensityMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
val = dpot*den
elif mode == 'sigma':
phi = libgrid.get_PotentialMap_On_Spherical_1d_Grid(nb,nr,rmax,eps=eps,f=f,fm=fm)
rho = libgrid.get_DensityMap_On_Spherical_1d_Grid(nb,nr,rmax,f=f,fm=fm)
# dr
r1 = libgrid.get_r_Of_Spherical_1d_Grid(nr,rmax,offr=0,f=f,fm=fm)
r2 = libgrid.get_r_Of_Spherical_1d_Grid(nr,rmax,offr=1,f=f,fm=fm)
dr = r2-r1
val = libdisk.get_1d_Sigma_From_Rho_Phi(rho=rho,phi=phi,r=r,dr=dr)
x = r
y = val
# use log
if logxy != None:
x,y = Graph_UseLog(x,y,logxy)
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,logxy)
g.ctype(colors[file])
# plot points
g.connect(x,y)
g.ctype(0)
g.xlabel('Radius')
g.ylabel('Val')
# -- end ---
Graph_End(g,ps)
Event Timeline
Log In to Comment