Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F100631338
gplotcyl_r-val_wc-003
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
Sat, Feb 1, 07:59
Size
8 KB
Mime Type
text/x-python
Expires
Mon, Feb 3, 07:59 (2 d)
Engine
blob
Format
Raw Data
Handle
24002915
Attached To
rGTOOLS Gtools
gplotcyl_r-val_wc-003
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 *
import Ptools as pt
import string
import sys
import os
from optparse import OptionParser
from Gtools import *
from Gtools import vanderwaals as vw
import copy
unit_velocity = 0.879
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("--momentum",
action="store",
dest="momentum",
type="int",
default = 0,
help="momentum to plot 0,1,2 ",
metavar=" INT")
parser.add_option("--val",
action="store",
dest="val",
type="string",
default = "nb.mass",
help="value to plot : nb.Vr() ",
metavar=" STRING")
parser.add_option("--add_legend",
action="store_true",
dest="add_legend",
default = False,
help="add legend")
(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
#######################################
# MakePlot
#######################################
def MakePlot(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
momentum = options.momentum
val = options.val
#######################################
# init graph
#######################################
labelfont = 16
ax = pt.gca()
pt.setp(ax.get_xticklabels(),fontsize=labelfont)
pt.setp(ax.get_yticklabels(),fontsize=labelfont)
palette = pt.GetPalette()
colors = pt.SetColorsForFiles(files,palette)
#######################################
# LOOP
#######################################
valcmd = "val=%s"%(val)
print valcmd
crv_lst = []
i = 0
# 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
################
nbs = nb
################
# warm gas
nb = nbs.select('gas')
nb = nb.selectc((nbs.u>0))
# extract x and y and normalize
# radius
r = log_filter(nb.rxy(),0,rmax,rc)
# azimuth
t = (nb.phi_xy()+pi)/(2*pi)
nr = nx
nt = ny
exec(valcmd)
# compute the value
vec_sigma = Extract1dMeanFrom2dMap(r,t,nb.mass,val,nr,nt,nmin,momentum)
# cmopute x and y
k = (arange(nx)/float(nx)).astype(float)
xw = log_filter_inv(k,0,rmax,rc)
yw = vec_sigma
################
# cold gas
nb = nbs.select('gas')
nb = nb.selectc((nbs.u<0))
# extract x and y and normalize
# radius
r = log_filter(nb.rxy(),0,rmax,rc)
# azimuth
t = (nb.phi_xy()+pi)/(2*pi)
nr = nx
nt = ny
exec(valcmd)
# compute the value
vec_sigma = Extract1dMeanFrom2dMap(r,t,nb.mass,val,nr,nt,nmin,momentum)
# cmopute x and y
k = (arange(nx)/float(nx)).astype(float)
xc = log_filter_inv(k,0,rmax,rc)
yc = vec_sigma
################
# stars
nb = nbs.select('disk')
# extract x and y and normalize
# radius
r = log_filter(nb.rxy(),0,rmax,rc)
# azimuth
t = (nb.phi_xy()+pi)/(2*pi)
nr = nx
nt = ny
exec(valcmd)
# compute the value
nmin = 5
vec_sigma = Extract1dMeanFrom2dMap(r,t,nb.mass,val,nr,nt,nmin,momentum)
# cmopute x and y
k = (arange(nx)/float(nx)).astype(float)
xs = log_filter_inv(k,0,rmax,rc)
ys = vec_sigma
yw = yw*100*unit_velocity
yc = yc*100*unit_velocity
ys = ys*100*unit_velocity
xw = compress(yw>0,xw)
yw = compress(yw>0,yw)
xc = compress(yc>0,xc)
yc = compress(yc>0,yc)
xs = compress(ys>0,xs)
ys = compress(ys>0,ys)
x = concatenate((xw,xc,xs))
y = concatenate((yw,yc,ys))
if file == files[0]:
xmin,xmax,ymin,ymax = pt.SetLimits(options.xmin,options.xmax,options.ymin,options.ymax,x,y,options.log)
'''
#g.ctype(colors[file])
# plot points
g.ltype(1-i)
g.ctype(192)
g.connect(xw,yw)
g.ctype(64)
if not i==0:
g.connect(xc,yc)
g.ltype(0)
'''
if i==0:
ls = ':'
else:
ls = '-'
crv_lst.append( pt.plot(xw,yw,'r',ls=ls) )
if not i==0:
crv_lst.append( pt.plot(xc,yc,'b',ls=ls) )
i = i + 1
pt.SetAxis(xmin,xmax,ymin,ymax,options.log)
pt.xlabel(r'$\rm{Radius}\,\left[ kpc \right]$',fontsize=labelfont)
pt.ylabel(r'$\rm{Radial\,\,Velocity\,\,Dispersion\,\,height}\,\left[ km/s \right]$',fontsize=labelfont)
# add legend
if options.add_legend:
labels = [r'$\rm{N=1 : visible\,\,gas}$',
r'$\rm{N=3 : visible\,\,gas}$',
r'$\rm{N=3 : dark\,\,gas}$']
pt.legend(crv_lst,labels,'upper right',shadow=True)
####################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
#pt.figure(figsize=(8*2,6*2))
#pt.figure(dpi=10)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
Event Timeline
Log In to Comment