Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/plots/plot_scalars/pplot.py
0bbffe913903master

Last Change
Raw File

pplot.py
View Options

#!/usr/bin/env python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
from pNbody import ctes
from pNbody import thermodyn
import string
from scipy import optimize
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)
parse = pt.add_units_options(parser)
parser.add_option("--x",
action="store",
dest="x",
type="string",
default = 'r',
help="x value to plot",
metavar=" STRING")
parser.add_option("--y",
action="store",
dest="y",
type="string",
default = 'T',
help="y value to plot",
metavar=" STRING")
parser.add_option("--z",
action="store",
dest="z",
type="string",
default = None,
help="z value to plot",
metavar=" STRING")
parser.add_option("--legend",
action="store_true",
dest="legend",
default = False,
help="add a legend")
parser.add_option("--colorbar",
action="store_true",
dest="colorbar",
default = False,
help="add a colorbar")
parser.add_option("--nx",
action="store",
type="int",
dest="nx",
default = 64,
help="number of bins in x")
parser.add_option("--ny",
action="store",
type="int",
dest="ny",
default = 64,
help="number of bins in y")
parser.add_option("--nxh",
action="store",
type="int",
dest="nxh",
default = 0,
help="number of bins in x")
parser.add_option("--nopoints",
action="store_true",
dest="nopoints",
default = False,
help="do not plot points")
parser.add_option("--map",
action="store_true",
dest="map",
default = False,
help="plot map instead of points")
parser.add_option("--accumulate",
action="store_true",
dest="accumulate",
default = False,
help="integrate histogramm")
parser.add_option("--density",
action="store_true",
dest="density",
default = False,
help="compute density")
parser.add_option("--data",
action="store",
type="string",
dest="data",
default = None,
help="data")
parser.add_option("--mc",
action="store",
type="int",
dest="mc",
default = 0,
help="number montecarlo point")
parser.add_option("--forceComovingIntegrationOn",
action="store_true",
dest="forceComovingIntegrationOn",
default = False,
help="force the model to be in in comoving integration")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
def get_value_and_label(nb,val):
if val == 'R':
label = r'$\rm{Radius}\,\left[ kpc \right]$'
val = nb.Rxyz()
return val,label
elif val == 'T':
label = r'$\rm{Temperature}\,\left[ \rm{K} \right]$'
val = nb.T()
return val,label
elif val == 'logT':
label = r'$\rm{log\,Temperature}\,\left[ \rm{K} \right]$'
val = nb.T()
val = log10(val)
return val,label
elif val == 'rho':
label = r'$\rm{Density}\,\left[ \rm{atom/cm^3} \right]$'
val = nb.Rho()
# unit conversion
fact_g = units.PhysCte(1,nb.localsystem_of_units.UnitDic['kg'])
fact_g = fact_g.into(units.cgs)
fact_cm = units.PhysCte(1,nb.localsystem_of_units.UnitDic['m'])
fact_cm = fact_cm.into(units.cgs)
# factor conversion CU->cgs
fact_density = fact_g/fact_cm**3
# g -> H atom
fact_density = fact_density/ctes.PROTONMASS.into(units.cgs)
val = val*fact_density
return val,label
elif val == 'logrho':
label = r'$\rm{log\,Density}\,\left[ \rm{atom/cm^3} \right]$'
val = nb.Rho()
# unit conversion
fact_g = units.PhysCte(1,nb.localsystem_of_units.UnitDic['kg'])
fact_g = fact_g.into(units.cgs)
fact_cm = units.PhysCte(1,nb.localsystem_of_units.UnitDic['m'])
fact_cm = fact_cm.into(units.cgs)
# factor conversion CU->cgs
fact_density = fact_g/fact_cm**3
# g -> H atom
fact_density = fact_density/ctes.PROTONMASS.into(units.cgs)
val = val*fact_density
val = log10(val)
return val,label
elif val == 'Age':
label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
val = nb.Age()
return val,label
elif val == 'age':
label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Msol,units.Unit_Gyr,units.Unit_K])
val = nb.age() * nb.localsystem_of_units.convertionFactorTo(out_units.UnitTime)
return val,label
elif val == 'aFe':
label = r'$\left[ \rm{Fe}/\rm{H} \right]$'
val = nb.metals[:,0]
return val,label
elif val == 'aMgFe':
label = r'$\left[ \rm{Fe}/\rm{H} \right]$'
val = nb.metals[:,1]/(nb.metals[:,0]+1e-20)
return val,label
elif val == 'Fe':
label = r'$\left[ \rm{Fe}/\rm{H} \right]$'
val = nb.Fe()
return val,label
elif val == 'Mg':
label = r'$\left[ \rm{Mg}/\rm{H} \right]$'
val = nb.Mg()
return val,label
elif val == 'MgFe':
label = r'$\left[ \rm{Mg}/\rm{Fe} \right]$'
val = nb.MgFe()
return val,label
else:
label = r'%s'%val
print "val = %s"%val
exec("val = %s"%val)
return val, label
#################################
def ComputeHisto(x,m,v,xmin,xmax,nx):
#################################
'''
Compute an histogram
'''
shape = (nx,)
# normalize x
x = (x-xmin)/(xmax-xmin)
x = x.astype(float32)
m = m.astype(float32)
v = v.astype(float32)
# compute the map
h = mkmap1d(x,m,v,shape).astype(float)
# compute x
x = linspace(xmin,xmax,nx)
return x,h
#################################
def ComputeMeanHisto(x,m,v,xmin,xmax,nx):
#################################
"""
Compute the mean histogram
"""
x,h0 = ComputeHisto(x,m,ones(len(m)),xmin,xmax,nx)
x,h1 = ComputeHisto(x,m,v,xmin,xmax,nx)
h1 = where(h0==0,0,h1)
h0 = where(h0==0,1,h0)
h = h1/h0
return x,h
def get_histo(nb,val1,val2,val1min,val1max,n):
if (val2=='m') or (val2=='Lum'):
mode= 'extensive'
else:
mode= 'intensive'
print val2
sys.exit()
val1,l1 = get_value_and_label(nb,val1)
val2,l2 = get_value_and_label(nb,val2)
if (mode=='extensive'):
m = ones(nb.nbody)
x,vx=ComputeHisto(val1,m,val2,val1min,val1max,n)
else:
m = nb.mass
x,vx=ComputeMeanHisto(val1,m,val2,val1min,val1max,n)
xlabel = l1
ylabel = l2
return x,vx,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)
# read 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 "---------------------------------------------------------"
################
# get values
################
if not opt.nopoints:
x,xlabel = get_value_and_label(nb,opt.x)
y,ylabel = get_value_and_label(nb,opt.y)
if opt.z!=None:
z,zlabel = get_value_and_label(nb,opt.z)
data = pt.DataPoints(x,y,color=z,label=file,tpe='points')
else:
data = pt.DataPoints(x,y,color=colors.get(),label=file,tpe='points')
datas.append(data)
################
# get mean values
################
if opt.nxh>0:
v,h,xlabel,ylabel = get_histo(nb,opt.x,opt.y,opt.xmin,opt.xmax,opt.nxh)
if not opt.accumulate:
x = v
y = h
else:
x = v
y = add.accumulate(h)
if opt.density:
V = 4./3.*pi*(v[1:])**3 - 4./3.*pi*(v[:-1])**3
y = y[:-1]/V
x = x[:-1]
data = pt.DataPoints(x,y,color=colors.get(),label=file,tpe='both')
datas.append(data)
norm = None
# now, plot
if not opt.map :
for d in datas:
if d.tpe=='points' or d.tpe=='both':
if opt.log!=None:
if string.find(opt.log,'z')!=-1:
#norm = pt.colors.LogNorm(opt.zmin,opt.zmax)
norm = pt.colors.LogNorm()
else:
norm = None
cmap = pt.GetColormap('rainbow4',revesed=True)
pt.scatter(d.x,d.y,c=d.color,s=5,linewidths=0,marker='o',vmin=opt.zmin,vmax=opt.zmax,norm=norm,cmap=cmap)
if d.tpe=='line' or d.tpe=='both':
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)
if opt.map :
# now, plot
for d in datas:
x = d.x
y = d.y
z = zeros(len(x))
# set log
#if log!=None:
# if string.find(opt.log,'x') != -1:
# x = log10(x)
# if string.find(opt.log,'y') != -1:
# y = log10(y)
c = isreal(x)*isreal(y)
x = compress(c,x)
y = compress(c,y)
# re-scale between 0 and 1
xs = 1 - (x-xmax)/(xmin-xmax)
ys = 1 - (y-ymax)/(ymin-ymax)
pos = transpose(array([xs,ys,z],float32))
m = ones(len(x),float32)
data = mapping.mkmap2d(pos,m,m,(opt.nx,opt.ny))
data = transpose(data)/sum(ravel(data))
#########################
# compute and draw level
#########################
rdata = ravel(data)
rdata = rdata/sum(rdata)
zmin = min(rdata)
zmax = max(rdata)
def levfct(x,flev):
da = where(rdata<x,0.,rdata)
return sum(da)-flev
# here, we cut at 90% of the flux
levelmax = zmax
levelmin = optimize.bisect(levfct, a=zmin, b=zmax, args = (0.90,), xtol = 1e-10, maxiter = 100)
cmap = pt.GetColormap('heat',revesed=True)
#im = pt.imshow(data, interpolation='bilinear', origin='lower',cmap=cmap, extent=(xmin,xmax,ymin,ymax),aspect='auto')
im = pt.imshow(data, interpolation='bilinear', origin='lower',cmap=cmap, extent=(xmin,xmax,ymin,ymax),aspect='auto',vmin=levelmin,vmax=levelmax)
#cset = pt.contour(data,array([levelmin]),origin='lower',linewidths=1,extent=(xmin,xmax,ymin,ymax))
#########################
# compute monte carlo
#########################
if opt.mc > 0:
from pNbody import montecarlolib
datamc = transpose(data)
pos = montecarlolib.mc2d(datamc.astype(float32),opt.mc,random.random()*1000)
# transform into physical coord
x = pos[:,0]
y = pos[:,1]
x = x* (xmax-xmin) + xmin
y = y* (ymax-ymin) + ymin
#pt.scatter(x,y,marker='o',s=1)
xe = 0.1*ones(len(y))
ye = 0.2*ones(len(y))
pt.errorbar(x,y,xerr=xe,yerr=ye,fmt=None,ecolor='g',lw=0.1)
# 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)
if opt.legend:
pt.LegendFromDataPoints(datas)
if opt.colorbar:
#lev = array([-2,0,15])
#l_f = pt.LogFormatter(10, labelOnlyBase=False)
#cb = pt.colorbar(ticks=lev, format = l_f)
cb = pt.colorbar()
cb.set_label(zlabel)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
c4science · Help