gExtractCylindrical
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Feb 1, 10:46

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(Accx2+Accy2)

	# general case
	#x,y,z=Grt.get_xyz(offr=0.5)
	#r = sqrt(x2+y2+z**2)
	#Ar = -(Accxx + Accyy + 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(Accx2+Accy2)
	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.36opt.GSdend/kappat * Q
	stats['srQ1'] = srQ1

	lambdaj = pisrQ1*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)

gExtractCylindricalNo OneTemporaryActions

File Metadata

gExtractCylindricalView Options

Event Timeline

gExtractCylindrical
No OneTemporary
Actions

gExtractCylindrical
View Options