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