Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F108733241
gplotcyl_r-vcirc_multiphase
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, Apr 18, 05:17
Size
7 KB
Mime Type
text/x-python
Expires
Sun, Apr 20, 05:17 (2 d)
Engine
blob
Format
Raw Data
Handle
25642406
Attached To
rGTOOLS Gtools
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 numpy import *
from pNbody import *
from pNbody import cosmo
from pNbody.libutil import *
from pNbody import libdisk
import Ptools as pt
import string
import sys
import os
from optparse import OptionParser
from Gtools import *
from Gtools import vanderwaals as vw
from Gtools import io
from scipy import interpolate
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
velocity_factor = 100.
#############################
# 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])
# axes properties
labelfont = 16
ax = pt.gca()
pt.setp(ax.get_xticklabels(),fontsize=labelfont)
pt.setp(ax.get_yticklabels(),fontsize=labelfont)
#######################################
# 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
color = {'wg':'r','cg':'b','disk':'y','bulge':'g','halo':'c','stars':'y',None:'k'}
else:
ncomp = 1
comps = [None]
lcomp = 0
color = {'wg':'r','cg':'b','disk':'y','bulge':'g','halo':'c','stars':'y',None:'k'}
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,0)/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,0)/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)) * velocity_factor
# smooth curve
n = len(x)
s1 = n-sqrt(2*n)
s2 = n+sqrt(2*n)
tck = interpolate.fitpack.splrep(x,y2,s=0.001,k=2)
yi2 = interpolate.fitpack.splev(x,tck)
yi = sqrt(yi2)
# write output
if write:
io.write_crv2('%s.crv2'%(comp),x,yi2)
# 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 = pt.SetLimits(xmin,xmax,ymin,ymax,x,y)
#g.ctype(colors[file])
#g.ltype(lcomp)
# plot points
#g.connect(x,y)
pt.plot(x,y,color[comp])
#pt.plot(x,yi,'--')
# add label
'''
pt.legend((
'total',
'visible gas',
'dark gas',
'disk',
'bulge',
'halo'
),'upper right',shadow=True)
'''
pt.axis([0,50,0,250])
pt.xlabel('Radius [kpc]',fontsize=labelfont)
pt.ylabel('Circular Velocity [km/s]',fontsize=labelfont)
# -- end ---
if ps!="None":
if ps==None:
pt.show()
else:
pt.savefig(ps)
Event Timeline
Log In to Comment