nbody-convol-0.py
800e93ab2944master

nbody-convol-0.py
View Options

			#!/usr/bin/env python


			from pNbody import *
			from pNbody import ic
			from pNbody.nbodymodule import *
			import Ptools as pt


			n = 100000
			eps = 1
			Rmax = 10
			nb = ic.plummer(n,1,1,1,eps,Rmax)

			# compute density between 0-N-1
			N = 128
			rmin = -10.
			rmax = 10.

			shape = (N,)
			val = ones(nb.pos.shape).astype(float32)
			mass = nb.mass.astype(float32)
			r = nb.x()
			r = (r-rmin)/(rmax-rmin)
			r = r.astype(float32)

			mr = mkmap1d(r,mass,val,shape).astype(float)


			# set the kernel

			e = 1e-2
			e2 = e*e
			def Kernel(r):
			return -1/sqrt(e2+r*r)


			# set the grid
			K = 2*N
			M = zeros(2*K,float)
			G = zeros(2*K,float) # for the "manual convolution"


			# set M
			for i in xrange(-N,N):
			if i<0:
			M[i+N+N] = 0
			else:
			M[i+N+N]=mr[i]

			for i in xrange(0,N):
			g = Kernel(i)

			G[i+K] = g # 0 -> N-1
			G[i] = g # -2N -> -N-1
			G[K-1- i] = g # -1 -> -N
			G[2*K-1- i] = g # 2N-1 -> N


			# do the convolution
			Phi = zeros(2*K,float)

			for i in xrange(N,3*N): # -N, N (physique)
			for j in xrange(N,3*N): # -N, N (physique)
			Phi[i] = Phi[i] + G[i-j] * M[j]




			pt.subplot(3,1,1)
			pt.plot(M)
			pt.subplot(3,1,2)
			pt.plot(G)
			pt.subplot(3,1,3)
			pt.plot(Phi)
			pt.show()