test_density.py
de93a8c7a339master

test_density.py
View Options

			from numarray import *
			import SM

			from pNbody import *

			import time

			import treelib


			eps = 0.1

			# open a model
			nb = Nbody('disk.dat',ftype='gadget')
			#nb = nb.selectc(nb.rxyz()<2)
			nb = nb.select('gas')

			# create a tree object
			print "make tree"
			MyTree = treelib.Tree(npart=nb.npart,pos=nb.pos,mass=nb.mass)



			'''
			pos = array([[0,0,0]])
			pos = pos.astype(Float32)

			hsml = array([[10]])
			hsml = hsml.astype(Float32)


			density = MyTree.Density(pos,hsml)

			print density
			'''




			R = arange(0,20,0.1)
			#pos = concatenate((R,zeros(len(R)),zeros(len(R))))
			#pos.shape = (3,len(R))
			#pos = transpose(pos)
			#pos = pos.astype(Float32)
			pos = nb.pos
			R = nb.rxyz()
			#hsml = 1*ones(len(pos)).astype(Float32)
			print "init hsml"
			hsml = MyTree.InitHsml(33,3)
			print "compute density"
			density,hsml = MyTree.Density(pos,hsml,33,3)

			g = SM.plot()
			#g.limits(R,pots2)
			g.limits(R,density)
			g.box()
			g.points(R,density)
			g.show()