Homec4science
Diffusion Gtools (master) Browse

gtools/dev/graphs_london/model/gplotcyl_r-val_sr
8cfdf441e0b6master

Last Change
Raw File

gplotcyl_r-val_sr
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 = 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")
(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
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
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)
c4science · Help