% Check pointwise evaluation error of native & close-eval quadrature scheme in
% BIE Laplace & Helmholtz
% Barnett 3/15/12

clear; expt='f'; N = 60;   % choose expt and base # bdry pts
%s = segment.smoothstar(N, 0.3, 3, 0.2); % ......... choose shape........
s = segment.smoothnonsym(N, 0.3, 0.3, 3);
%a=.3;b=.3;w=3; a1=0.1; a2=0.1; t0=3*pi/4; % smoothnonsym w/ gaussian bump @ t0
%s = segment.radialfunc(N, {@(q) 1 + a*cos(w*(q+b*cos(q))) - a1*exp(-0.5*((q-t0)/a2).^2), @(q) -a*sin(w*(q+b*cos(q))).*w.*(1-b*sin(q)) + (a1/a2^2)*(q-t0).*exp(-0.5*((q-t0)/a2).^2) }); %don't need curvatue yet since Alpert acc ok given bump
% .......................................................................

d = domain(s, 1); k = 0; d.k = k;  % choose Laplace or Helmholtz
o = []; d.addlayerpot(s, 'D', o); % d.plot;
dx = 0.01; bb = d.boundingbox;[xx yy] = meshgrid(bb(1):dx:bb(2),bb(3):dx:bb(4));
zz = xx(:)+1i*yy(:); di = d.inside(zz); ii = find(di);
A = d.bas{1}.eval(pointset(zz(ii)));   % do grid eval by hand

if expt=='1'
  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'
  s.requadrature(3*N); % solve on finer grid
  p = bvp(d);  % now solve a BVP to get a more interesting density
  d.bas{1}.quad='a'; d.bas{1}.ord = 16; % Alpert for self-eval (if no kappa)
  if k==0, f = @(z) real(z).*imag(z); % note: better since not zero along the ray 3*pi/4 that bump is on!
    % was: real(z).^2 - imag(z).^2;  % u(x,y) = x^2 - y^2 exact, Lap
  else, f = @(z) exp(1i*k*real(exp(1i*pi/3)'*z)); end % u = plane wave
  maxtau = max(abs(f(s.x)));  % obviously this *isn't* tau max in strip!
  s.setbc(-1, 'D', [], @(t) f(s.Z(t)));
  p.solvecoeffs; %opts.comparefunc = f; p.showsolution(opts);
  co = p.co(2:3:end-1);  % extract the coincident nodes (3N->N downsample)
  s.requadrature(N); % for plotting later
end

fz = nan*xx; fz(ii) = f(zz(ii)); fz = reshape(fz, size(xx)); % fill fz, u...
u = nan*xx; u(ii) = A*co; u = reshape(u, size(xx));

% usual color image:
figure; imagesc(xx(1,:),yy(:,1),log10(abs(u-fz))); colorbar; caxis([-16 0]);
axis equal; set(gca,'ydir','normal'); axis tight;
%stop

% contour comparison:
figure; errs = 10.^(-16:2:0); %10.^(-16:2:0);
contourf(xx,yy,log10(abs(u-fz)), log10(errs)); caxis([-16 0]);
colorbar; axis equal; hold on; plot([s.x; s.x(1)], 'k.-');
taus = log(maxtau./errs)/N; %[0:.05:.3]; % use complex map imag shifts
o = []; o.levels = taus/(2*pi);
if exist('t0') % hack to get bump to plot imag dist ok
  n=10; no=10; o.to = [((1:n)-.5)/n, (t0+2*a2*(-no:2:no)/no)/(2*pi)] + 0.05i;
end
[c h mis gx gy] = d.showimagparam(o);
title(sprintf('log10 eval errors of native quadr, and predicted, N=%d',N));
xlabel x; ylabel y; axis tight; %axis(1.2*[-1 1 -1 1]);
%print -depsc2 native_eval_tau1_nonsym_N60.eps
%print -depsc2 native_eval_tauf_nonsymbump_N60.eps
% try for helmholtz now, see if same - yes it is - see above
