Homec4science
Diffusion Gtools (master) Browse

gtools/dev/graphs_london/model/gplot_r-sdens_01
8cfdf441e0b6master

Last Change
Raw File

gplot_r-sdens_01
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
import string
import sys
import os
from pNbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
import Ptools as pt
unit_mass = 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("--nbins",
action="store",
dest="nbins",
type="int",
default = 50,
help="number of bins",
metavar=" INT")
parser.add_option("--rmax",
action="store",
dest="rmax",
type="float",
default = 50.,
help="max radius of bins",
metavar=" INT")
parser.add_option("--mode",
action="store",
dest="mode",
type="string",
default = 'sr',
help="mode",
metavar=" STRING")
parser.add_option("--fSx",
action="store",
dest="fSx",
type="float",
default = 0.5,
help="fSx",
metavar=" FLOAT")
parser.add_option("--fSb",
action="store",
dest="fSb",
type="float",
default = 0.1,
help="fSb",
metavar=" FLOAT")
(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
nbins = options.nbins
rmax = options.rmax
mode = options.mode
localsystem = Set_SystemUnits_From_Options(options)
def get_gamma(R):
b=0.5
e = 1.
Xs = 10./(R**2+e**2) * options.fSx
Xb = ones(len(R)) * options.fSb
X = Xs + Xb
gamma = X**b/(1+X**b)
return gamma
#######################################
# 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
#######################################
# conversion 2.23e9 Msol/kpc^2 -> Msol/pc^2
fac_sdens = 2.23e9/(1000**2) * unit_mass
pt.rc('figure', figsize=(8,10))
# axes properties
#labelfont = 16
#ax = pt.gca()
#pt.setp(ax.get_xticklabels(),fontsize=labelfont)
#pt.setp(ax.get_yticklabels(),fontsize=labelfont)
pt.subplot(212)
#######################################
# LOOP
#######################################
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
################
nb_wg = nb.select('wg')
nb_cg = nb.select('cg')
nb_ba = nb.select('gas','stars','disk','bulge')
nb_vb = nb.select('wg','stars','disk','bulge')
nb_bu = nb.select('bulge')
nb_di = nb.select('disk')
r = getr(nr=31,nt=32,rm=rmax)
r = None
r,sdens = nb_ba.msdens(r=r,nb=nbins,rm=rmax)
# baryons
x0 = r.astype(float)
y0 = sdens.astype(float)*fac_sdens
r = getr(nr=31,nt=32,rm=rmax)
r = None
r,sdens = nb_wg.msdens(r=r,nb=nbins,rm=rmax)
# visible gas
x1 = r.astype(float)
y1 = sdens.astype(float)*fac_sdens
r = getr(nr=31,nt=32,rm=rmax)
r = None
r,sdens = nb_cg.msdens(r=r,nb=nbins,rm=rmax)
# dark gas
x2 = r.astype(float)
y2 = sdens.astype(float)*fac_sdens
r = getr(nr=31,nt=32,rm=rmax)
r = None
r,sdens = nb_vb.msdens(r=r,nb=nbins,rm=rmax)
# visible baryons
x3 = r.astype(float)
y3 = sdens.astype(float)*fac_sdens
r = getr(nr=31,nt=32,rm=rmax)
r = None
r,sdens = nb_bu.msdens(r=r,nb=nbins,rm=rmax)
# bulge
x4 = r.astype(float)
y4 = sdens.astype(float)*fac_sdens
r = getr(nr=31,nt=32,rm=rmax)
r = None
r,sdens = nb_di.msdens(r=r,nb=nbins,rm=rmax)
# disk
x5 = r.astype(float)
y5 = sdens.astype(float)*fac_sdens
x = concatenate((x0,x1,x2,x3))
y = concatenate((y0,y1,y2,y3))
# use log
#if log != None:
# x0,y0 = Graph_UseLog(x0,y0,log)
# x1,y1 = Graph_UseLog(x1,y1,log)
# x2,y2 = Graph_UseLog(x2,y2,log)
# x ,y = Graph_UseLog(x,y,log)
# g.ticksize(-1,0,-1,-1)
if file == files[0]:
xmin,xmax,ymin,ymax = pt.SetLimits(xmin,xmax,ymin,ymax,x,y)
# compute dark gas
gamma = get_gamma(x2)
y2 = y1 * (1-gamma)/gamma
#pt.plot(x0,y0,'k') # baryons
#pt.plot(x2,y2,'b') # dark gas
#pt.plot(x2,y1+y2,'k--') # total gas
pt.plot(x3,y3,'k') # visible baryons
pt.plot(x5,y3+y2,'k--') # baryons
pt.plot(x1,y1,'r') # visible gas
pt.plot(x4,y4,'g') # bulge
pt.plot(x5,y5,'c') # disk
pt.plot(x2,y2,'b') # dg
i = i + 1
pt.legend((
'total visible baryons',
'total baryons',
'visible gas',
'bulge',
'disk',
'dark gas'
),'upper right',shadow=True)
# add MW gas obs (Levine,Blitz, Heiles 2006)
#x = array([12,15,16,18,20])
#y = array([7,3,2.2,1.5,0.3])
#pt.plot(x,y,'o')
# ngc 5033
'''
x = array([0,8,16,24,32])
y = 10**array([1,1,0.85,0.5,0.25])
pt.plot(x,y,'o')
'''
#pt.xlabel(r'$\textrm{Radius}\,\left[ \rm{kpc} \right]$',fontsize=labelfont)
#pt.ylabel(r'$\textrm{Surface Density}\,\left[\rm{M}_\odot/\rm{pc}^2\right]$',fontsize=labelfont)
pt.xlabel('Radius [kpc]')
pt.ylabel('Surface Density [Ms/pc^3]')
pt.semilogy()
#pt.axis([xmin,xmax,ymin,ymax])
pt.axis([0,50,0.1,1e4])
#################################################
g = pt.subplot(211)
R = arange(0,51,1)
Rmin = min(R)
Rmax = max(R)
gamma = get_gamma(R)
pt.plot(R,gamma,'k',linewidth=2)
pt.xlabel('Radius [kpc]')
pt.ylabel('Visible / Dark gas')
#xmin,xmax,ymin,ymax = pt.SetLimits(None,None,0,1,R,gamma)
pt.axis([0,50,0,1])
#pt.grid()
# make the shaded regions
ix = R
iy = gamma
verts = [(Rmin,0)] + zip(ix,iy) + [(Rmax,0)]
poly = pt.Polygon(verts, facecolor='r', edgecolor='r')
g.add_patch(poly)
verts = [(Rmin,1)] + zip(ix,iy) + [(Rmax,1)]
poly = pt.Polygon(verts, facecolor='b', edgecolor='b')
g.add_patch(poly)
# -- end ---
if ps==None:
pt.show()
else:
pt.savefig(ps)
c4science · Help