% Two ways to put the edge of a round drum on a square grid.
%
% The grid covers the square [-1, 1]^2 with spacing h = 2/N, and the drum is
% the unit disk.  Every node strictly inside the disk is an unknown; every
% other node is held at zero.  What differs is what the five-point Laplacian
% does at an unknown whose neighbour lies outside.
%
%   'staircase'  nothing: the outside neighbour is simply zero, as if the rim
%                ran through it.  The drum the grid sees is the union of the
%                cells around the inside nodes, a disk with a jagged edge, and
%                its radius is wrong by up to h.
%   'sw'         Shortley and Weller (1938): if the rim crosses the grid line
%                at a fraction theta of the way to the outside neighbour, use
%                the three-point second difference for unequal spacing,
%
%                  u_xx ~ 2/h^2 [ u_W / (1 + theta) - u_P / theta ]
%                         (with u = 0 at the crossing),
%
%                which is consistent to first order at the rim and gives
%                second-order frequencies.  The price is the 1/theta: a node
%                very close to the rim has a very large coefficient, and an
%                explicit time step has to respect it.
%
% The result is the stencil, times h^2, one number per node and direction:
% C at the node itself and E, W, Nn, S for its neighbours at (i, j+1), (i, j-1),
% (i+1, j), (i-1, j), where i counts y and j counts x.  `in` marks the
% unknowns.  The driver runs leapfrog with it.

function [C, E, W, Nn, S, in] = diskstencil(N, method)
h = 2 / N;
x = -1 + (0:N) * h;
[X, Y] = meshgrid(x, x);        % X(i, j) = x(j), Y(i, j) = x(i)
in = X .^ 2 + Y .^ 2 < 1 - 1e-12;
n1 = N + 1;
C = zeros(n1, n1);
E = zeros(n1, n1);
W = zeros(n1, n1);
Nn = zeros(n1, n1);
S = zeros(n1, n1);
sw = strcmp(method, 'sw');
for i = 2:n1 - 1
    for j = 2:n1 - 1
        if ~in(i, j)
            continue;
        end
        % x direction: distances to the east and west, as fractions of h
        te = 1;
        tw = 1;
        if ~in(i, j + 1) && sw
            te = (sqrt(1 - Y(i, j) ^ 2) - X(i, j)) / h;
        end
        if ~in(i, j - 1) && sw
            tw = (X(i, j) + sqrt(1 - Y(i, j) ^ 2)) / h;
        end
        tn = 1;
        ts = 1;
        if ~in(i + 1, j) && sw
            tn = (sqrt(1 - X(i, j) ^ 2) - Y(i, j)) / h;
        end
        if ~in(i - 1, j) && sw
            ts = (Y(i, j) + sqrt(1 - X(i, j) ^ 2)) / h;
        end
        % the unequal-spacing second difference in each direction; with every
        % theta equal to 1 this is the ordinary five-point stencil
        E(i, j) = 2 / (te * (te + tw)) * in(i, j + 1);
        W(i, j) = 2 / (tw * (te + tw)) * in(i, j - 1);
        Nn(i, j) = 2 / (tn * (tn + ts)) * in(i + 1, j);
        S(i, j) = 2 / (ts * (tn + ts)) * in(i - 1, j);
        C(i, j) = -2 / (te * tw) - 2 / (tn * ts);
    end
end
end