Homec4science
Diffusion Gtools (master) Browse

gtools/dev/gplot_r-sdens_cw
8cfdf441e0b6master

Last Change
Raw File

gplot_r-sdens_cw
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 Ptools as pt
import string
import sys
import os
from pNbody.libutil import histogram
from optparse import OptionParser
from Gtools import *
from Gtools import vanderwaals as vw
unit_mass = 0.774
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("--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
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])
#######################################
# 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)
fig = pt.gcf()
fig.subplots_adjust(left=0.11)
fig.subplots_adjust(top=1)
fig.subplots_adjust(bottom=0.1)
fig.subplots_adjust(right=1)
fig.subplots_adjust(wspace=0.15)
fig.subplots_adjust(hspace=0.1)
fac_sdens = 2.23e9/(1000**2) * unit_mass
#######################################
# LOOP
#######################################
i = 0
lss = ['-','--',':']
labels = ['N=1','N=3','N=5']
# 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')
r = getr(nr=31,nt=32,rm=rmax)
r = None
r,sdens = nb_ba.msdens(r=r,nb=nbins,rm=rmax)
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)
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)
x2 = r.astype(float)
y2 = sdens.astype(float)*fac_sdens
x = concatenate((x0,x1,y1))
y = concatenate((y0,x2,y2))
if file == files[0]:
xmin,xmax,ymin,ymax = pt.SetLimits(options.xmin,options.xmax,options.ymin,options.ymax,x,y,options.log)
pt.plot(x0,y0,c='k',ls=lss[i],label=r"$\rm{%s : total\,\,baryons}$"%labels[i])
pt.plot(x1,y1,c='r',ls=lss[i],label=r"$\rm{%s : visible\,\,gas}$"%labels[i])
pt.plot(x2,y2,c='b',ls=lss[i],label=r"$\rm{%s : dark\,\, gas}$"%labels[i])
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{Surface\,\,Density}\,\left[ M_\odot/pc^2 \right]$',fontsize=labelfont)
# add legend
if options.add_legend:
#labels = [r'$N=1$ : visible gas',r'$N=3$ : visible gas',r'$N=3$ : dark gas']
#pt.legend(crv_lst,'upper right',shadow=True)
pt.legend(loc='upper right',shadow=True,labelsep=0.0001)
####################################
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