Page MenuHomec4science
Files F102364828

gplotcyl_r-val_wc-003
No OneTemporary
Actions

File Metadata

Created
Wed, Feb 19, 23:01

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

c4science · Help