% Fig for CE paper. Driver for layerpot.eval QBX slow (non-FMM)
% Barnett 9/25/12 - 10/3/12

clear; expt='f'; verb=1; figs=1; N = 130; %330 %120;  % base # pts
s = segment.smoothnonsym(N, 0.3, 0.3, 3);
zsing = [];%1+.5i; % data for expt=f: singularity loc, if empty, use entire func
intp = 't'; % t=trig poly (all k, standard); n=Nystrom (k=0 only; beta not odd)
d = domain(s, 1); k = 0; d.k = k;  % choose Laplace or Helmholtz
o = []; d.addlayerpot(s, 'D', o); l = d.bas{1};  % DLP
alsing = 2*pi*min(abs(imag(s.invertZparam(zsing)))); % alpha sing wrt (0,2pi)

dx = 0.01; bb = d.boundingbox;  % choose eval grid...
[xx yy] = meshgrid(bb(1):dx:bb(2),bb(3):dx:bb(4));
zz = xx(:)+1i*yy(:); tic; di = d.inside(zz); ii = find(di);
fprintf('inside %d pts, %d verts = %.3g s\n',numel(zz),numel(s.approxv),toc);

% fill tau density ... (choose expt)
if expt=='1'            % tau=1 case
  if k>0, error('Helmholtz has no tau=1 known potential!'); end
  f = @(z) -1 + 0*z;               % known solution -1 inside
  maxtau = 1;
  co = ones(N,1);                    % unit density
elseif expt=='f'         % solve a BVP
  if k==0    % Laplace data
    if ~isempty(zsing), f = @(z) log(abs(z-zsing));
    else, %f = @(z) real(z).*imag(z); % note: better since not zero along the ray 3*pi/4 that bump is on!
      f = @(z) cos(real(z)).*exp(+imag(z)); % note -imag(z) worse by 1 digit
    end % was: real(z).^2 - imag(z).^2;  % u(x,y) = x^2 - y^2 exact, Lap
  else      % Helmholtz data
    if ~isempty(zsing), f = @(z) 3*besselh(0,k*abs(z-zsing));  % u = pt src
    else, f = @(z) exp(1i*k*real(exp(1i*pi/3)'*z)); % u = plane wave
    end
  end
  maxtau = max(abs(f(s.x)));  % obviously this *isn't* tau max in strip!
  d.bas{1}.quad='m'; % Kress for self-eval (if kappa known & N<1400)
  %d.bas{1}.quad='a'; d.bas{1}.ord = 16; % Alpert for self-eval (if no kappa)
  p = bvp(d); s.setbc(-1, 'D', [], @(t) f(s.Z(t))); % BVP
  p.fillquadwei; p.sqrtwei=0*p.sqrtwei+1; % NB unit sqrtwei!
  tic; p.solvecoeffs; co=p.co; %opts.comparefunc = f; p.showsolution(opts);
  fprintf('BVP time fill+solve size %d lin sys = %.3g s\n',numel(co),toc)
  zp = .3-.1i;  % check at interior pt
  fprintf('err at distant pt = %.3g\n',p.pointsolution(pointset(zp))-f(zp))
  if verb>1, figure; semilogy(abs(fft(co)));hold on;title('\tau fourier decay');
    if ~isempty(zsing), m=1:N/2; plot(m,exp(-alsing*m),'r-'); end, drawnow; end
  if N<500, Nf = 3*N; p.updateN(Nf); % solve on finer grid to check
    p.sqrtwei = 1+0*p.sqrtwei; p.solvecoeffs;
    coa = p.co(2:3:end); p.updateN(N); % downsample back to N
    fprintf('est rel L2 err in tau = %.3g\n',norm(co-coa)/norm(coa))
  end
end
% native eval....
fz = nan*xx; fz(ii) = f(zz(ii)); fz = reshape(fz, size(xx)); % fill fz, u...
u = nan*xx;
if N<500   % switch for dense (crude eval matrix fill) or FMM
  tic; A = d.bas{1}.eval(pointset(zz(ii)));   % do grid eval by hand
  fprintf('time for filling interior di eval matrix A = %.3g s\n',toc);
  u(ii) = A*co; 
else
  tic; u(ii) = d.bas{1}.evalFMM(co,pointset(zz(ii)));
    fprintf('time for FMM interior di eval  = %.3g s\n',toc);
end
u = reshape(u, size(xx)); unat = u;
maxf = max(abs(fz(:)));  % max value of solution (scale to this)

if figs, figure; imagesc(xx(1,:),yy(:,1),log10(abs(unat-fz)/maxf)); % fig (a)
caxis([-16 0]); axis equal tight xy off; colormap(goodbw); colorbar;
%text(-1.3,1.1,'(a)','fontsize',14);%$$\log_{10} |u-f|$$','interpreter','latex')
set(gcf,'paperposition', [0 0 6 5]);
print -depsc2 rfn_k0_entire_N130_native.eps
end

% QBX geom params: ====
delta = 2.5;  % how many "h" away center is from boundary
al = delta/N; % imag strip width wrt [0,1]-scaled seg
Bdelta = 2.5; B = ceil(N/(2*Bdelta)); % how many boxes (Bdelta=boxheight/h)
als = 5.0/N; % imag strip width wrt [0,1]-scaled seg, for acc to 1e-13
% end params ======
t0 = (1:B)/B+1i*al; z0 = s.Z(t0);  % all preimages and centers
v = s.Z(s.t + 1i*als); % imag-shifted N bdry pts. Now clip off loopy ones:
vsd = abs(repmat(s.x,[1 N]) - repmat(v.',[N 1])); % find dists v to pO, O(N^2)
[vsdm, ivsd] = min(vsd,[],1); clear vsd; % get indices of nearest pts
oi.toz = s.t(ivsd);  % nearest on-pO initial guesses to help invertZ
vs = s.invertZparam(v,oi); % slow if wobbly curve
v = v(find(abs(imag(vs))>0.99*als)); % snip off too-close self-int loops
fc = ~utils.inpolywrapper(zz(ii), v); % for ii pts, true if in strip
ic = ii(find(fc)); zzc = zz(ic);      % pts to correct

% QBX convergence params:
%Ms = 2:2:18; bes = 2:.3:7;  % convergence sweeps (M=p, beta=Nf/N>1)
Ms = 10; bes=4;  % just test one param set. (Nyst interp: bes neq odd # !)
ei = nan(numel(Ms), numel(bes)); e2 = ei;  % rel error data (Linfty, L2)
for j=1:numel(bes), Nf = ceil(bes(j)*N); % ........... double loop over p, beta
  sf = utils.copy(s); sf.requadrature(Nf);  % fine segment
  lpf = layerpot(sf, 'D'); lpf.doms = d;    % fine layerpot
  if intp=='t', cof = quadr.interptrig(co, Nf); % fine grid density via trig interp - okay.
  elseif intp=='n',
    o.quad='m'; Af = layerpot.D(k,s,sf,o); % self-eval no jump coarse-to-fine
    % for coincident pts need to correct Nystrom matrix:
    for i=1:Nf, dd = min(abs(sf.t(i)-s.t)); di = find(abs(sf.t(i)-s.t)==dd);
      if dd<1e-15, Af(i,di) = -s.kappa(di)*s.speed(di)/(4*pi*N); end, end
    rhsf = f(sf.x); cof = -2*(rhsf-Af*co); % Nyst interp density
  end
  for i=1:numel(Ms), M = Ms(i);
    for b=1:B % loop over boxes, get indices in `box', localexp-eval
      [box ib] = l.QBXbox(real(t0(b)), [1/(2*B) als], zz);
      LEo.besselcode = 'm'; LEo.real = 0; LEo.rescale_rad = 1.5*s.dist(z0(b));
      LE = regfbbasis(z0(b), M, LEo); LE.doms = d;
      S2L = layerpot.localfromDLP(sf, LE); coLE = S2L * cof;
      if ~isempty(find(ib)), u(ib) = LE.eval(pointset(zz(ib))) * coLE; end
    end
    maxe = max(abs(u(ii)-fz(ii))); ei(i,j) = maxe/maxf;
    e2(i,j) = norm(u(ii)-fz(ii))/norm(fz(ii));
    fprintf('p=%d beta=%g: rel L-infty, L2 errs = %.3g %.3g\n',M,bes(j),ei(i,j),e2(i,j))
  end
end % .........................................
if numel(ei)>1, figure; imagesc(bes,Ms,log10(ei)); colormap(goodbw); colorbar;
  axis xy; xlabel('\beta'); ylabel('p'); %title('rel L-infty err');
  caxis([-14 0]);
  set(gcf,'paperposition', [0 0 6 5]);
  %print -depsc2 rfn_k0_entire_N130_qbx_conv.eps
end

if verb, figure; imagesc(xx(1,:),yy(:,1),log10(abs(u-fz)/maxf)); % rel err
  axis equal tight xy off; caxis([-16 log10(maxe/maxf)]);
colormap(goodbw); colorbar; if ~isempty(zsing), hold on; plot(zsing, 'k*'); end
b=B; o.nospill=1; box = l.QBXbox(real(t0(b)), [1/(2*B) als], [], o);
hold on; plot([box.v; box.v(1)],'-','color',[.5 .5 .5]);% show a box and center
plot(z0(b), '.', 'markersize',20, 'color', [.5 .5 .5]);
if figs, %text(-1.3,1.1,'(b)','fontsize',14);
  set(gcf,'paperposition', [0 0 6 5]);
  print -depsc2 rfn_k0_entire_N130_qbx.eps
end

end
if verb>2, figure; imagesc(xx(1,:),yy(:,1),real(u)); axis equal tight xy;
title('u'); end % potential

if 0, gam = nan(B,1); figure; s.plot; hold on;  % plot geom ........
for b=1:B
  box = l.QBXbox(real(t0(b)), [1/(2*B) 1/(1*B)]); % note not always 2x1 boxes
  R = max(abs(box.v - z0(b))); % max radius of box about z0
  gam(b) = s.QBXgammaconst(t0(b), R);
  plot([box.v; box.v(1)], 'b-'); plot(z0(b), 'k*');
end
figure; plot(1:B, gam, '+-'); title('gamma consts for z0''s');
end % ..........

