% Leapfrog in a rigid rectangular room: the seven-point Laplacian in three
% dimensions, stable for lam <= 1/sqrt(3).
%
% The grid has n(1) by n(2) by n(3) nodes, with the walls on the outermost
% nodes.  A rigid wall means du/dn = 0, which is imposed by mirroring the
% first node inside each wall into a ghost node just outside it.  That
% mirror is the discrete version of the method of images: the grid's
% solution in the room is exactly the grid's free-space solution for the
% source and all its mirror images, so the only error left is the scheme's
% own dispersion.
%
% The room starts silent.  At every step the source node receives src(n),
% already scaled by the driver to dt^2 f(t_n) / h^3, which is a point source
% of strength f(t) spread over one cell.  The result is the pressure at the
% receiver node after each step.
%
% The field is stored as one long vector and the neighbours of node q are
% q +- 1, q +- sy and q +- sz.  numbl runs this form about seven times faster
% than the same loop written with three subscripts.

function y = room3d(n, lam, M, is, ir, src)
nx = n(1) + 2;
ny = n(2) + 2;
nz = n(3) + 2;
sy = nx;
sz = nx * ny;
u = zeros(nx * ny * nz, 1);
up = u;
un = u;
l2 = lam ^ 2;
qs = is(1) + 1 + is(2) * sy + is(3) * sz;       % padded, so every index is one further in
qr = ir(1) + 1 + ir(2) * sy + ir(3) * sz;
y = zeros(M, 1);
for step = 1:M
    % ghosts: the mirror image of the first interior node, on all six faces
    for c = 1:nz
        for j = 1:ny
            b = (c - 1) * sz + (j - 1) * sy;
            u(b + 1) = u(b + 3);
            u(b + nx) = u(b + nx - 2);
        end
    end
    for c = 1:nz
        for i = 1:nx
            b = (c - 1) * sz + i;
            u(b) = u(b + 2 * sy);
            u(b + (ny - 1) * sy) = u(b + (ny - 3) * sy);
        end
    end
    for j = 1:ny
        for i = 1:nx
            b = (j - 1) * sy + i;
            u(b) = u(b + 2 * sz);
            u(b + (nz - 1) * sz) = u(b + (nz - 3) * sz);
        end
    end
    for c = 2:nz - 1
        for j = 2:ny - 1
            b = (c - 1) * sz + (j - 1) * sy;
            for i = 2:nx - 1
                q = b + i;
                un(q) = 2 * u(q) - up(q) + l2 * (u(q + 1) + u(q - 1) + u(q + sy) + u(q - sy) ...
                    + u(q + sz) + u(q - sz) - 6 * u(q));
            end
        end
    end
    un(qs) = un(qs) + src(step);
    up = u;
    u = un;
    y(step) = u(qr);
end
end