Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F110057111
gplotcyl_r-val
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
Thu, Apr 24, 12:45
Size
5 KB
Mime Type
text/x-python
Expires
Sat, Apr 26, 12:45 (2 d)
Engine
blob
Format
Raw Data
Handle
25781529
Attached To
rGTOOLS Gtools
gplotcyl_r-val
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 = 1.
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 = 2,
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")
(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
momentum = options.momentum
val = options.val
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)
#g.lweight(1.0)
#g.expand(1.0)
# axes properties
#labelfont = 16
#ax = pt.gca()
#pt.setp(ax.get_xticklabels(),fontsize=labelfont)
#pt.setp(ax.get_yticklabels(),fontsize=labelfont)
pt.InitPlot(files,options)
#######################################
# LOOP
#######################################
valcmd = "val=%s"%(val)
print valcmd
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
################
# 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)
x = log_filter_inv(k,0,rmax,rc)
y = vec_sigma*100*unit_velocity
if file == files[0]:
xmin,xmax,ymin,ymax = pt.SetLimits(xmin,xmax,ymin,ymax,x,y)
color = pt.pcolors[file]
pt.plot(x,y,c=color)
i = i + 1
# add legend
#pt.legend((
#r'$\textrm{stars}$',
#r'$\textrm{visible gas}$',
#r'$\textrm{dark gas}$',
#),'upper right',shadow=True)
pt.xlabel(r'$\rm{Radius}\,\left[ kpc \right]$',fontsize=pt.labelfont)
pt.ylabel(r'$\rm{Velocity}\,\left[\rm{km}/\rm{s}\right]$',fontsize=pt.labelfont)
pt.axis([xmin,xmax,ymin,ymax])
# -- end ---
if ps==None:
pt.show()
else:
pt.savefig(ps)
Event Timeline
Log In to Comment