% plain old test of Nystrom convergence rate, multiple types of error,
% compare to the mystery of decay in fourier coeffs of soln. Interior BVPs only.
% Laplace/Helm, Dir/Neu.
%
% 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
%
% Barnett 10/3/12, added Neu BCs & fixed Nystrom interp for Neu case, 10/27/12
clear; bc = 'D'; % 'D' or 'N'
s = segment.smoothnonsym(100, 0.3, 0.3, 3); sf = utils.copy(s);
zsing = [];%1+.5i; % data for expt=f: singularity loc, if empty, use entire func
d = domain(s, 1); k = 0; d.k = k;  % choose Laplace or Helmholtz % k=10
o = [];
if bc=='D', d.addlayerpot(s, 'D', o); else, d.addlayerpot(s, 'S', o); end
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
    %fx = @(z) 3*besselh(1,k*abs(z-zsing));
    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
  % ****** fix for besselh
  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
s.setbc(-1, bc, [], data); % BVP
Ns = 60:20:340; %330;
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; if bc=='N', 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);
  zp = .3*exp(1i*2*pi*(1:10)'/10);  % check at some deep interior pts
  off = 0; if bc=='N', off = p.pointsolution(pointset(zp(1)))-f(zp(1)); end
  ei(i) = max(abs(p.pointsolution(pointset(zp))-off-f(zp))); % Neu allow offset
  p.updateN(3*N); p.fillquadwei; p.sqrtwei=0*p.sqrtwei+1;
  p.fillbcmatrix; if bc=='N', p.A(:,2) = p.A(:,2) + 1; end % Atkinson p*=sf.x(2)
  p.fillrighthandside; p.linsolve;
  if bc=='N', 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;
