% 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 = NaN(nm, nn);
    T = NaN(nm, nn);
    Q = NaN(nm, nn);
    EN = NaN(nm, nn);
    for k = 1:nm
        vw = solve_bie(k, 16, a, nlobe, p.par, px(1), py(1));   % warm-up
        for j = 1:nn
            n = ns(j);
            % Keep the fastest of several runs: three at least whenever a
            % solve is cheap enough to allow it, many more when it is very
            % cheap, one when it is expensive. A single measurement of a
            % twenty-millisecond solve varies by a factor of two with whatever
            % else the machine is doing, and the numbers here are quoted in the
            % prose. Chapter 11 times itself the same way.
            el = Inf;
            tot = 0;
            rep = 0;
            while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
                t0 = tic;
                v = solve_bie(k, n, a, nlobe, p.par, px, py);
                dt = toc(t0);
                el = min(el, dt);
                tot = tot + dt;
                rep = rep + 1;
            end
            q = solve_bie(k, n, a, nlobe, p.par, qx, qy);
            E(k, j) = max(abs(v - ve));
            EN(k, j) = max(abs(q - qe));
            T(k, j) = el * 1000;
        end
        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);