% Driver for Chapter 9: a room, by finite differences.
%
% The editable script supplies
%
%     y = room3d(n, lam, M, is, ir, src)
%
% leapfrog in a rigid box on a grid of n(1) by n(2) by n(3) nodes, driven at
% node is and heard at node ir.
%
% The source is a point source of strength f(t), the derivative of a
% Gaussian, which puts no net volume of air into the box.  For rigid walls
% the exact pressure is the image sum
%
%   u(t) = sum over images  f(t - d/c) / (4 pi c^2 d),
%
% which the driver evaluates directly; Chapter 8's image method with beta = 1
% is the same formula, and here it is not a model but the solution.
%
% mode 'converge'  the error at the receiver against the grid, and the time
%                  per node per step, from which the driver extrapolates what
%                  a full-band simulation would cost.
% mode 'hear'      a longer run on two grids, against the exact response.
p = wv_params();
out = struct();
c = 343;
L = p.L(:)';
lam = p.lam;
sg = p.sigma;
t0p = 5 * sg;

if strcmp(p.mode, 'converge')
    Ns = p.Ns(:)';
    E = NaN(1, numel(Ns));
    TM = NaN(1, numel(Ns));
    NS = NaN(1, numel(Ns));
    cells = NaN(1, numel(Ns));
    for k = 1:numel(Ns)
        [y, ye, dt, ms, n] = drv_run(Ns(k), L, lam, p.T, p.xs, p.xr, sg, t0p, c);
        E(k) = sqrt(sum((y - ye) .^ 2) / sum(ye .^ 2));
        TM(k) = ms;
        cells(k) = prod(n);
        NS(k) = ms * 1e6 / (prod(n) * numel(y));      % nanoseconds per node per step
        if k == numel(Ns)
            out.t = (1:numel(y)) * dt * 1000;
            out.y = y';
            out.ye = ye';
        end
        if k == 1
            out.t1 = (1:numel(y)) * dt * 1000;
            out.y1 = y';
        end
    end
    out.Ns = Ns;
    out.err = E;
    out.ms = TM;
    out.cells = cells;
    out.ns = NS;
    out.orders = wv_orders(E, Ns);
    % points per wavelength at the pulse's upper frequency, taken as the
    % point where its spectrum is 40 dB below its peak
    f = linspace(1, 20000, 20000);
    S = f .* exp(-0.5 * (2 * pi * f * sg) .^ 2);
    S = S / max(S);
    fhi = f(find(S > 0.01, 1, 'last'));
    out.fhi = fhi;
    out.ppw = c / fhi ./ (L(1) ./ Ns);
    % what a full-band run of Chapter 8's room would cost at the measured
    % rate: 10 points per shortest wavelength at 20 kHz, one second of sound.
    % The rate is the best of the grids from N = 32 up, where the fixed costs
    % no longer matter and a single slow run (a garbage collection, another
    % tab) does not decide the answer.
    rate = min(NS(3:end));
    out.rate = rate;
    V8 = 81;
    h8 = c / 20000 / 10;
    dt8 = lam * h8 / c;
    upd = V8 / h8 ^ 3 / dt8;
    out.fullband = upd * rate * 1e-9;              % seconds
    out.fullcells = V8 / h8 ^ 3;
    h2 = c / 2000 / 10;                             % the same up to 2 kHz only
    out.lowband = V8 / h2 ^ 3 / (lam * h2 / c) * rate * 1e-9;
else
    Ns = p.Ns(:)';
    fs = p.fs;
    te = (1:round(p.T * fs))' / fs;
    ye = drv_exact(te, L, p.xs, p.xr, sg, t0p, c);
    sc = max(abs(ye));
    out.clipExact = wv_audio(ye, fs, sc);
    clips = cell(1, numel(Ns));
    for k = 1:numel(Ns)
        [y, ~, dt, ~, ~] = drv_run(Ns(k), L, lam, p.T, p.xs, p.xr, sg, t0p, c, true);
        clips{k} = wv_audio(y, 1 / dt, sc);
        fsk(k) = 1 / dt;
    end
    out.clips = clips;
    out.Ns = Ns;
    out.fsk = fsk;
    nw = round(0.25 * fs);
    [fr, db] = wv_spectrum(ye(1:nw), fs, 3000, 600);
    out.fr = fr;
    out.dbe = db;
end
wv_write(out);

function [y, ye, dt, ms, n] = drv_run(N, L, lam, T, xs, xr, sg, t0p, c, noexact)
h = L(1) / N;
n = round(L / h) + 1;
dt = lam * h / c;
M = round(T / dt);
t = (1:M)' * dt;
% step k computes u at t(k) from the two levels before it, so the source term
% it adds is the one at the middle level, f(t(k) - dt)
ts = t - dt;
f = drv_pulse(ts, sg, t0p);
src = dt ^ 2 * f / h ^ 3;
is = round(xs / h) + 1;
ir = round(xr / h) + 1;
room3d(n, lam, 3, is, ir, src(1:3));             % warm-up
tic0 = tic;
y = room3d(n, lam, M, is, ir, src);
ms = toc(tic0) * 1000;
if nargin > 9 && noexact
    ye = [];
else
    ye = drv_exact(t, L, xs, xr, sg, t0p, c);
end
end

function f = drv_pulse(t, sg, t0p)
% the derivative of a Gaussian: no DC, so no permanent pressure in the box
s = (t - t0p) / sg;
f = -s .* exp(-0.5 * s .^ 2);
end

function y = drv_exact(t, L, xs, xr, sg, t0p, c)
% the image sum for rigid walls, every image close enough to be heard by
% the end of t
y = zeros(size(t));
tend = t(end);
rmax = c * (tend + 0.001);
nl = ceil(rmax ./ (2 * L)) + 1;
dtt = t(2) - t(1);
w = 7 * sg;                                    % the pulse is below 1e-10 beyond this
for p = 0:1
    for q = 0:1
        for r = 0:1
            for l = -nl(1):nl(1)
                dx = (1 - 2 * p) * xs(1) + 2 * l * L(1) - xr(1);
                for m = -nl(2):nl(2)
                    dy = (1 - 2 * q) * xs(2) + 2 * m * L(2) - xr(2);
                    for k = -nl(3):nl(3)
                        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
                        tc = d / c + t0p;
                        i0 = max(1, floor((tc - w - t(1)) / dtt) + 1);
                        i1 = min(numel(t), ceil((tc + w - t(1)) / dtt) + 1);
                        a = 1 / (4 * pi * c * c * d);
                        for i = i0:i1
                            s = (t(i) - tc) / sg;
                            y(i) = y(i) - a * s * exp(-0.5 * s * s);
                        end
                    end
                end
            end
        end
    end
end
end