% 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), entry by entry. A
% double loop is about twice as fast here as the same expression written with
% implicit expansion, and several times faster than forming the coordinate
% differences with outer products against ones(1,n); it also writes each
% entry once instead of passing over four n-by-n temporaries.
A = zeros(n, n);
c = h / (2*pi);
for j = 1:n
    xj = x(j);
    yj = y(j);
    nxj = nx(j);
    nyj = ny(j);
    sj = sp(j);
    for i = 1:n
        dx = x(i) - xj;
        dy = y(i) - yj;
        A(i, j) = c * (nxj*dx + nyj*dy) / (dx*dx + dy*dy) * sj;
    end
    % the diagonal limit, not a singularity, plus the jump relation's -1/2
    A(j, j) = -0.5 - c * kappa(j) * sj / 2;
end

sigma = A \ lp_bexact(x, y);

% evaluate the layer potential at the target points, by the same loop
nt = numel(px);
v = zeros(nt, 1);
for k = 1:nt
    xk = px(k);
    yk = py(k);
    s = 0;
    for j = 1:n
        dx = xk - x(j);
        dy = yk - y(j);
        s = s + (nx(j)*dx + ny(j)*dy) / (dx*dx + dy*dy) * sp(j) * sigma(j);
    end
    v(k) = c * s;
end
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