Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/plots/plot_profiles/plot_spherical_profile.py
0bbffe913903master

Last Change
Raw File

plot_spherical_profile.py
View Options

#!/usr/bin/env python
'''
Plot value computed in a spherical grid, as a function of the radius
Yves Revaz
Wed Jan 9 09:29:46 CET 2008
revison Fri Nov 25 11:01:42 CET 2011
'''
import Ptools as pt
from pNbody import *
from pNbody import myNumeric
from pNbody import libdisk
from pNbody import libgrid
from pNbody import profiles
from pNbody import cosmo
import string
import sys
import os
from numpy import *
from optparse import OptionParser
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_reduc_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_select_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_display_options(parser)
parser = pt.add_info_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser = pt.add_legend_options(parser)
parse = pt.add_units_options(parser)
parser.add_option("--forceComovingIntegrationOn",
action="store_true",
dest="forceComovingIntegrationOn",
default = False,
help="force the model to be in in comoving integration")
parser.add_option("--fct",
action="store",
dest="fct",
type="string",
default = None,
help="transformation function ex : 'r**3' should be of the form M(r)",
metavar=" STR")
parser.add_option("--fctm",
action="store",
dest="fctm",
type="string",
default = None,
help="inverse transformation function ex : 'r**(1/3.)'",
metavar=" STR")
parser.add_option("--val",
action="store",
type="string",
dest="val",
default = 'density',
help="value to plot")
parser.add_option("--rmax",
action="store",
dest="rmax",
type="float",
default = 50.,
help="max radius of bins",
metavar=" FLOAT")
parser.add_option("--nr",
action="store",
dest="nr",
type="int",
default = 32,
help="number of bins in r",
metavar=" INT")
parser.add_option("--HubbleParam",
action="store",
dest="HubbleParam",
type="float",
default = None,
help="force the Hubble parameter value during units conversion",
metavar=" FLOAT")
parser.add_option("--atime",
action="store",
dest="atime",
type="float",
default = None,
help="force the atime value during units conversion",
metavar=" FLOAT")
(options, args) = parser.parse_args()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#######################################
# compute values
#######################################
def get_val(nb,val,G,donotclean=False):
if val=='density':
# g/cm^3
#out_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_g,units.Unit_s,units.Unit_K])
#ylabel = r'$\rm{Density}\,\left[ \rm{g/cm^3} \right]$'
#Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
#Unit_atom.set_symbol('atom')
#out_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_s,units.Unit_K])
#ylabel = r'$\rm{Density}\,\left[ \rm{atom/cm^3} \right]$'
# Msol/kpc^3
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Msol,units.Unit_s,units.Unit_K])
ylabel = r'$\rm{Density}\,\left[ \rm{kpc/M_\odot} \right]$'
x = G.get_r()
y = G.get_DensityMap(nb) * nb.localsystem_of_units.convertionFactorTo(out_units.UnitDensity)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
elif val=='mass':
x = G.get_r()
y = G.get_MassMap(nb)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Mass}$'
elif val=='imass':
x = G.get_r()
y = G.get_MassMap(nb)
y = add.accumulate(y)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Cumulated\,Mass}$'
elif val=='sdens':
x,z = G.get_rz()
y = G.get_SurfaceDensityMap(nb)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Surface\,Density}$'
elif val=='massproj':
x,t = G.get_rt()
y = G.get_MassMap(nb)
y = sum(y,axis=1)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Mass}$'
elif val=='imassproj':
x,t = G.get_rt()
y = G.get_MassMap(nb)
y = sum(y,axis=1)
y = add.accumulate(y)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Mass}$'
elif val=='massproj2':
x = G.get_r()
y = G.get_MassMap(nb.rxy(),nb.mass)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Mass}$'
elif val=='imassproj2':
x = G.get_r()
y = G.get_MassMap(nb.rxy(),nb.mass)
y = add.accumulate(y)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Mass}$'
if not donotclean:
x,y=pt.CleanVectorsForLogX(x,y)
x,y=pt.CleanVectorsForLogY(x,y)
return x,y,xlabel,ylabel
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
datas = []
# define local units
unit_params = pt.do_units_options(opt)
#######################
# set the grid scaling
#######################
if opt.fct!=None and opt.fctm!=None:
print "this does not works yes..."
sys.exit()
#print "f =lambda r:(%s)"%opt.fct
#exec("f =lambda r:(%s)"%opt.fct)
#print "fm =lambda r:(%s)"%opt.fctm
#exec("fm =lambda r:(%s)"%opt.fctm)
else:
rc = 1
f = lambda r:log(r/rc+1.)
fm = lambda r:rc*(exp(r)-1.)
#################
# loop over files
#################
for file in files:
nb = Nbody(file,ftype=opt.ftype)
################
# units
################
# define local units
nb.set_local_system_of_units(params=unit_params)
# define output units
# nb.ToPhysicalUnits()
if opt.forceComovingIntegrationOn:
nb.setComovingIntegrationOn()
################
# apply options
################
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_select_options(nb,opt)
nb = pt.do_center_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
nb = pt.do_info_options(nb,opt)
nb = pt.do_display_options(nb,opt)
################
# some info
################
print "---------------------------------------------------------"
nb.localsystem_of_units.info()
nb.ComovingIntegrationInfo()
print "---------------------------------------------------------"
##############################
# convert to physical units
##############################
nb.toPhysicalUnits(a=opt.atime,h=opt.HubbleParam)
##########################################
# Spherical_1d_Grid
##########################################
if opt.val=='density':
G = libgrid.Spherical_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color=colors.next(),label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='mass':
G = libgrid.Spherical_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='imass':
G = libgrid.Spherical_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='bf':
G = libgrid.Spherical_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
xs,ys,xlabel,ylabel=get_val(nb_stars,'mass',G,donotclean=True)
xh,yh,xlabel,ylabel=get_val(nb_halo,'mass',G,donotclean=True)
x = xs
y = ys/(yh+ys)
x,y=pt.CleanVectorsForLogX(x,y)
x,y=pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
ylabel = r'$\rm{Baryonic\,Fraction}$'
elif opt.val=='ibf':
G = libgrid.Spherical_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
xs,ys,xlabel,ylabel=get_val(nb_stars,'imass',G,donotclean=True)
xh,yh,xlabel,ylabel=get_val(nb_halo,'imass',G,donotclean=True)
x = xs
y = ys/(yh+ys)
x,y=pt.CleanVectorsForLogX(x,y)
x,y=pt.CleanVectorsForLogY(x,y)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
ylabel = r'$\rm{Cumulated\,Baryonic\,Fraction}$'
elif opt.val=='sigma':
G = libgrid.Spherical_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
x = G.get_r()
y = G.get_SigmaValMap(nb,nb.Vz()) * velocity_factor
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Velocity\,Dispertion}$'
##########################################
# Cylindrical_2drt_Grid
##########################################
elif opt.val=='sdens':
# use a cylindrical grid
G = libgrid.Cylindrical_2drz_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,zmin=-10000,zmax=10000,nz=3,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='massproj':
# use a cylindrical grid
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='imassproj':
# use a cylindrical grid
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='sigmaproj':
G = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,nt=1,g=f,gm=fm)
x,t = G.get_rt()
y = G.get_SigmaValMap(nb,nb.Vz())
y = sum(y,axis=1)
y = y*velocity_factor
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Velocity\,Dispertion}$'
##########################################
# Generic_1d_Grid
##########################################
'''
we get the same results than massproj, imassproj and sigmaproj
'''
elif opt.val=='massproj2':
G = libgrid.Generic_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='imassproj2':
G = libgrid.Generic_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
x ,y ,xlabel,ylabel=get_val(nb,opt.val,G)
#xs,ys,xlabel,ylabel=get_val(nb_stars,opt.val,G)
#xh,yh,xlabel,ylabel=get_val(nb_halo,opt.val,G)
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
#datas.append(pt.DataPoints(xs,ys,color='r',label=None,tpe='points'))
#datas.append(pt.DataPoints(xh,yh,color='b',label=None,tpe='points'))
elif opt.val=='sigmaproj2':
G = libgrid.Generic_1d_Grid(rmin=0,rmax=opt.rmax,nr=opt.nr,g=f,gm=fm)
x = G.get_r()
y,m,mv,mv2 = G.get_SigmaMap(nb.rxy(),nb.mass,nb.Vz())
y = y*velocity_factor
datas.append(pt.DataPoints(x, y, color='k',label=None,tpe='points'))
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Velocity\,Dispertion}$'
###############
# plot
###############
for d in datas:
pt.plot(d.x,d.y,color=d.color)
# set limits and draw axis
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.grid(False)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
c4science · Help