Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F122165889
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
Wed, Jul 16, 07:12
Size
18 KB
Mime Type
text/x-python
Expires
Fri, Jul 18, 07:12 (2 d)
Engine
blob
Format
Raw Data
Handle
27443275
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")
parser.add_option("--PJeans",
action="store_true",
dest="PJeans",
default = False,
help="add Jeans Pressure Floor")
parser.add_option("--JeansHsml",
action="store",
type="float",
dest="JeansHsml",
default = None,
help="resolution for the Jeans pressure floor")
parser.add_option("--NJ",
action="store",
type="float",
dest="NJ",
default = 10.,
help="Jeans factor for jeans pressure floor")
parser.add_option("--Softening",
action="store",
type="float",
dest="Softening",
default = None,
help="Gravitational softening for the pressure floor")
parser.add_option("--SofteningMaxPhys",
action="store",
type="float",
dest="SofteningMaxPhys",
default = None,
help="Gravitational softening (maxphys) for the pressure floor")
(options, args) = parser.parse_args()
pt.check_files_number(args)
files = args
return files,options
def get_value_and_label(nb,val,args=None):
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 == 'logTJeans':
label = r'$\rm{log\,Jeans\,Temperature}\,\left[ \rm{K} \right]$'
Hsml = None
if args!=None:
if args.has_key("JeansHsml"):
Hsml = args["JeansHsml"]
Softening = None
if args!=None:
if args.has_key("Softening"):
Softening = args["Softening"]
SofteningMaxPhys = None
if args!=None:
if args.has_key("SofteningMaxPhys"):
SofteningMaxPhys = args["SofteningMaxPhys"]
val = nb.TJeans(Hsml=Hsml,Softening=Softening,SofteningMaxPhys=SofteningMaxPhys)
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)
if opt.PJeans:
'''
args={"JeansHsml":opt.JeansHsml,"Softening":opt.Softening,"SofteningMaxPhys":opt.SofteningMaxPhys}
yj,Plabel = get_value_and_label(nb,"logTJeans",args=args)
data = pt.DataPoints(x,yj,color=colors.get(),label=file,tpe='points')
datas.append(data)
'''
# some constants
gamma = nb.unitsparameters.get('gamma')
xi = nb.unitsparameters.get('xi')
ionisation = nb.unitsparameters.get('ionisation')
mu = thermodyn.MeanWeight(xi,ionisation)
mh = ctes.PROTONMASS.into(nb.localsystem_of_units)
k = ctes.BOLTZMANN.into(nb.localsystem_of_units)
G = ctes.GRAVITY.into(nb.localsystem_of_units)
NJ = opt.NJ # Jeans Mass factor
T = nb.T()
if opt.JeansHsml==None and opt.Softening==None and opt.SofteningMaxPhys==None:
Hsml = nb.SphRadius()
else:
if opt.Softening!=None and opt.SofteningMaxPhys!=None:
print
print " using Softening = %g and SofteningMaxPhys = %g"%(opt.Softening,opt.SofteningMaxPhys)
Hsml = nb.ComputeSofteningCosmo(opt.Softening,opt.SofteningMaxPhys)
print Hsml
print
else:
Hsml = opt.JeansHsml
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])
rho_plot = nb.Rho(units=out_units.UnitDensity)
rho = nb.Rho()
#####################################
# compute the crittical temperature
#####################################
TJ = (mu*mh)/k * 4./pi * G/gamma * NJ**(2./3.) * Hsml**(2) * rho
x = log10(rho_plot)
y = log10(TJ)
data = pt.DataPoints(x,y,color='y',label=file,tpe='points')
datas.append(data)
#####################################
# compute the crittical density
#####################################
rhocJ = pi/4.0 * k/(mu*mh) * gamma/G * NJ**(-2./3.) * Hsml**(-2) * T
rhocJ = rhocJ * nb.localsystem_of_units.convertionFactorTo(out_units.UnitDensity)
x = log10(rhocJ)
y = log10(T)
#data = pt.DataPoints(x,y,color='r',label=file,tpe='points')
#datas.append(data)
#####################################
# unstable regions
#####################################
c = rho_plot > rhocJ
x = log10(compress(c,rho_plot))
y = log10(compress(c,T))
data = pt.DataPoints(x,y,color='b',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