% shape and integrand on pO for FACM2012

clear; expt='1'; N = 60;   % base # bdry pts
s = segment.smoothnonsym(N, 0.3, 0.3, 3);
z0 = 0.6-0.5i;

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

f = @(z) -1 + 0*z;               % known solution -1 inside
maxtau = 1;
co = ones(N,1);                    % unit density

fz = nan*xx; fz(ii) = f(zz(ii)); if expt=='1', fz(find(~di)) = 0; end %tau=1 out
fz = reshape(fz, size(xx)); % fill fz, u...
u = nan*xx; u(ii) = A*co; u = reshape(u, size(xx));

if 0
figure; axes('position',[0.1 0.3 0.6 0.7]);
errs = 10.^(-16:2:0);
contourf(xx,yy,log10(abs(u-fz)), log10(errs)); caxis([-16 0]); colorbar;
axis equal; set(gca,'ydir','normal'); axis tight; axis off;
hold on; plot([s.x; s.x(1)], 'k.-');
plot(z0, 'k.','markersize', 20); text(real(z0),imag(z0)+.15, 'z','fontsize',26);

axes('position',[0.1 0.07 0.6 0.2]);
s.requadrature(N); I = d.bas{1}.eval(pointset(z0));
plot(s.t*2*pi, I./s.w, '.', 'markersize',10); hold on; % integrand
s.requadrature(1e3); % for plotting smooth curve
If = d.bas{1}.eval(pointset(z0));
plot(s.t*2*pi, If./s.w, '-'); axis([0 2*pi -1 .4]); xlabel('s');

print -depsc2 talk_native_eval_tau1_nonsym_N60.eps
end

figure; errs = 10.^(-16:2:0); %10.^(-16:2:0);
contourf(xx,yy,log10(abs(u-fz)), log10(errs)); caxis([-16 0]);
%set(gca,'fontsize',14);
%colormap(1-gray(256));   % for paper
colorbar; axis equal; hold on; plot([s.x; s.x(1)], 'k.-');
errs = 10.^(-14:2:0);   % if need different prediction conts: for (a)
%errs = 10.^(-8:2:0);   % if need different prediction contours
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
%o.outside = 1;  % switch this: 1 for contours out&in, 0 for just inside
o.bb = bb; [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 off; %axis(bb); %axis(1.2*[-1 1 -1 1]);


print -depsc2 talk_native_eval_tau1_nonsym_N60_pred.eps

