% Check pointwise evaluation error of native quadrature scheme in BIE
% Barnett 1/29/12

clear all classes; expt='1'; N = 50;
%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) });
% ........

d = domain(s, 1); d.k = 0; %d.plot;
o = []; d.addlayerpot(s, 'D', o);
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'
  f = @(z) -1 + 0*z;               % known solution -1 inside
  co = ones(N,1);                    % unit density
elseif expt=='f'
  p = bvp(d);  % now solve a BVP to get a more interesting density
  d.bas{1}.quad='a'; d.bas{1}.ord = 16; % make alpert for self-eval
  f = @(z) real(z).^2 - imag(z).^2;    % u(x,y) = x^2 - y^2   exact soln
  s.setbc(-1, 'D', [], @(t) f(s.Z(t)));
  p.solvecoeffs; %opts.comparefunc = f; p.showsolution(opts);
  co = p.co;
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));

figure; errs = 10.^(-8: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(1./errs)/N; %[0:.05:.3]; % use complex map imag shifts
s.requadrature(1000);
for i=1:numel(taus), t = [s.t; s.t(1)] + 1i*taus(i)/(2*pi);
  plot(s.Z(t), 'b-', 'linewidth', 2); hold on; end
axis equal; plot([s.x; s.x(1)], 'k+-');
title(sprintf('log10 eval errors of native quadr, and predicted, N=%d',N));
xlabel x; ylabel y; axis(1.2*[-1 1 -1 1]);
%print -depsc2 native_eval_nonsym_N50.eps

% try for helmholtz now, see if same!
