Homec4science
Diffusion Gtools (master) Browse

gtools/dev/gplotcylrQ.py
b0423933b748master

Last Change
Raw File

gplotcylrQ.py
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 *
from pNbody import libdisk
import string
import sys
import os
from optparse import OptionParser
from Gtools import *
from Gtools import vanderwaals as vw
import Ptools as pt
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("--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("--eps",
action="store",
dest="eps",
type="float",
default = 0.28,
help="smoothing length",
metavar=" FLOAT")
parser.add_option("--opt",
action="store",
dest="opt",
type="string",
default = 'gas',
help="options",
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
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
center = options.center
select = options.select
cmd = options.cmd
display = options.display
rmax = options.rmax
nx = options.nx
ny = options.ny
nmin = options.nmin
eps = options.eps
opt = options.opt
localsystem = Set_SystemUnits_From_Options(options)
#######################################
# 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
#######################################
ctypes = [0,192,64,192,64]
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)
################
# get values
################
nr = nx
nt = ny
G = 1.0
# compute radius
R = arange(nr)*rmax/float(nr)
#R = R + (R[1]-R[0])/2
# compute potential in a spherical grid
Phi = nb.getPotentialInCylindricalGrid(eps,z=0,Rmax=rmax,nr=nr,nt=nt,UseTree=True)
#------------------- do the selection only here, because before, we need all part.
# 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')
nb_sdens = nb.select('gas','disk','stars')
if opt =='wg':
nb_sigma = nb.select('gas')
nb_sigma = nb.selectc(nb.u>0)
elif opt =='cg':
nb_sigma = nb.select('gas')
nb_sigma = nb.selectc(nb.u<0)
elif opt =='gas':
nb_sigma = nb.select('gas')
elif opt =='expd':
nb_sigma = nb.select('disk')
elif opt =='wg+disk':
nb_gas = nb.select('gas')
nb_gas = nb.selectc(nb.u>0)
nb_disk = nb.select('disk','stars')
nb_sigma = nb_gas + nb_disk
elif opt =='disk':
nb_sigma = nb.select('gas','disk','stars')
# observed
nb_sdens_obs = nb.select('wg','disk','stars')
nb_sigma_obs = nb_sdens_obs
#-------------------
# compute surface density (only for gas and disk particles)
Sdens = nb_sdens.getSurfaceDensityInCylindricalGrid(Rmax=rmax,nr=nr,nt=nt)
Sdens_obs = nb_sdens_obs.getSurfaceDensityInCylindricalGrid(Rmax=rmax,nr=nr,nt=nt)
# compute radial velocity dispersion (only for gas and disk particles)
SigmaR = nb_sigma.getRadialVelocityDispersionInCylindricalGrid(Rmax=rmax,nr=nr,nt=nt)
SigmaR_obs = nb_sigma_obs.getRadialVelocityDispersionInCylindricalGrid(Rmax=rmax,nr=nr,nt=nt)
# compute number of particles per cell (only for gas and disk particles)
Num = nb_sigma.getNumberParticlesInCylindricalGrid(Rmax=rmax,nr=nr,nt=nt)
Num_obs = nb_sigma_obs.getNumberParticlesInCylindricalGrid(Rmax=rmax,nr=nr,nt=nt)
# mean azimuthal potential
Phi = sum(Phi,0)/nt
# mean azimuthal surface density
Sdens = sum(Sdens,0)/nt
Sdens_obs = sum(Sdens_obs,0)/nt
# mean azimuthal surface density
SigmaR = get1dMeanFrom2dMap(SigmaR,Num,nmin=32)
SigmaR_obs = get1dMeanFrom2dMap(SigmaR_obs,Num_obs,nmin=32)
# compute frequencies
dPhi = libdisk.Diff(Phi,R)
d2Phi = libdisk.Diff(dPhi,R)
kappa = libdisk.Kappa(R,dPhi,d2Phi)
vcirc = libdisk.Vcirc(R,dPhi)
qtoom = libdisk.QToomre(G,R,SigmaR,kappa,Sdens)
qtoom_obs = libdisk.QToomre(G,R,SigmaR_obs,kappa,Sdens_obs)
m = (R * kappa**2) / (2*pi*G*Sdens*3)
X = (R * kappa**2) / (2*pi*G*Sdens*4)
x = R
y = qtoom
y_obs = qtoom_obs
# min, max
xmin,xmax,ymin,ymax = pt.SetLimits(options.xmin,options.xmax,options.ymin,options.ymax,x,y,options.log)
# plot
pt.plot(x, y, '-', linewidth=1,c=colors[file])
pt.plot(x, y_obs, ':', linewidth=1,c=colors[file])
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{Q}$',fontsize=labelfont)
pt.grid()
# add legend
if options.add_legend:
times = [r'$0.0\rm{Gyr}$',r'$0.9\rm{Gyr}$',r'$1.8\rm{Gyr}$',r'$2.3\rm{Gyr}$',r'$3.2\rm{Gyr}$',r'$4.1\rm{Gyr}$']
pt.legend(times,'lower 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)
c4science · Help