% shape and integrand on pO for FACM2012, LE coeff

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

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));
z0 = 0.45-0.3i;
R = 0.5; M = 18; Nf = 4*N;
iz0 = abs(zz-z0) < R;
LEo.besselcode='m'; LEo.real=0; LEo.rescale_rad=R;
LE = regfbbasis(z0, M, LEo); LE.doms = d;  % local exp
sf = utils.copy(s); sf.requadrature(Nf); cof = quadr.interptrig(co, Nf);
S2L = layerpot.localfromDLP(sf, LE); coLE = S2L * cof; % LE coeffs
u(iz0) =  LE.eval(pointset(zz(iz0))) * coLE;   % add local-exp potential
u(~di) = nan;


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));
m = 10; plot(sf.t*2*pi, real(S2L(m,:))./sf.w, '-'); xlabel('s');
axis tight;

print -depsc2 talk_native_eval_tau1_nonsym_N60_LE.eps

