Homec4science
Diffusion Gtools (master) Browse

gtools/dev/gplotcyl_r-vcirc_multiphase
b0423933b748master

Last Change
Raw File

gplotcyl_r-vcirc_multiphase
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
from Gtools import io
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")
(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
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 is needed
if multicomp:
ncomp = len(nb.npart)
comps = [None] + ['wg','cg','disk','bulge','halo','stars']
lcomp = 0
else:
ncomp = 1
comps = [None]
lcomp = 0
for comp in comps:
if comp != None:
nb_sub = nb.select(comp)
lcomp = lcomp+1
else:
nb_sub = nb
lcomp = 0
if nb_sub.nbody == 0:
continue
nb_sub.Tree = None
#print nb_sub.nbody
if useacc:
# compute acceleration
Accx,Accy,Accz = nb_sub.getAccelerationInCylindricalGrid(eps=eps,z=0,Rmax=rmax,nr=nr,nt=nt,UseTree=True)
Ar = sqrt(Accx**2+Accy**2)
# mean azimuthal values
Ar = sum(Ar)/nt
# circular velocity
Vc2 = Ar*R
else:
# compute potential in a spherical grid
nb_sub.Tree = None
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)
Vc2 = libdisk.Vcirc2(R,dPhi)
x = R.astype(Float32)
y2 = Vc2.astype(Float32)
y = sqrt(fabs(y2))
# write output
if write:
io.write_crv2('%s.crv2'%(comp),x,y2)
# 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)
# add label
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)
c4science · Help