% Two boundary methods for the interior Dirichlet problem on a smooth curve.
%
% 1. NYSTROM, double layer. Look for the solution as the potential of a dipole
%    density sigma spread over the curve,
%
%        u(x) = (1/2pi) INT_Gamma  (n_y . (x - y)) / |x - y|^2  sigma(y) ds(y),
%
%    which is harmonic inside for any sigma. Letting x approach the curve from
%    inside picks up a jump, and the boundary condition becomes a Fredholm
%    equation of the second kind,
%
%        -sigma(x)/2 + INT_Gamma K(x, y) sigma(y) ds(y) = g(x).
%
%    Second kind is the important phrase: the operator is identity plus
%    compact, its eigenvalues cluster away from zero, and the discrete system
%    is well conditioned no matter how fine the discretisation. Compare the
%    finite difference matrices, whose condition number grows like N^2.
%
%    The kernel K looks singular but is not: as y approaches x it tends to
%    -kappa(x)/(4pi), where kappa is the curvature. So the integrand is smooth
%    and periodic, and the plain trapezoid rule, which is spectrally accurate
%    for smooth periodic integrands, converges geometrically. That single fact
%    is what makes this method fast.
%
% 2. FUNDAMENTAL SOLUTIONS. Place point sources on a circle of radius par
%    outside the curve and fit their combined field to the boundary data by
%    least squares. Every basis function is harmonic, so the equation is
%    satisfied exactly and only the boundary condition is approximated. It is
%    three lines of code and it also converges geometrically, but the matrix is
%    exponentially ill conditioned in the number of sources, and the accuracy
%    it reaches depends on where the sources are put.

function names = method_names()
names = {'double-layer Nystrom', 'fundamental solutions'};
end

function v = solve_bie(mid, n, a, nlobe, par, px, py)
if mid == 1
    v = nystrom(n, a, nlobe, px, py);
else
    v = mfs(n, a, nlobe, par, px, py);
end
end

function v = nystrom(n, a, nlobe, px, py)
t = 2*pi*(0:n-1)' / n;
h = 2*pi / n;
C = lp_curve(a, nlobe, t);
x = C(:,1);
y = C(:,2);
xp = C(:,3);
yp = C(:,4);
xpp = C(:,5);
ypp = C(:,6);

sp = sqrt(xp.^2 + yp.^2);            % speed |gamma'(t)|
nx = yp ./ sp;                       % outward normal, counterclockwise curve
ny = -xp ./ sp;
kappa = (xp .* ypp - yp .* xpp) ./ sp.^3;

% the Nystrom matrix: A = -I/2 + h * K * diag(speed)
DX = x * ones(1, n) - ones(n, 1) * x';
DY = y * ones(1, n) - ones(n, 1) * y';
R2 = DX.^2 + DY.^2;
K = ((ones(n, 1) * nx') .* DX + (ones(n, 1) * ny') .* DY) ./ (2*pi*R2);
K(1:n+1:end) = -kappa / (4*pi);      % the diagonal limit, not a singularity
A = -0.5*eye(n) + h * K .* (ones(n, 1) * sp');

sigma = A \ lp_bexact(x, y);

% evaluate the layer potential at the target points
np = numel(px);
EX = px(:) * ones(1, n) - ones(np, 1) * x';
EY = py(:) * ones(1, n) - ones(np, 1) * y';
KE = ((ones(np, 1) * nx') .* EX + (ones(np, 1) * ny') .* EY) ./ (2*pi*(EX.^2 + EY.^2));
v = (h * KE .* (ones(np, 1) * sp')) * sigma;
end

function v = mfs(n, a, nlobe, par, px, py)
% n sources on a circle of radius par, fitted at 3n boundary points.
ts = 2*pi*(0:n-1)' / n;
sx = par * cos(ts);
sy = par * sin(ts);

nb = 3 * n;
tb = 2*pi*(0:nb-1)' / nb;
C = lp_curve(a, nlobe, tb);
bx = C(:,1);
by = C(:,2);

A = source_matrix(bx, by, sx, sy);
c = A \ lp_bexact(bx, by);
v = source_matrix(px(:), py(:), sx, sy) * c;
end

function A = source_matrix(x, y, sx, sy)
% log|z - z_j| for each source, plus a constant column
np = numel(x);
n = numel(sx);
A = zeros(np, n + 1);
for j = 1:n
    A(:, j) = 0.5 * log((x - sx(j)).^2 + (y - sy(j)).^2);
end
A(:, n+1) = 1;
end