Page MenuHomec4science
Files F72773464

plotprofiles.py
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 16, 23:46

plotprofiles.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
'''
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
from scipy.optimize import leastsq
from scipy.optimize import bisection
#from mpfit import *
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.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("--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("--fit",
action="store_true",
dest="fit",
default = False,
help="fit profile")
parser.add_option("--val",
action="store",
type="string",
dest="val",
default = 'density',
help="value to plot")
(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':
x = G.get_r()
y = G.get_DensityMap(nb)
xlabel = r'$\rm{Radius}\,\left[ kpc \right]$'
ylabel = r'$\rm{Density}$'
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 = []
velocity_factor = 207.
# read files
for file in files:
nb = Nbody(file,ftype=opt.ftype)
# 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)
#nb_stars = nb.select('stars')
#nb_halo = nb.select('halo')
'''
################
# set the grid
################
if fct!=None:
print "f =lambda r:(%s)"%fct
exec("f =lambda r:(%s)"%fct)
else:
f = None
if fctm!=None:
print "fm =lambda r:(%s)"%fctm
exec("fm =lambda r:(%s)"%fctm)
else:
fm = None
'''
rc = 1
f = lambda r:log(r/rc+1.)
fm = lambda r:rc*(exp(r)-1.)
################
# get values
################
##########################################
# 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='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=='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)
#pt.scatter(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)

Event Timeline

c4science · Help