% Driver for Chapter 7: letting sound out of a grid.
%
% The editable script supplies
%
%     U = openbox(u0, uold, lam, M, method, W, keep)
%
% which runs leapfrog on a square region with some treatment of its edges,
% and returns the region's field at the nodes `keep` after every step.
%
% The reference is exact without being a formula.  An explicit scheme moves
% information by at most one node per step, so the same leapfrog scheme on a
% grid that extends M nodes further in every direction cannot feel its own
% edges anywhere in the region during M steps: inside the region, it is the
% scheme's free-space solution, to rounding error.  The driver runs that
% reference with its own code.  Whatever differs from it in the region is
% what the edge treatment sent back, and nothing else.
%
% mode 'snap'   one method: snapshots, and the reflection it produced.
% mode 'sweep'  every method at several layer widths: reflection against
%               cost.
% mode 'hear'   a longer run with a click, heard at one point, for a rigid
%               box and for two absorbing edges.
p = wv_params();
out = struct();
N = p.N;
n0 = N + 1;
lam = p.lam;
h = p.L / N;
c = 343;
dt = lam * h / c;
x = (0:N)' * h;
[X, Y] = meshgrid(x, x);
sg = p.sigma;
u0 = exp(-((X - p.xs) .^ 2 + (Y - p.ys) .^ 2) / (2 * sg ^ 2));
% released from rest: uold = u0 + (lam^2/2) Lap_h u0, the same for every run
lap = zeros(n0, n0);
lap(2:N, 2:N) = u0(3:n0, 2:N) + u0(1:N - 1, 2:N) + u0(2:N, 3:n0) + u0(2:N, 1:N - 1) - 4 * u0(2:N, 2:N);
um = u0 + 0.5 * lam ^ 2 * lap;

if strcmp(p.mode, 'hear')
    M = round(p.T / dt);
    ir = round(p.xr / h) + 1;
    jr = round(p.yr / h) + 1;
    meths = p.methods;
    Yh = zeros(M, numel(meths));
    for k = 1:numel(meths)
        U = openbox(u0, um, lam, M, meths{k}, p.W, (jr - 1) * n0 + ir);
        Yh(:, k) = U(1, :)';
    end
    s = max(abs(Yh(:)));
    fs = 1 / dt;
    clips = cell(1, numel(meths));
    for k = 1:numel(meths)
        clips{k} = wv_audio(Yh(:, k), fs, s);
    end
    out.clips = clips;
    out.methods = meths;
    out.fs = fs;
    out.t = (1:M) * dt * 1000;
    nkeep = min(M, round(0.08 / dt));
    out.tt = out.t(1:nkeep);
    out.traces = wv_rows(Yh(1:nkeep, :)' / s);
else
    M = round(p.T / dt);
    Uref = drv_reference(u0, um, lam, M);
    rn = sum(Uref(:) .^ 2);
    if strcmp(p.mode, 'snap')
        t0 = tic;
        U = openbox(u0, um, lam, M, p.method, p.W, (1:n0 * n0)');
        out.ms = toc(t0) * 1000;
        E = U - Uref;
        out.refl = 10 * log10(sum(E(:) .^ 2) / rn);
        ks = min(M, round(p.tsnap / dt));
        out.tsnap = ks * dt * 1000;
        out.ref = wv_rows(reshape(Uref(:, ks), n0, n0));
        out.num = wv_rows(reshape(U(:, ks), n0, n0));
        out.err = wv_rows(reshape(E(:, ks), n0, n0));
        out.emax = max(abs(E(:)));
        out.umax = max(abs(Uref(:, ks)));
        out.urefmax = max(abs(Uref(:)));
        % the reflected energy in the region, step by step
        out.t = (1:M) * dt * 1000;
        out.eng = 10 * log10(max(sum(E .^ 2, 1) / max(sum(Uref .^ 2, 1)), 1e-30));
        out.L = p.L;
    else
        meths = p.methods;
        Ws = p.Ws(:)';
        R = NaN(numel(meths), numel(Ws));
        TM = NaN(numel(meths), numel(Ws));
        for k = 1:numel(meths)
            for j = 1:numel(Ws)
                if (strcmp(meths{k}, 'rigid') || strcmp(meths{k}, 'mur')) && j > 1
                    continue;
                end
                t0 = tic;
                U = openbox(u0, um, lam, M, meths{k}, Ws(j), (1:n0 * n0)');
                TM(k, j) = toc(t0) * 1000;
                E = U - Uref;
                R(k, j) = 10 * log10(sum(E(:) .^ 2) / rn);
            end
        end
        out.methods = meths;
        out.Ws = Ws;
        out.refl = wv_rows(R);
        out.ms = wv_rows(TM);
    end
end

wv_write(out);

function U = drv_reference(u0, uold, lam, M)
% The same scheme on a grid M + 2 nodes larger on every side, with its far
% edges held at zero: nothing from those edges reaches the region in M steps.
n0 = size(u0, 1);
e = M + 2;
n1 = n0 + 2 * e;
u = zeros(n1, n1);
uo = zeros(n1, n1);
u(e + 1:e + n0, e + 1:e + n0) = u0;
uo(e + 1:e + n0, e + 1:e + n0) = uold;
un = zeros(n1, n1);
l2 = lam ^ 2;
U = zeros(n0 * n0, M);
for n = 1:M
    % only the nodes that can have become nonzero need updating
    lo = max(2, e + 1 - n);
    hi = min(n1 - 1, e + n0 + n);
    for j = lo:hi
        for i = lo:hi
            un(i, j) = 2 * u(i, j) - uo(i, j) + l2 * (u(i + 1, j) + u(i - 1, j) + u(i, j + 1) + u(i, j - 1) - 4 * u(i, j));
        end
    end
    uo = u;
    u = un;
    U(:, n) = reshape(u(e + 1:e + n0, e + 1:e + n0), n0 * n0, 1);
end
end