Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65307980
pplot
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
Sun, Jun 2, 19:41
Size
14 KB
Mime Type
text/x-python
Expires
Tue, Jun 4, 19:41 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18046800
Attached To
rGTOOLS Gtools
pplot
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)
parser = 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("--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")
parser.add_option("--forceComovingIntegrationOff",
action="store_true",
dest="forceComovingIntegrationOff",
default = False,
help="force the model not 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]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Rxyz(units=out_units.UnitLength)
return val,label
if val == 'logR':
label = r'$\rm{log Radius}\,\left[ kpc \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Rxyz(units=out_units.UnitLength)
val = log10(val)
return val,label
if val == 'Rxy':
label = r'$\rm{Radius}\,\left[ kpc \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Rxy(units=out_units.UnitLength)
return val,label
if val == 'Hsml':
label = r'$\rm{SPH\,Smoothin\,Length}\,\left[ kpc \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.SphRadius(units=out_units.UnitLength)
return val,label
if val == 'logHsml':
label = r'$\rm{log SPH\,Smoothin\,Length}\,\left[ kpc \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = log10(nb.SphRadius(units=out_units.UnitLength))
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 == 'logTJean':
label = r'$\rm{log\,Jeans\,Temperature}\,\left[ \rm{K} \right]$'
val = nb.TJeans()
val = log10(val)
return val,label
elif val == 'logResolvedSPHMass':
label = r'$\rm{log\,Mass}\,\left[ \rm{Msol} \right]$'
out_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_s,units.Unit_K])
h = nb.SphRadius()
rho = nb.Rho()
val = 4./3.*pi*rho*h**3 * nb.localsystem_of_units.convertionFactorTo(out_units.UnitMass)
val = log10(val)
return val,label
elif val == 'rho':
label = r'$\rm{Density}\,\left[ \rm{atom/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])
val = nb.Rho(units=out_units.UnitDensity)
return val,label
elif val == 'logrho':
label = r'$\rm{log\,Density}\,\left[ \rm{atom/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])
val = nb.Rho(units=out_units.UnitDensity)
val = log10(val)
return val,label
elif val == 'Tcool':
label = r'$\rm{Cooling\,Time}\,\left[ \rm{Myr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Tcool(units=out_units.UnitTime)
return val,label
elif val == 'logTcool':
label = r'$\rm{log\,Cooling\,Time}\,\left[ \rm{Myr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Tcool(units=out_units.UnitTime)
val = log10(val)
return val,label
elif val == 'Tff':
label = r'$\rm{Free\,Fall\,Time}\,\left[ \rm{Myr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Tff(units=out_units.UnitTime)
return val,label
elif val == 'logTff':
label = r'$\rm{log\,Free\,Fall\,Time}\,\left[ \rm{Myr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Myr,units.Unit_K])
val = nb.Tff(units=out_units.UnitTime)
val = log10(val)
return val,label
elif val == 'Age':
label = r'$\rm{Age}\,\left[ \rm{Gyr} \right]$'
out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Ms,units.Unit_Gyr,units.Unit_K])
val = nb.StellarAge(units=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
# abundance ratio : X/Y
elif val.find("/")!=-1:
elt1,elt2=val.split('/')
label = r'$\left[ \rm{%s}/\rm{%s} \right]$'%(elt1,elt2)
val = nb.AbRatio(elt1,elt2)
return val,label
else:
label = r'%s'%val
print "val = %s"%val
exec("val = %s"%val)
return val, label
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
datas = []
# read files
for file in files:
nb = Nbody(file,ftype=opt.ftype)
################
# units
################
# define local units
unit_params = pt.do_units_options(opt)
nb.set_local_system_of_units(params=unit_params)
# define output units
# nb.ToPhysicalUnits()
if opt.forceComovingIntegrationOn:
nb.setComovingIntegrationOn()
if opt.forceComovingIntegrationOff:
nb.setComovingIntegrationOff()
################
# 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)
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)
pt.scatter(d.x,d.y,c=d.color,edgecolor=d.color,s=1,linewidths=1,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))
if sum(ravel(data))>0:
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
if sum(ravel(data))>0:
# 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)
Event Timeline
Log In to Comment