% Driver for the boundary integral chapter.
%
% The editable script supplies
%
%     names = method_names()
%     v = solve_bie(mid, n, a, nlobe, par, px, py)
%
% returning the computed solution at the points (px, py) inside the curve
% r(t) = 1 + a cos(nlobe t). n is the method's resolution: boundary nodes for
% the Nystrom method, source points for the method of fundamental solutions.
% par is the auxiliary parameter each method needs (the radius of the source
% circle, for the second one).
%
% The Dirichlet data is the trace of lp_bexact, so the exact solution is known
% everywhere inside and the error needs no reference computation.
p = lp_params();
a = p.a;
nlobe = p.nlobe;
names = method_names();
nm = numel(names);

out = struct();
out.names = names;

% the curve itself, for the pictures
tt = linspace(0, 2*pi, 401)';
C = lp_curve(a, nlobe, tt);
out.cx = C(:,1)';
out.cy = C(:,2)';

if strcmp(p.mode, 'field')
    % a Cartesian grid over the bounding box, masked to the interior
    m = 96;
    rmax = 1 + abs(a) + 0.05;
    xs = linspace(-rmax, rmax, m)';
    [X, Y] = meshgrid(xs, xs);
    th = atan2(Y, X);
    rr = sqrt(X.^2 + Y.^2);
    rb = 1 + a * cos(nlobe * th);
    inside = rr < rb - 1e-9;
    px = X(inside);
    py = Y(inside);

    t0 = tic;
    v = solve_bie(p.mid, p.n, a, nlobe, p.par, px, py);
    out.ms = toc(t0) * 1000;
    ve = lp_bexact(px, py);
    out.maxerr = max(abs(v - ve));

    U = NaN(m, m);
    E = NaN(m, m);
    U(inside) = v;
    E(inside) = abs(v - ve);
    out.x = xs';
    out.y = xs';
    out.U = lp_rows(U);
    out.Err = lp_rows(E);
    out.unknowns = p.n;
else
    % errors on a polar lattice: four rings, the last of them close enough to
    % the boundary to show what close evaluation costs
    nth = 60;
    th = 2*pi*(0:nth-1)'/nth;
    rb = 1 + a * cos(nlobe * th);
    rho = [0.25; 0.5; 0.75];
    px = [rho(1)*rb .* cos(th); rho(2)*rb .* cos(th); rho(3)*rb .* cos(th)];
    py = [rho(1)*rb .* sin(th); rho(2)*rb .* sin(th); rho(3)*rb .* sin(th)];
    ve = lp_bexact(px, py);
    % a ring at 98 per cent of the radius, where a layer potential is hardest
    % to evaluate accurately
    qx = 0.98 * rb .* cos(th);
    qy = 0.98 * rb .* sin(th);
    qe = lp_bexact(qx, qy);

    ns = p.ns(:)';
    nn = numel(ns);
    E = zeros(nm, nn);
    T = zeros(nm, nn);
    Q = zeros(nm, nn);
    EN = zeros(nm, nn);
    ET = drv_methodsb(1, NaN(3*nm, nn), ns, a, nlobe, p.par, px, py, ve, qx, qy, qe, nm);
    for k = 1:nm
        E(k, :) = ET(3*k-2, :);
        EN(k, :) = ET(3*k-1, :);
        T(k, :) = ET(3*k, :);
        Q(k, :) = lp_orders(E(k, :), ns);
    end
    out.Enear = lp_rows(EN);
    out.ns = ns;
    out.E = lp_rows(E);
    out.T = lp_rows(T);
    out.orders = lp_rows(Q);
end

lp_write(out);

function ET = drv_methodsb(k, ET, ns, a, nlobe, par, px, py, ve, qx, qy, qe, nm)
% Recursion rather than a loop, as in the other drivers: a call made from
% inside a loop to a routine that itself loops is compiled through numbl's slow
% path, and these are timed solves.
if k > nm
    return;
end
vw = solve_bie(k, 16, a, nlobe, par, px(1), py(1));    % warm-up
R = drv_resb(1, NaN(3, numel(ns)), k, ns, a, nlobe, par, px, py, ve, qx, qy, qe);
ET(3*k-2, :) = R(1, :);
ET(3*k-1, :) = R(2, :);
ET(3*k, :) = R(3, :);
ET = drv_methodsb(k+1, ET, ns, a, nlobe, par, px, py, ve, qx, qy, qe, nm);
end

function R = drv_resb(j, R, mid, ns, a, nlobe, par, px, py, ve, qx, qy, qe)
if j > numel(ns)
    return;
end
n = ns(j);
t0 = tic;
v = solve_bie(mid, n, a, nlobe, par, px, py);
el = toc(t0);
if el < 0.005
    el = drv_bestb(4, el, mid, n, a, nlobe, par, px, py);
end
q = solve_bie(mid, n, a, nlobe, par, qx, qy);
R(1, j) = max(abs(v - ve));
R(2, j) = max(abs(q - qe));
R(3, j) = el * 1000;
R = drv_resb(j+1, R, mid, ns, a, nlobe, par, px, py, ve, qx, qy, qe);
end

function el = drv_bestb(k, el, mid, n, a, nlobe, par, px, py)
if k == 0
    return;
end
t0 = tic;
v = solve_bie(mid, n, a, nlobe, par, px, py);
e2 = toc(t0);
if e2 < el
    el = e2;
end
el = drv_bestb(k-1, el, mid, n, a, nlobe, par, px, py);
end