Homec4science
Diffusion Gtools (master) Browse

gtools/scripts/gplotsph_r-val
8cfdf441e0b6master

Last Change
Raw File

gplotsph_r-val
View Options

#!/usr/bin/env python
'''
Extract and plot energy and mass values contained in the
output Gadget file called by default "energy.txt".
Yves Revaz
Fri Nov 28 14:32:09 CET 2008
'''
from numpy import *
from pNbody import *
import string
import sys
import os
import copy as docopy
from pNbody.libutil import histogram
from pNbody import libgrid
from optparse import OptionParser
from Gtools import *
from Gtools import io
import Ptools as pt
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_limits_options(parser)
parser = pt.add_log_options(parser)
parser = pt.add_ftype_options(parser)
parser = pt.add_reduc_options(parser)
parser = pt.add_center_options(parser)
parser = pt.add_display_options(parser)
parser = pt.add_select_options(parser)
parser = pt.add_cmd_options(parser)
parser = pt.add_info_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("--eps",
action="store",
dest="eps",
type="float",
default = 0.28,
help="smoothing length",
metavar=" FLOAT")
parser.add_option("--nmin",
action="store",
dest="nmin",
type="float",
default = 32,
help="min number of particles in a cell to accept value",
metavar=" INT")
parser.add_option("--G",
action="store",
dest="G",
type="float",
default = 1.,
help="Gravitational constant value",
metavar=" FLOAT")
parser.add_option("--ErrTolTheta",
action="store",
dest="ErrTolTheta",
type="float",
default = 0.5,
help="Error tolerance theta",
metavar=" FLOAT")
parser.add_option("--AdaptativeSoftenning",
action="store_true",
dest="AdaptativeSoftenning",
default = False,
help="AdaptativeSoftenning")
parser.add_option("--statfile",
action="store",
type="string",
dest="statfile",
default = 'stat.h5',
help="stat output file")
parser.add_option("--mode",
action="store",
type="string",
dest="mode",
default = 'den',
help="mode")
(options, args) = parser.parse_args()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
#######################################
# MakePlot
#######################################
def MakePlot(dirs,opt):
# some inits
palette = pt.GetPalette()
colors = pt.SetColorsForFiles(files,palette)
labels = []
#######################################
# LOOP
#######################################
# read files
for file in files:
######################################
# open file and apply option
######################################
nb = Nbody(file,ftype=opt.ftype)
nb = pt.do_select_options(nb,opt)
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
nb = pt.do_center_options(nb,opt)
nb = pt.do_info_options(nb,opt)
nb = pt.do_display_options(nb,opt)
######################################
# computes values
######################################
# create cylindrical rt grid
rc = 1.
g = lambda r:log(r/rc+1.)
gm = lambda r:rc*(exp(r)-1.)
g = None;gm=None
Gr = libgrid.Spherical_1d_Grid(rmin=0,rmax=opt.Rmax,nr=opt.nR,g=g,gm=gm)
# build the tree
opt.ComputePotential=False
if opt.ComputePotential:
print "ComputeTree"
nbt.getTree(force_computation=True,ErrTolTheta=opt.ErrTolTheta)
nb.getTree(force_computation=True,ErrTolTheta=opt.ErrTolTheta)
# radius
R = Gr.get_r()
stats = {}
stats['nR'] = opt.nR
stats['Rmax'] = opt.Rmax
stats['eps'] = opt.eps
stats['ErrTolTheta'] = opt.ErrTolTheta
stats['AdaptativeSoftenning'] = opt.AdaptativeSoftenning
stats['nmin'] = opt.nmin
stats['G'] = opt.G
stats['R'] = R
if opt.mode == 'den':
den = Gr.get_DensityMap(nb)
stats['Den'] = den
x = R
y = den
if opt.mode == 'sv' or opt.mode == 'sigma':
v = sqrt(nb.vx()**2+nb.vy()**2+nb.vz()**2)
sv = Gr.get_SigmaValMap(nb,v)
stats['sv'] = sv
x = R
y = sv
if opt.mode == 'sr':
sr = Gr.get_SigmaValMap(nb,nb.Vr())
stats['sr'] = sr
x = R
y = sr
if opt.mode == 'sz':
sz = Gr.get_SigmaValMap(nb,nb.Vz())
stats['sz'] = sz
x = R
y = sz
if opt.mode == 'st':
st = Gr.get_SigmaValMap(nb,nb.Vt())
stats['sz'] = st
x = R
y = st
###################################
# save hdf
##################s#################
pt.io.write_hdf5(opt.statfile,stats)
###################################
# now, plot
###################################
xmin,xmax,ymin,ymax = pt.SetLimits(opt.xmin,opt.xmax,opt.ymin,opt.ymax,x,y,opt.log)
pt.plot(x,y)
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.show()
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
#pt.figure(figsize=(8*2,6*2))
#pt.figure(dpi=10)
pt.pcolors
#fig = pt.gcf()
#fig.subplots_adjust(left=0.1)
#fig.subplots_adjust(right=1)
#fig.subplots_adjust(bottom=0.12)
#fig.subplots_adjust(top=0.95)
#fig.subplots_adjust(wspace=0.25)
#fig.subplots_adjust(hspace=0.02)
MakePlot(files,opt)
pt.EndPlot(files,opt)
c4science · Help