Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F122393819
gplot_r-sdens_00
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
Thu, Jul 17, 15:30
Size
6 KB
Mime Type
text/x-python
Expires
Sat, Jul 19, 15:30 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
27455703
Attached To
rGTOOLS Gtools
gplot_r-sdens_00
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")
(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)
#######################################
# 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
mode
# axes properties
labelfont = 16
ax = pt.gca()
pt.setp(ax.get_xticklabels(),fontsize=labelfont)
pt.setp(ax.get_yticklabels(),fontsize=labelfont)
#######################################
# 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,log)
#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(x1,y1,'r') # visible gas
pt.plot(x4,y4,'g') # bulge
pt.plot(x5,y5,'c') # disk
i = i + 1
'''
pt.legend((
r'$\textrm{total baryons}$',
r'$\textrm{visible gas}$',
r'$\textrm{dark gas}$',
r'$\textrm{visible baryons}$'
),'upper right',shadow=True)
'''
pt.legend((
'total visible baryons',
'visible gas',
'bulge',
'disk'
),'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])
# -- end ---
if ps==None:
pt.show()
else:
pt.savefig(ps)
Event Timeline
Log In to Comment