% Three explicit schemes for the square membrane.
%
%   'fd2'   leapfrog in time, the five-point Laplacian in space: second order
%           in both, stable for lam <= 1/sqrt(2).  Chapter 4's scheme.
%   'fd4'   leapfrog in time, the fourth-order Laplacian in space, whose 1-D
%           stencil is (-1, 16, -30, 16, -1)/12: fourth order in space but
%           still second order in time, and stable only for lam <= 0.612.
%   'fd44'  fourth order in both, by the modified-equation trick (Dablain
%           1986): the leapfrog update's leading time error is
%           (dt^2/12) u_tttt = (dt^2/12) c^4 Lap^2 u, so add it back, using
%           the five-point Laplacian applied twice for Lap^2, since the term
%           is already multiplied by dt^2.  One more stencil pass per step,
%           and the time step is 1/sqrt(2) again.
%
% The wide stencils need values one node beyond the rim.  With u = 0 on the
% rim, the wave equation makes every even derivative of u vanish there too,
% so the solution continues across the rim as an odd reflection, and that is
% exactly what the ghost nodes are filled with.  For the sine modes of the
% square the reflection is not an approximation at all.
%
% The contract is Chapter 4's: u and uold are (N+1)-by-(N+1), rim included;
% the result is the displacement at node (ir, jr) after each of M steps, and
% the whole field after each of the first K.

function [y, U] = drum2d(u, uold, lam, M, scheme, ir, jr, K)
n1 = size(u, 1);
g = n1 + 4;                       % two ghost layers on every side
P = zeros(g, g);
Q = zeros(g, g);
R = zeros(g, g);
W = zeros(g, g);
P(3:n1 + 2, 3:n1 + 2) = u;
Q(3:n1 + 2, 3:n1 + 2) = uold;
y = zeros(M, 1);
U = zeros(n1 * n1, K);
if strcmp(scheme, 'fd2')
    s = 0;
elseif strcmp(scheme, 'fd4')
    s = 1;
else
    s = 2;
end
l2 = lam ^ 2;
l4 = lam ^ 4 / 12;
a = n1 + 2;                       % the index of the far rim
for n = 1:M
    % odd reflection across the rim, two layers deep
    for j = 1:g
        P(2, j) = -P(4, j);
        P(1, j) = -P(5, j);
        P(a + 1, j) = -P(a - 1, j);
        P(a + 2, j) = -P(a - 2, j);
    end
    for i = 1:g
        P(i, 2) = -P(i, 4);
        P(i, 1) = -P(i, 5);
        P(i, a + 1) = -P(i, a - 1);
        P(i, a + 2) = -P(i, a - 2);
    end
    if s == 0
        for j = 4:a - 1
            for i = 4:a - 1
                R(i, j) = 2 * P(i, j) - Q(i, j) ...
                    + l2 * (P(i + 1, j) + P(i - 1, j) + P(i, j + 1) + P(i, j - 1) - 4 * P(i, j));
            end
        end
    else
        if s == 2
            % W = five-point Laplacian of P (times h^2), zero on the rim,
            % reflected oddly one layer out
            for j = 4:a - 1
                for i = 4:a - 1
                    W(i, j) = P(i + 1, j) + P(i - 1, j) + P(i, j + 1) + P(i, j - 1) - 4 * P(i, j);
                end
            end
            for j = 1:g
                W(2, j) = -W(4, j);
                W(a + 1, j) = -W(a - 1, j);
            end
            for i = 1:g
                W(i, 2) = -W(i, 4);
                W(i, a + 1) = -W(i, a - 1);
            end
        end
        for j = 4:a - 1
            for i = 4:a - 1
                L4 = (-P(i + 2, j) + 16 * P(i + 1, j) - 30 * P(i, j) + 16 * P(i - 1, j) - P(i - 2, j) ...
                      - P(i, j + 2) + 16 * P(i, j + 1) - 30 * P(i, j) + 16 * P(i, j - 1) - P(i, j - 2)) / 12;
                v = 2 * P(i, j) - Q(i, j) + l2 * L4;
                if s == 2
                    v = v + l4 * (W(i + 1, j) + W(i - 1, j) + W(i, j + 1) + W(i, j - 1) - 4 * W(i, j));
                end
                R(i, j) = v;
            end
        end
    end
    Q = P;
    P = R;
    y(n) = P(ir + 2, jr + 2);
    if n <= K
        U(:, n) = reshape(P(3:a, 3:a), n1 * n1, 1);
    end
end
end