% Fig nysttauconv.
% plain old test of Nystrom convergence rate, multiple types of error,
% compare to the mystery of decay in fourier coeffs of soln.
%
% 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.
%
% Barnett 10/17/12
clear; figs = 1;
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 = []; d.addlayerpot(s, 'D', o); l = d.bas{1};  % DLP
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));
  else, f = @(z) cos(real(z)).*exp(+imag(z)); %was f = @(z) real(z).*imag(z);
  end
else      % Helmholtz data
  if ~isempty(zsing), f = @(z) 3*besselh(0,k*abs(z-zsing));  % u = pt src
  else, f = @(z) exp(1i*k*real(exp(1i*pi/3)'*z)); end        % u = plane wave
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); s.setbc(-1, 'D', [], @(t) f(s.Z(t))); % BVP
Ns = 60:10:360;
for i=1:numel(Ns), N = Ns(i);  % ............ convergence
  p.updateN(N); p.fillquadwei; p.sqrtwei=0*p.sqrtwei+1; % NB unit sqrtwei!
  p.solvecoeffs; 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
  ei(i) = max(abs(p.pointsolution(pointset(zp))-f(zp)));
  p.updateN(3*N); p.fillquadwei; p.sqrtwei=0*p.sqrtwei+1;
  p.solvecoeffs; coa = p.co(2:3:end); % "exact" tau soln
  et(i) = norm(coa-co)/norm(coa);
  eti(i) = norm(p.co-quadr.interptrig(co,3*N))/norm(p.co); % est trig interp err
  Af = p.A; Af(diagind(Af)) = diag(Af) + 1/2; % remove the jump term on fine seg
  Af=Af(:,2:3:end); cof = -2*(p.rhs-3*Af*co); % cols of Nyst mat, Nyst interp
  % (note factor of 3 is to change to weights on N rather than on 3N)
  etn(i)=norm(cof-p.co)/norm(p.co);  % Nyst interpolant 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

if figs
  figure;
  semilogy(Ns,ei,'+-b',Ns,et,'o-g',Ns,eti,'x-r',Ns,etn,'*-c',Ns,ef,'d-m');
  semilogy(Ns,ei,'+-b',Ns,et,'o-g',Ns,etn,'x-r',Ns,eti,'*-c',Ns,ef,'d-m');
  axis([60 340 1e-16 1e-3]);
  legend('interior pt error','\tau_j l^2 error','\tau Nystrom interp L^2', '\tau trig poly interp L^2','Fourier coeff decay');
  hold on; plot(Ns,1e-1*exp(-[Ns;Ns/2]*alsing),'k--');
  xlabel('N'); ylabel('error');text(120,1e-4,'(a)','fontsize',14);
  set(gcf,'paperposition',[0 0 5 3]);
  print -depsc2 rfn_k0_zsing_nysttauconv.eps
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;
