Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73244863
gplotcyl_r-vcirc
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
Fri, Jul 19, 23:07
Size
7 KB
Mime Type
text/x-python
Expires
Sun, Jul 21, 23:07 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
19171831
Attached To
rGTOOLS Gtools
gplotcyl_r-vcirc
View Options
#!/usr/bin/env python
'''
Plot velocity dispertion in z as a function of the radius
Yves Revaz
Mon Jun 4 16:05:18 CEST 2007
'''
from numarray import *
from pNbody import *
from pNbody import cosmo
from pNbody.libutil import *
from pNbody import libdisk
import SM
import string
import sys
import os
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_gas_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("--rc",
action="store",
dest="rc",
type="float",
default = 7.,
help="critial radius for filter",
metavar=" FLOAT")
parser.add_option("--nx",
action="store",
dest="nx",
type="int",
default = 32,
help="number of bins in x",
metavar=" INT")
parser.add_option("--ny",
action="store",
dest="ny",
type="int",
default = 64,
help="number of bins in y",
metavar=" INT")
parser.add_option("--nmin",
action="store",
dest="nmin",
type="float",
default = 32,
help="min number of particles in a cell to accept value",
metavar=" INT")
parser.add_option("--eps",
action="store",
dest="eps",
type="float",
default = 0.28,
help="smoothing length",
metavar=" FLOAT")
parser.add_option("--notree",
action="store_true",
dest="notree",
default = False,
help="do not use tree to compute potential or forces")
parser.add_option("--useacc",
action="store_true",
dest="useacc",
default = False,
help="use acceleration instead of potential to compute velocity")
parser.add_option("--multicomp",
action="store_true",
dest="multicomp",
default = False,
help="decompose model into components")
parser.add_option("--write",
action="store_true",
dest="write",
default = False,
help="write output")
parser.add_option("--ref",
action="store",
dest="ref",
default = None,
type="string",
help="reference rotation curve")
(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
ftype = options.ftype
gamma = options.gamma
mu = options.mu
av = options.av
bv = options.bv
center = options.center
select = options.select
cmd = options.cmd
display = options.display
rmax = options.rmax
rc = options.rc
nx = options.nx
ny = options.ny
nmin = options.nmin
eps = options.eps
notree = options.notree
useacc = options.useacc
multicomp = options.multicomp
write = options.write
ref = options.ref
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:
# 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
################
nr = nx
nt = ny
G = 1.0
# compute radius
R = arange(nr)*rmax/float(nr)
R = R + (R[1]-R[0])/2
# here, we loop on components
ncomp = len(nb.npart)
comps = [None] + range(ncomp)
for comp in comps:
if comp != None:
nb_sub = nb.select(comp)
lcomp = comp+1
else:
nb_sub = nb
lcomp = 0
if nb_sub.nbody == 0:
continue
print "comp = ",comp
if useacc:
# compute acceleration
Accx,Accy,Accz = nb_sub.getAccelerationInCylindricalGrid(eps=eps,z=0,Rmax=rmax,nr=nr,nt=nt,UseTree=not notree)
Ar = sqrt(Accx**2+Accy**2)
# mean azimuthal values
Ar = sum(Ar)/nt
# circular velocity
Vc = sqrt(Ar*R)
else:
# compute potential in a spherical grid
Phi = nb_sub.getPotentialInCylindricalGrid(eps,z=0,Rmax=rmax,nr=nr,nt=nt,UseTree=not notree)
# mean azimuthal potential
Phi = sum(Phi)/nt
# compute frequencies
dPhi = libdisk.Diff(Phi,R)
Vc = libdisk.Vcirc(R,dPhi)
x = R.astype(Float32)
y = Vc.astype(Float32)
# write output
if write:
io.write_crv('%s.crv'%(comp),x,y)
# use log
if log != None:
x,y = Graph_UseLog(x,y,log)
g.ticksize(-1,0,-1,-1)
if (file == files[0]) and (comp==None):
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])
g.ltype(lcomp)
# plot points
g.connect(x,y)
if ref!= None:
x,y=io.read_crv(ref)
# use log
if log != None:
x,y = Graph_UseLog(x,y,log)
g.ticksize(-1,0,-1,-1)
g.ltype(0)
g.ctype(192)
g.connect(x,y)
g.ctype(0)
g.ctype(0)
if log == 'xy' or log == 'yx':
g.xlabel('log Radius')
g.ylabel('log Circular Velocity')
elif log == 'x':
g.xlabel('log Radius')
g.ylabel('Circular Velocity')
elif log == 'y':
g.xlabel('Radius')
g.ylabel('log Circular Velocity')
else:
g.xlabel('Radius')
g.ylabel('Circular Velocity')
# -- end ---
Graph_End(g,ps)
Event Timeline
Log In to Comment