Homec4science
Diffusion Gtools (master) Browse

gtools/old.scripts/gplot-sph_r-val
b0423933b748master

Last Change
Raw File

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)
c4science · Help