Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F91156430
gplotcyl_r-val_svd
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, Nov 8, 11:30
Size
6 KB
Mime Type
text/x-python
Expires
Sun, Nov 10, 11:30 (2 d)
Engine
blob
Format
Raw Data
Handle
22205593
Attached To
rGTOOLS Gtools
gplotcyl_r-val_svd
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)
#######################################
# 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
################
nbs = nb
################
# warm gas
nb = nbs.select('wg')
# 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('cg')
# 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)
# compute 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','stars')
# 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)
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
x = concatenate((xw,xc,xs))
y = concatenate((yw,yc,ys))
if file == files[0]:
xmin,xmax,ymin,ymax = pt.SetLimits(xmin,xmax,ymin,ymax,x,y)
pt.plot(xs,ys,'k')
pt.plot(xw,yw,'r')
pt.plot(xc,yc,'b')
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=labelfont)
pt.ylabel(r'$\rm{Velocity}\,\left[\rm{km}/\rm{s}\right]$',fontsize=labelfont)
pt.axis([xmin,xmax,ymin,ymax])
# -- end ---
if ps==None:
pt.show()
else:
pt.savefig(ps)
Event Timeline
Log In to Comment