test_Potential.py
021ca11fcae0master

test_Potential.py
View Options

			#!/usr/bin/env python

			from mpi4py import MPI

			from pNbody import ic

			from numpy import *
			from PyGadget import gadget

			import Ptools as pt

			import sys
			import time

			random.seed(MPI.COMM_WORLD.Get_rank()) # 2 for two points

			n = 10000
			rc = 2
			rmax = 100

			nb = ic.plummer(n,1,1,1,eps=rc,rmax=rmax,ftype='gadget')
			nb.rename('plummer.dat')
			nb.write()

			gadget.InitMPI() # init MPI
			gadget.InitDefaultParameters() # init default parameters
			gadget.Info()


			params = {}

			params['ErrTolTheta'] = 0.7

			params['DesNumNgb'] = 50
			params['MaxNumNgbDeviation'] = 1

			params['UnitLength_in_cm'] = 3.085e+21
			params['UnitMass_in_g'] = 4.435693e+44
			params['UnitVelocity_in_cm_per_s'] = 97824708.2699

			params['SofteningGas'] = rc/100.
			params['SofteningHalo'] = rc/100.
			params['SofteningDisk'] = rc/100.
			params['SofteningBulge'] = rc/100.
			params['SofteningStars'] = rc/100.
			params['SofteningBndry'] = rc/100.
			params['SofteningGasMaxPhys'] = rc/100.
			params['SofteningHaloMaxPhys'] = rc/100.
			params['SofteningDiskMaxPhys'] = rc/100.
			params['SofteningBulgeMaxPhys'] = rc/100.
			params['SofteningStarsMaxPhys'] = rc/100.
			params['SofteningBndryMaxPhys'] = rc/100.

			gadget.SetParameters(params)



			params = gadget.GetParameters()



			gadget.LoadParticles(array(nb.npart),nb.pos,nb.vel,nb.mass,nb.num,nb.tpe)

			r = arange(0,100,0.1)
			y = zeros(len(r))
			z = zeros(len(r))
			pos = transpose(array([r,y,z]))
			#pos = array([[0,0,0]])
			eps = rc/100.

			pot = gadget.Potential(pos,eps)



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

			val = pot


			pt.scatter(r,val,s=1)
			pt.show()