% The image-source method for a rectangular room (Allen and Berkley 1979).
%
% A wall that reflects sound specularly behaves, for the listener, like a
% window onto a mirror-image room with a mirror-image source in it.  Mirror
% the mirror rooms in their turn and the whole of space fills with a lattice
% of image sources, eight to each cell of size 2 Lx by 2 Ly by 2 Lz:
%
%   x = (1 - 2p) xs + 2 l Lx,   p in {0, 1}, l any integer,
%
% and the same in y and z.  The path from image to listener crosses the
% walls x = 0 and x = Lx |l - p| and |l| times, so the image's pressure is
%
%   beta^(number of reflections) / (4 pi d),   arriving d / c after the source,
%
% with beta the pressure reflection factor of the walls.  For rigid walls,
% beta = 1, this is not a model: it is the exact solution of the wave
% equation in the room, by the method of images.  For beta < 1 it is a model,
% one in which every wall reflects the same fraction at every angle and every
% frequency, which no real wall does.
%
% Each arrival is a delta at a time that falls between samples.  With taps = 0
% it is put on the nearest sample, as Allen and Berkley did; otherwise it is
% spread over `taps` samples as a Hann-windowed sinc, which is the delta
% band-limited to the Nyquist frequency, cut off at the window.
%
% L = [Lx Ly Lz], xs and xr the source and the receiver; the result is the
% impulse response, sampled at fs, up to tmax seconds.

function h = images(L, xs, xr, beta, fs, tmax, taps)
% the two ways of placing an arrival are two separate loops, because numbl
% will not compile a loop in which a variable is assigned on only one branch
% of an `if`, and a loop nest with the choice inside it is such a loop
if taps == 0
    h = img_nearest(L, xs, xr, beta, fs, tmax);
else
    h = img_sinc(L, xs, xr, beta, fs, tmax, taps);
end
end

function h = img_nearest(L, xs, xr, beta, fs, tmax)
c = 343;
rmax = c * tmax;
h = zeros(ceil(tmax * fs) + 2, 1);
nl = ceil(rmax / (2 * L(1))) + 1;
nm = ceil(rmax / (2 * L(2))) + 1;
nk = ceil(rmax / (2 * L(3))) + 1;
% beta^n from a table: numbl compiles an indexed read, but not a power whose
% exponent it cannot prove is a whole number
bp = beta .^ (0:2 * (nl + nm + nk) + 6);
for p = 0:1
    for q = 0:1
        for r = 0:1
            for l = -nl:nl
                dx = (1 - 2 * p) * xs(1) + 2 * l * L(1) - xr(1);
                ax = abs(l - p) + abs(l);
                for m = -nm:nm
                    dy = (1 - 2 * q) * xs(2) + 2 * m * L(2) - xr(2);
                    ay = abs(m - q) + abs(m);
                    for k = -nk:nk
                        dz = (1 - 2 * r) * xs(3) + 2 * k * L(3) - xr(3);
                        d = real(sqrt(dx * dx + dy * dy + dz * dz));
                        if d <= rmax
                            a = bp(ax + ay + abs(k - r) + abs(k) + 1) / (4 * pi * d);
                            jn = floor(d / c * fs + 0.5) + 1;
                            h(jn) = h(jn) + a;
                        end
                    end
                end
            end
        end
    end
end
end

function h = img_sinc(L, xs, xr, beta, fs, tmax, taps)
c = 343;
rmax = c * tmax;
h = zeros(ceil(tmax * fs) + taps + 2, 1);
nl = ceil(rmax / (2 * L(1))) + 1;
nm = ceil(rmax / (2 * L(2))) + 1;
nk = ceil(rmax / (2 * L(3))) + 1;
bp = beta .^ (0:2 * (nl + nm + nk) + 6);
half = taps / 2;
for p = 0:1
    for q = 0:1
        for r = 0:1
            for l = -nl:nl
                dx = (1 - 2 * p) * xs(1) + 2 * l * L(1) - xr(1);
                ax = abs(l - p) + abs(l);
                for m = -nm:nm
                    dy = (1 - 2 * q) * xs(2) + 2 * m * L(2) - xr(2);
                    ay = abs(m - q) + abs(m);
                    for k = -nk:nk
                        dz = (1 - 2 * r) * xs(3) + 2 * k * L(3) - xr(3);
                        d = real(sqrt(dx * dx + dy * dy + dz * dz));
                        if d > rmax
                            continue;
                        end
                        a = bp(ax + ay + abs(k - r) + abs(k) + 1) / (4 * pi * d);
                        tau = d / c * fs;
                        j0 = floor(tau);
                        for j = max(0, j0 - half + 1):j0 + half
                            x = j - tau;
                            if abs(x) < 1e-12
                                s = 1;
                            else
                                s = sin(pi * x) / (pi * x);
                            end
                            h(j + 1) = h(j + 1) + a * s * (0.5 + 0.5 * cos(2 * pi * x / taps));
                        end
                    end
                end
            end
        end
    end
end
end