test02-direct.py
229daad692b3master

test02-direct.py
View Options

			#!/usr/bin/env python


			from pNbody import *
			from pNbody import ic
			import time
			import treelib

			import Ptools as pt

			ErrTolTheta = 0.7
			eps = 0.1

			# open a model
			#nb = Nbody("./Schaller/galaxy.dat",ftype='gadget')

			random.seed(1)
			e = 5.
			nb = ic.plummer(524288,1,1,1,e,200,ftype='gadget')
			#nb = ic.homosphere(100000,1,1,1,ftype='gadget')



			nb.rename('snap.dat')
			nb.write()

			MaxPartPerNode=10

			# create a tree object
			t1 = time.time()
			print "start tree"
			MyTree = treelib.Tree(npart=array(nb.npart),pos=nb.pos,vel=nb.vel,mass=nb.mass,ErrTolTheta=ErrTolTheta,MaxPartPerNode=MaxPartPerNode)
			print "tree done."
			t2 = time.time()
			print "time = ",t2-t1

			#pot = MyTree.AllDirectPotential(eps)
			t1 = time.time()
			#pot = MyTree.AllPotentialNpart(eps)
			pot = MyTree.AllPotential(eps)
			t2 = time.time()
			print "time = ",t2-t1


			Direct=False
			if Direct:
			t1 = time.time()
			potdirect = MyTree.AllDirectPotential(eps)
			t2 = time.time()
			print "time = ",t2-t1


			sys.exit()

			################################
			# plot
			################################

			rp = nb.rxyz()

			################################
			# 1
			pt.subplot(2,1,1)
			pt.scatter(rp,pot)

			# analytic form
			r = sort(rp)
			pot = -1/sqrt(r2 + e2)
			pt.plot(r,pot)



			################################
			# 2
			if Direct:
			pt.subplot(2,1,2)
			pt.scatter(rp,potdirect)

			# analytic form
			r = sort(rp)
			pot = -1/sqrt(r2 + e2)
			pt.plot(r,pot)





			pt.show()