% plain old test of Nystrom convergence rate, multiple types of error,
% compare to the mystery of decay in fourier coeffs of soln. Int/ext BVPs.
% Laplace/Helm, Dir/Neu. (Except not the ext Neu case).
% The "exact" values at each N are computed via the 3N soln, inefficient!
%
% Discussion: when zsing chosen close to bdry, tau has a sing at same place.
% Thus it controls convergence of Nystrom, via PTR convergence.
%
% Issue is Nystrom interpolation vs trig-poly interp: turns out Nystrom
% converges twice as fast as trig poly, for k=0!
% But for k>0 Nyst is same rate due to product quadrature, so Nyst interp not
% useful to be implemented here.
% For Neu int, implemented Atkinson removal of nonuniqueness & fitting of offset
% For ext scatt BVP, uses plane wave is set zsing=[];  Ext Neu not coded up.
%
% Barnett 10/3/12, added Neu BCs & fixed Nystrom interp for Neu case, 10/27/12
clear; bc = 'D'; side = 'e'; % 'D' or 'N';  'i' interior BVP, 'e' exterior BVP
s = segment.smoothnonsym(100, 0.3, 0.3, 3); sf = utils.copy(s);
zsing = []; % data for expt=f: singularity loc, if empty, use entire func
if side=='i', zsing = 1+.5i; else, zsing = 0.6+.2i; end
if side=='i', d = domain(s, 1); else d = domain([],[], s, -1); end
k = 30; d.k = k;  % choose Laplace or Helmholtz % k=10
o = [];
if side=='i'
  if bc=='D', d.addlayerpot(s, 'D', o); else, d.addlayerpot(s, 'S', o); end
else, d.addlayerpot(s, [-1i*k 1], o); end % CFIE
l = d.bas{1};  % DLP or SLP
alsing = 2*pi*min(abs(imag(s.invertZparam(zsing)))); % alpha sing wrt (0,2pi)
if k==0    % Laplace data
  if ~isempty(zsing), f = @(z) log(abs(z-zsing));
    fx = @(z) real(z-zsing)./abs(z-zsing).^2; fy = @(z) imag(z-zsing)./abs(z-zsing).^2;
  else, f = @(z) cos(real(z)).*exp(+imag(z)); % needs generic Taylor exp
    fx = @(z) -sin(real(z)).*exp(+imag(z)); fy = f;
  end
else      % Helmholtz data
  if ~isempty(zsing), f = @(z) 3*besselh(0,k*abs(z-zsing));  % u = pt src
    f = @(z) 3*besselh(0,k*abs(z-zsing)); % u = radiating pt src
    fx=@(z) -3*k*real(z-zsing)./abs(z-zsing).*besselh(1,k*abs(z-zsing));
    fy=@(z) -3*k*imag(z-zsing)./abs(z-zsing).*besselh(1,k*abs(z-zsing));
  else, ang=exp(1i*pi/3); f = @(z) exp(1i*k*real(ang'*z));  % u = plane wave
    fx = @(z) 1i*k*real(ang)*f(z); fy = @(z) 1i*k*imag(ang)*f(z);
  end
end
d.bas{1}.quad='m'; % Kress for self-eval (if kappa known & N<1400)
%d.bas{1}.quad='a'; d.bas{1}.ord = 16; % Alpert for self-eval (if no kappa)
p = bvp(d);
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
if side=='i', s.setbc(-1, bc, [], data); else s.setbc(1, bc, [], data); end %BVP
Ns = 60:20:360; % 340
noexf = isempty(zsing) & side=='e'; % no exact f known, so use converged val
for i=1:numel(Ns), N = Ns(i);  % ............ convergence
  p.updateN(N); p.fillquadwei; p.sqrtwei=0*p.sqrtwei+1; % NB unit sqrtwei!
  p.fillbcmatrix;
  atk = (k==0 & bc=='N'); % use Atkinson book handling nonunique int Neu
  if atk, p.A(:,1) = p.A(:,1) + 1; end % Atkinson  p* = s.x(1)
  p.fillrighthandside; p.linsolve; co=p.co; w = p.segs.w'; % one segment only
  fco = abs(fft(co)); ef(i) = fco(N/2+1)/max(fco);
  rad = .3; if side=='e', rad = 3; end
  zp = rad*exp(1i*2*pi*(1:10)'/10);  % check at some pts deep into domain
  up = p.pointsolution(pointset(zp));
  off = 0; if atk, off = up(1)-f(zp(1)); end % Neu soln offset
  p.updateN(3*N); p.fillquadwei; p.sqrtwei=0*p.sqrtwei+1;
  p.fillbcmatrix; if atk, p.A(:,2) = p.A(:,2) + 1; end % Atkinson p*=sf.x(2)
  p.fillrighthandside; p.linsolve; upe = p.pointsolution(pointset(zp));
  if ~noexf, ei(i) = max(abs(up-off-f(zp))); % pointwise err rel to exact
    else ei(i) = max(abs(up-upe)); end % ptwise err rel to 3N values
  if atk, p.A(:,2) = p.A(:,2) - 1; end %undo, but doesn't matter either way
  coa = p.co(2:3:end); % "exact" tau @ nodes
  et(i) = norm(coa-co)/norm(coa);
  eti(i) = norm(p.co-quadr.interptrig(co,3*N))/norm(p.co); % est trig interp err
  if k==0, Af = p.A; if bc=='D', JR=-1; else JR=1; end % diag term to kill
    Af(diagind(Af)) = diag(Af) - JR/2; % axe jump term fine seg
    Af=Af(:,2:3:end); cof = 2*JR*(p.rhs-3*Af*co); % cols of Nyst mat, Nyst interp - only works for Laplace, not k>0
  % (note factor of 3 is to change to weights on N rather than on 3N)
    etn(i)=norm(cof-p.co)/norm(p.co); else, etn(i)=nan; end %Nyst interp est err
  fprintf('N=%d errs (pts, tau, taui, taun) %.3g %.3g %.3g %.3g  fou_decay %.3g\n',N,ei(i),et(i),eti(i),etn(i),ef(i))
end                           % ..............
figure; semilogy(Ns,[ei;et;eti;etn;ef], '+-'); v=axis; v(3)=1e-16; axis(v);
legend('int pts err','tau_j L2 err','tau interp L2 err', 'tau Nyst-int L2','tau fou decay');
if ~isempty(zsing), hold on; plot(Ns,1e-2*exp(-[Ns;Ns/2]*alsing),'--');end


%figure; semilogy(abs(fft(co)));hold on;title('\tau fourier decay');
%if ~isempty(zsing), m=1:N/2; plot(m,exp(-alsing*m),'r-'); end, drawnow;
