1% This script tests a finite difference scheme on the 3D coulomb 2% convolution problem. 3 4clear all; 5close all; 6 7% Size of cube 8nsize = 20; 9% Dimensions of cube 10L = 5.0; 11% points in real space 12x = linspace(-L/2,L/2,nsize); 13y = linspace(-L/2,L/2,nsize); 14z = linspace(-L/2,L/2,nsize); 15% make rho 16fprintf(1,'making rho ...\n\n'); 17rho = zeros(nsize,nsize,nsize); 18for ri = 1:nsize 19 for rj = 1:nsize 20 for rk = 1:nsize 21 rho(ri,rj,rk) = exp(-5*(x(ri)*x(ri) + y(rj)*y(rj) + z(rk)*z(rk))); 22 end 23 end 24end 25% compute indicies 26[I,J,K] = ind2sub([nsize nsize nsize], [1:nsize*nsize*nsize]); 27 28% return value 29rv = zeros(nsize,nsize,nsize); 30 31% main loop 32% if you think of this convolution as a matrix-vector multiply, 33% mi and mj index into the matrix which has the size of [npoints, npoints] 34% where npoints is the total number of points in the cube 35fprintf(1,'running main loop ...\n\n'); 36for mi = 1:nsize*nsize*nsize 37 fprintf(1,'mi = %d\n', mi); 38 for mj = 1:nsize*nsize*nsize 39 % compute indices to express mi and mj into a 3-D index 40 % so we can associate ki and kj with variables x and x' 41 ki = [I(mi) J(mi) K(mi)]; 42 kj = [I(mj) J(mj) K(mj)]; 43 % compute distance 44 distx = x(ki(1)) - x(kj(1)); 45 disty = x(ki(2)) - x(kj(2)); 46 distz = x(ki(3)) - x(kj(3)); 47 dist = sqrt(distx^2 + disty^2 + distz^2); 48 temp = rho(kj(1), kj(2), kj(3))/dist; 49 rv(ki) = rv(ki) + temp; 50 end 51end 52