% Check pointwise evaluation error of native quadrature scheme in BIE
% Laplace & Helmholtz. Recall Neu int BVP for Lapl is ill-posed.
% Barnett 1/29/12, 3/9/12. SLP Neu BVP added. 9/27/12 subfigures

clear; sf = 'b';  % subfigure choose 'a','b', etc...
expt='1'; bc = 'D'; N = 60; % choose expt, BC (D or N), base # bdry pts
k = 0;  % choose Laplace or Helmholtz
s = segment.smoothnonsym(N, 0.3, 0.3, 3);  % a b d
if sf=='c' 
  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
end
if sf=='b', expt='f'; bc = 'D';
  rs = 1; %refinement factor for BIE soln (expt='f' only)
end
if sf=='d', expt='g'; k = 2; end  % GRF for Helmholtz

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

if k==0, f = @(z) real(z).*imag(z);   % set up an exact known PDE soln
  % note: better since not zero along the ray 3*pi/4 that bump is on!
  fx = @(z) imag(z); fy = @(z) real(z);  % partials
  % was: real(z).^2 - imag(z).^2;  % u(x,y) = x^2 - y^2 exact, Lap
else, %kvec = k*exp(1i*pi/3); f = @(z) exp(1i*real(kvec'*z)); % u=plane wave
  %fx = @(x) 1i*real(kvec)*f(x); fy = @(x) 1i*imag(kvec)*f(x); 
  zsing = 1+1i; R = @(z) abs(z - zsing); ampl = 3; % point source
  f = @(z) ampl*bessely(0,k*R(z));
  df = @(z) -ampl*k*(z-zsing)./R(z).*bessely(1, k*R(z)); % grad u as C# (f re)
  fx = @(z) real(df(z)); fy = @(z) imag(df(z)); % hack only for f real
end

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'
  if rs>1, s.requadrature(3*N); end % 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)
  maxtau = max(abs(f(s.x)));  % obviously this *isn't* tau max in near strip!
  if bc=='D', data = @(t) f(s.Z(t)); % choose Dir or Neu data of f
  else, data = @(t) fx(s.Z(t)).*real(s.Zn(t)) + fy(s.Z(t)).*imag(s.Zn(t)); end
  s.setbc(-1, bc, [], data);
  p.solvecoeffs; %opts.comparefunc = f; p.showsolution(opts);
  if rs>1, co = p.co(2:3:end-1); else, co = p.co; end % extract the coincident nodes (3N->N downsample)
  s.requadrature(N); % for plotting later
elseif expt=='g'                    % GRF
  co = [fx(s.x).*real(s.nx) + fy(s.x).*imag(s.nx); -f(s.x)]; % GRF data
  maxtau = mean(abs(co));  % not true
end

fz = nan*xx; fz(ii) = f(zz(ii));
if expt~='f', 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, % usual color image:
figure; imagesc(xx(1,:),yy(:,1),log10(abs(u-fz))); colorbar; caxis([-16 0]);
axis equal xy tight;
end

% contour comparison:
figure; errs = 10.^(-14:2:0); %10.^(-16:2:0);
contourf(xx,yy,log10(abs(u-fz)), log10(errs)); caxis([-16 0]);
%colormap(goodbw); %set(gca,'fontsize',14);
c = goodbw(-9); colormap(c(1:end-1,:)); % needed to get discrete levels in cbar
colorbar; axis equal; hold on; plot([s.x; s.x(1)], 'k.-');
%if sf=='b', errs = 10.^(-8:2:0); end  % if need different prediction contours
taus = log(maxtau./errs)/N; %[0:.05:.3]; % use complex map imag shifts
oi = []; oi.outside = (sf=='d');  % 1 for contours out&in, 0 for just inside
oi.levels = taus/(2*pi);if exist('t0')==1 % hack get bump to plot imag dist ok
  n=10; no=10; oi.to = [((1:n)-.5)/n, (t0+2*a2*(-no:2:no)/no)/(2*pi)] + 0.04i; %0.05i;
end
oi.bb = bb; [c h mis gx gy] = d.showimagparam(oi);
axis xy tight off; %xlabel x; ylabel y;
if exist('zsing')==1, hold on; plot(zsing,'k.','markersize',20), end

if 0   % EPS OUTPUT
  set(gcf,'paperposition',[0 0 6 5]);
if sf=='a',
  print -depsc2 native_eval_tau1_nonsym_N60.eps
elseif sf=='b'
  print -depsc2 native_eval_tauf_nonsym_N60.eps
elseif sf=='c'
  print -depsc2 native_eval_tau1_nonsymbump_N60.eps
elseif sf=='d'
  print -depsc2 native_eval_grf_nonsym_N60_inout.eps
end
end

