nbody-convol-0.py
No OneTemporary
Actions

Subscribers

None

File Metadata

	#!/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()