Homec4science
Diffusion Gtools (master) Browse

gtools/Gtools/gscripts/gExtractCylindrical
ecf771675681master

Last Change
Raw File

gExtractCylindrical
View Options

#!/usr/bin/python
'''
Extract and plot energy and mass values contained in the
output Gadget file called by default "energy.txt".
Yves Revaz
ven jun 9 10:43:59 CEST 2006
'''
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 = pt.add_units_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 = 64,
help="number of bins in R",
metavar=" INT")
parser.add_option("--nt",
action="store",
dest="nt",
type="int",
default = 64,
help="number of bins in t",
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("--xmode",
action="store",
dest="xmode",
type="float",
default = 2.,
help="mode for X-Toomre parameter",
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("--ComputeLambdaJeans",
action="store_true",
dest="ComputeLambdaJeans",
default = False,
help="Compute LambdaJ eans")
parser.add_option("--statfile",
action="store",
type="string",
dest="statfile",
default = 'stat.dmp',
help="stat output file")
parser.add_option("--disk",
action="store",
type="string",
dest="disk",
default = "('gas','disk','stars')",
help="stat output file")
parser.add_option("--components",
action="store",
dest="components",
type="string",
default = "('gas','halo','disk','bulge','stars')",
help="list of components",
metavar=" TUPLE")
parser.add_option("--x",
action="store",
type="string",
dest="x",
default = 'R',
help="x")
parser.add_option("--y",
action="store",
type="string",
dest="y",
default = 'vct',
help="y")
parser.add_option("--mode",
action="store",
type="string",
dest="mode",
default = 'all',
help="mode")
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()
if len(args) == 0:
print "you must specify a filename"
sys.exit(0)
files = args
return files,options
def get1dMeanFrom2dMap(mat_val,mat_num,nmin=32,axis=0):
m1 = sum(mat_num,axis)
m0 = sum(ones(mat_val.shape),axis)
vec_num = where((m0!=0),m1/m0,axis)
c = (mat_num>nmin)
m1 = sum(mat_val*c,axis)
m0 = sum(ones(mat_val.shape)*c,axis)
vec_sigma = where((m0!=0),m1/m0,0)
return vec_sigma
#######################################
# MakePlot
#######################################
def MakePlot(dirs,opt):
#######################################
# deal with mode
#######################################
if opt.mode=='all':
opt.ComputePotential = True
opt.ComputeSurfaceDensity = True
opt.ComputeDispertions = True
opt.ComputeStability = True
elif opt.mode=='Sden':
opt.ComputePotential = False
opt.ComputeSurfaceDensity = True
opt.ComputeDispertions = False
opt.ComputeStability = False
elif opt.mode=='velocities':
opt.ComputePotential = False
opt.ComputeSurfaceDensity = False
opt.ComputeDispertions = True
opt.ComputeStability = False
elif opt.mode=='stability':
opt.ComputePotential = True
opt.ComputeSurfaceDensity = True
opt.ComputeDispertions = True
opt.ComputeStability = True
#######################################
# LOOP
#######################################
# read files
for file in files:
######################################
# open file and apply option
######################################
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()
################
# apply options
################
nb = pt.do_reduc_options(nb,opt)
nb = pt.do_cmd_options(nb,opt)
###############
# create total
###############
nbt = docopy.deepcopy(nb)
###############
# create disk
###############
nbd = docopy.deepcopy(nb)
if opt.disk!=None:
opt.disk = eval(opt.disk, dict(__builtins__=None))
print "disk = ",opt.disk
nbd = nbd.select(opt.disk)
######################################################
# compute the total density in the plane
if opt.ComputeLambdaJeans:
nbtt = docopy.deepcopy(nb)
nbtt.set_tpe(0)
# use a grid
# create cylindrical rt grid
rc = 1.
g = lambda r:log(r/rc+1.)
gm = lambda r:rc*(exp(r)-1.)
g = None;gm=None
Grt = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.Rmax,nr=opt.nR,nt=opt.nt,z=0,g=g,gm=gm)
# radius
R,t = Grt.get_rt()
x = R
y = zeros(len(x))
z = zeros(len(x))
pos = transpose(array([x,y,z]))
pos = pos.astype(float32)
pos = Grt.get_Points()
nbg = Nbody(pos=pos)
rho, rsp = nbtt.ComputeDensityAndHsml(pos = pos)
Rho0 = Grt.get_MeanValMap(nbg,rho)
Rho0 = sum(Rho0,1)/opt.nt
######################################################
# here, we should loop over component
if opt.components != None:
opt.components = eval(opt.components, dict(__builtins__=None))
print "components = ",opt.components
for component in opt.components:
print "-----------------------------"
print "-- component = ",component
print "-----------------------------"
nb = nbt.select(component)
opt.statfile = "%s.dmp"%(component)
if nb.nbody>0:
###############
# other info
###############
#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
Grt = libgrid.Cylindrical_2drt_Grid(rmin=0,rmax=opt.Rmax,nr=opt.nR,nt=opt.nt,z=0,g=g,gm=gm)
# build the tree
if opt.ComputePotential:
print "ComputeTree"
nbt.getTree(force_computation=True,ErrTolTheta=opt.ErrTolTheta)
nb.getTree(force_computation=True,ErrTolTheta=opt.ErrTolTheta)
# radius
R,t = Grt.get_rt()
stats = {}
stats['nR'] = opt.nR
stats['nt'] = opt.nt
stats['Rmax'] = opt.Rmax
stats['eps'] = opt.eps
stats['ErrTolTheta'] = opt.ErrTolTheta
stats['AdaptativeSoftenning'] = opt.AdaptativeSoftenning
stats['xmode'] = opt.xmode
stats['nmin'] = opt.nmin
stats['G'] = opt.G
stats['R'] = R
stats['t'] = t
###################################
# Surface density
###################################
if opt.ComputeSurfaceDensity:
print "ComputeSurfaceDensity"
Sden = Grt.get_SurfaceDensityMap(nb)
Sden = sum(Sden,1)/opt.nt
Sdend = Grt.get_SurfaceDensityMap(nbd)
Sdend = sum(Sdend,1)/opt.nt
stats['Sden'] = Sden
stats['Sdend'] = Sdend
###################################
# rotation curve and frequencies
###################################
if opt.ComputePotential:
print "ComputePotential"
# radial acceleration (total)
Accx,Accy,Accz = Grt.get_AccelerationMap(nbt,eps=opt.eps,UseTree=True,AdaptativeSoftenning=opt.AdaptativeSoftenning)
# assuming cylindrical model
Ar = sqrt(Accx**2+Accy**2)
# general case
#x,y,z=Grt.get_xyz(offr=0.5)
#r = sqrt(x**2+y**2+z**2)
#Ar = -(Accx*x + Accy*y + Accz*z)/r
Ar = sum(Ar,1)/opt.nt
dPhit = Ar
d2Phit = libgrid.get_First_Derivative(dPhit,R)
kappat = libdisk.Kappa(R,dPhit,d2Phit)
omegat = libdisk.Omega(R,dPhit)
vct = libdisk.Vcirc(R,dPhit)
#nut = libdisk.Nu(z,Phit)
# radial acceleration (selection)
Accx,Accy,Accz = Grt.get_AccelerationMap(nb,eps=opt.eps,UseTree=True,AdaptativeSoftenning=opt.AdaptativeSoftenning)
Ar = sqrt(Accx**2+Accy**2)
Ar = sum(Ar,1)/opt.nt
dPhi = Ar
d2Phi = libgrid.get_First_Derivative(dPhi,R)
kappa = libdisk.Kappa(R,dPhi,d2Phi)
omega = libdisk.Omega(R,dPhi)
vc = libdisk.Vcirc(R,dPhi)
#nu = libdisk.Nu(z,Phi)
#stats['Phi'] = Phi
stats['dPhi'] = dPhi
stats['d2Phi'] = d2Phi
stats['kappa'] = kappa
stats['omega'] = omega
#stats['nu'] = nu
#stats['Phit'] = Phit
stats['dPhit'] = dPhit
stats['d2Phit'] = d2Phit
stats['kappat'] = kappat
stats['omegat'] = omegat
#stats['nut'] = nut
stats['vct'] = vct
stats['vc'] = vc
if opt.ComputeDispertions:
print "ComputeDispertions"
###################################
# number of points per cell
###################################
nn = Grt.get_NumberMap(nb)
nm = sum(nn,1)/opt.nt
###################################
# mean velocities
###################################
vr = Grt.get_MeanValMap(nb,nb.Vr())
vr = get1dMeanFrom2dMap(vr,nn,nmin=opt.nmin,axis=1)
vm = Grt.get_MeanValMap(nb,nb.Vt())
vm = get1dMeanFrom2dMap(vm,nn,nmin=opt.nmin,axis=1)
vz = Grt.get_MeanValMap(nb,nb.Vz())
vz = get1dMeanFrom2dMap(vz,nn,nmin=opt.nmin,axis=1)
if nb.u!=None:
u = Grt.get_MeanValMap(nb,nb.u)
u = get1dMeanFrom2dMap(u,nn,nmin=opt.nmin,axis=1)
else:
u = None
###################################
# velocity dispertions
###################################
sr = Grt.get_SigmaValMap(nb,nb.Vr())
sr = get1dMeanFrom2dMap(sr,nn,nmin=opt.nmin,axis=1)
sp = Grt.get_SigmaValMap(nb,nb.Vt())
sp = get1dMeanFrom2dMap(sp,nn,nmin=opt.nmin,axis=1)
sz = Grt.get_SigmaValMap(nb,nb.Vz())
sz = get1dMeanFrom2dMap(sz,nn,nmin=opt.nmin,axis=1)
stats['nn'] = nn
stats['nm'] = nm
stats['vr'] = vr
stats['vm'] = vm
stats['vz'] = vz
stats['sr'] = sr
stats['sp'] = sp
stats['sz'] = sz
stats['u'] = u
###################################
# Q,X,A
###################################
if opt.ComputeStability:
print "ComputeStability"
Q = libdisk.QToomre(opt.G,R,sr,kappat,Sdend)
X = libdisk.XToomre(opt.G,R,kappat,Sdend,opt.xmode)
A = libdisk.AAraki(sr,sz)
stats['Q'] = Q
stats['A'] = A
stats['X'] = X
# jeans length
if opt.ComputeLambdaJeans:
# sigma_R for Q=1
Q = 1.
srQ1 = 3.36*opt.G*Sdend/kappat * Q
stats['srQ1'] = srQ1
lambdaj = pi*srQ1**2/opt.G/Rho0
stats['lambdaj'] = lambdaj
###################################
# add units
###################################
stats['localsystem_of_units'] = nb.localsystem_of_units
###################################
# save output
###################################
pt.io.write_dmp(opt.statfile,stats)
if __name__ == '__main__':
files,opt = parse_options()
MakePlot(files,opt)
c4science · Help