% high freq scatt demo, and figs for CE paper. Barnett 11/1/12 - 11/3/12
% Does: Nystrom solve (file high.mat),
%       QBX eval (file high_u.mat),
%       and computes correct interpolated ueg 'exact' u vals (file high_ueg.mat)
% Each is switchable out with loading the old file instead.
% See old/fig_highfreqnystconv1.m for development code scraps.

clear; verb = 1; % decide if to produce plots for CE paper, or not

if 0 % ===================== nystrom fill+solve? =============================
  N=9000; mmax = 40; wsc = .08; side = 'e'; % rand wiggly
  o.napproxv = 3*N; rand('twister',123);  % seed
  aj = wsc*(rand(1,mmax)-0.5); bj = wsc*(rand(1,mmax)-0.5); % `rough' shape
  s = segment.smoothfourier(N,aj,bj,o);
  k = 250;
  fprintf('diam, perim, in wavelengths = %.3g, %.3g\n',2.5*k/2/pi,sum(s.w)*k/2/pi)
  de = domain([], [], s, -1);   % exterior
  o.quad = 'm';   % Kress
  %o.quad='a'; o.ord=16; % Alpert order (only needed for FMM=1)
  de.addlayerpot(s, [-1i*k 1], o); % CFIE
  s.setbc(1, 'd', []);
  p = scattering(de, []); p.setoverallwavenumber(k);
  p.setincidentwave(-pi/5);  % if just angle given, it's a plane wave
  %pt = -.48+0.95i; % test pt (near, bright region) - too optimistic
  % pt = -1.5 + 1.5i; % far in bright region - but also too optimistic
  pt = 0.92-0.42i; % test pt (near fastest Zp part, not in Omega_bad), shadow

  %o.FMM = 1; o.meth = 'iter'; % 3 options for solve
  o.FMM = 0; o.meth = 'iter';
  %o.FMM = 0; o.meth = 'direct';
  % at pt = 0.92-0.42i... reference data
  %N=9000  (107 sec): 	u(pt) = 0.923200156202017 + -0.384308681893957 i
  %N=9200  (107 sec): 	u(pt) = 0.9232001562018376 + -0.3843086818937189 i
  Ns = 9000; % can do eg 9000:200:9200
  for i=1:numel(Ns), N=Ns(i); p.updateN(N);
    tic; p.solvecoeffs(o); up = p.pointsolution(pointset(pt));
    fprintf('N=%d  (%.3g sec): \tu(pt) = %.16g + %.16g i\n', N, toc, real(up), imag(up))
  end
  fprintf('ppw = %.3g\n', N/(sum(s.w)*k/2/pi))
  %figure; tic; o = []; o.FMM=1; o.dx = 0.002; p.showthreefields(o);
  %fprintf('solution field eval in %.3g sec\n', toc), hold on; plot(s.x,'w.');
  
  % whole thing - these pts used later too
  tic; o = []; o.FMM=1; o.bb = [-2 1.3 -1.4 2];
  o.dx = 0.001; [uN gx gy di] = p.showfullfield(o);
  fprintf('soln pot on %d pts eval in %.3g sec\n', numel(find(~isnan(di))),toc)
  p.A = []; % big & don't need to keep!
  save high.mat p uN gx gy di o N side % keep the useful things
else
  fprintf('loading coeffs & uN...\n');
  load high.mat; s = p.segs(1); k = p.k; de = p.doms(1);
end

% some diagnostic plots...
%figure; semilogy(abs(fft(co))); % only drops by 1e6 - how come???

if 0 % ========================= reference eval ==========================
  ii = find(~isnan(di)); [xx yy] = meshgrid(gx, gy); zz = xx(:)+1i*yy(:);
  ueg = refeval(p, zz, ii, side);
  ueg = reshape(ueg, size(xx));
  save high_ueg.mat ueg
else
  load high_ueg.mat
end

if 1 % ====================== QBX eval =========================================
  ii = find(~isnan(di)); [xx yy] = meshgrid(gx, gy); zz = xx(:)+1i*yy(:);
co = p.co; % do this once!
 % improved *local only* correction scheme, O(N), adapted from fig_orderNstrip
 % QBX geom params: ====
delta = 2.5;  % how many "h" away center is from boundary
al = delta/N; % imag dist wrt [0,1]-scaled seg
Bdelta = 1.5; B = ceil(N/(2*Bdelta)); % how many boxes (Bdelta=boxheight/h)
als = 5.0/N; % imag strip width wrt [0,1]-scaled seg, for acc to 1e-13
fdist = 1.5; % distance factor for fine pts to include in local correction
% end params ======
fprintf('new QBX eval w/ B=%d...\n', B), tic
if side=='e', al=-al; als = -als; end % put boxes outside curve
t0 = (1:B)/B+1i*al; z0 = s.Z(t0);  % all preimages and centers
v = s.Z(s.t + 1i*als); % imag-shifted N bdry pts. Now clip off loopy ones:
vsd = abs(repmat(s.x,[1 N]) - repmat(v.',[N 1])); % find dists v to pO, O(N^2)
[vsdm, ivsd] = min(vsd,[],1); clear vsd; % get indices of nearest pts
oi.toz = s.t(ivsd);  % nearest on-pO initial guesses to help invertZ
vs = s.invertZparam(v,oi); % slow if wobbly curve
v = v(find(abs(imag(vs))>0.99*abs(als))); % snip off too-close self-int loops
fc = ~utils.inpolywrapper(zz(ii), v); % for ii pts, true if in strip
if side=='e', fc = ~fc; end % remember inside the Gamma-curve for ext
ic = ii(find(fc)); zzc = zz(ic);      % pts to correct
% QBX conv params ------  don't code conv yet
M = 26; bes = 6; % p and beta
Nf = ceil(bes*N);   % "M" in paper
nnf = 0*z0;
Rmax = 1.6*max(s.speed)*abs(als); % max dist to search for pts near each center
sf = utils.copy(s); sf.requadrature(Nf);  % fine segment
l=p.bas{1}; lpf = layerpot(sf,[-1i*k 1]); lpf.doms = de; %fine layerpot CFIE
cof = quadr.interptrig(co, Nf);
fprintf('setup done %.3g sec\n', toc);
tic; uc = lpf.evalFMM(cof,pointset(zzc));    % fine FMM eval in strip
  fprintf('fine FMM in strip = %.3g s\n',toc);
% setup a dummy local lp and seg, whose x,nx,w will be overwritten for each b:
sl = utils.copy(sf); lpl = layerpot(sl,[-1i*k 1]); lpl.doms = de; % dummy lp
tic; %profile clear; profile on;
%pthis=-.05+1.3i; bthis = find(abs(z0-pthis)==min(abs(z0-pthis))); % subset of B
  for b=1:B%bthis-15:bthis+15%loop over boxes, get indices in box, localexp-eval
  iR = find(abs(zzc-z0(b))<Rmax); % strip pt indices
  [box iz0] = l.QBXbox(real(t0(b)), [1/(2*B) als*0.4], zzc(iR));
  %  plot([box.v; box.v(1)],'-','color',[.5 .5 .5]);
  zzR = zzc(iR(iz0)); %hold on; plot(zzR, '.'); drawnow % check the set of pts
  if numel(zzR)>0
  bs = max(abs(s.Zp(t0(mod(b+[-3:3]-1,B)+1))))*abs(als); % this wei box size 
  nf = find(abs(sf.x - z0(b))<bs*fdist); % dist criterion for local fine pts
  nnf(b) = numel(nf);
  sl.x = sf.x(nf); sl.nx = sf.nx(nf); sl.w = sf.w(nf); sl.speed=sf.speed(nf);
  coR = cof(nf);   % local fine coeffs
  %ucorr = zeros(numel(zzR),1); % skip the local correction for timing tests
  ucorr = -lpl.eval(pointset(zzR)) * coR;    % subtract off direct eval
  LEo.besselcode='m'; LEo.real=0; % NB 'm' is reliable, 'r' no faster.
  LEo.rescale_rad=0; % not needed; was max(abs(zzR-z0(b)));%0.7*bs;
  LE = regfbbasis(z0(b), M, LEo); LE.doms = de; o.fast = 1; %local exp, S2Lfast?
  if l.a(1)~=0, S2L = l.a(1)*layerpot.localfromSLP(sl, LE, o); % 2 kinds of LP
     if l.a(2)~=0, S2L = S2L + l.a(2)*layerpot.localfromDLP(sl, LE, o); end
   else, S2L = l.a(2)*layerpot.localfromDLP(sl, LE); end
   coLE = S2L * coR;
   ucorr = ucorr + LE.eval(pointset(zzR)) * coLE;   % add local-exp potential
  uc(iR(iz0)) = uc(iR(iz0)) + ucorr; % correct the fine FMM pot in local box
  end
end
u = uN; % use native outside bad strip
u(ic) = p.ui(zzc) + uc; % overwrite corrected strip back into main u array
% remembering it's scattering, so need inc wave too!
%profile off; profile viewer
fprintf('local-exp correction on %d pts in strip = %.3g s\n',numel(uc),toc);
%save high_u.mat u z0 nnf M bes delta Bdelta als

%bessel diagnosing:
%b=2805; LE.besselcode='r';A =  LE.eval(pointset(zzR));
%LE.besselcode='m';B =  LE.eval(pointset(zzR));

else, fprintf('loading QBX-evaled u...\n'); load high_u.mat
end

% ============================================================================
if 1 % diagnostic plots...
figure; imagesc(gx,gy,log10(abs(u-ueg)/max(abs(ueg(:))))); axis xy equal;
colormap(goodbw); colorbar; caxis([-16 -11]);
hold on; plot(z0,'*'); plot(s.x,'g.-');
bb = o.bb; %bb = [-.08 -.01 1.26 1.32];
axis(bb)
ig=find(di(:)==1 & real(zz)>=bb(1) & real(zz)<=bb(2) & imag(zz)>=bb(3) & imag(zz)<=bb(4));
errs = u(ig)-ueg(ig);
fprintf('max rel err in bb = %.4g\n', max(abs(errs))/max(abs(ueg(:))))
fprintf('mean rel err in bb = %.4g\n', mean(abs(errs))/max(abs(ueg(:))))
% zoom to troublespot...
j=find(abs(errs)==max(abs(errs))); z = zz(ig(j));
l=.05; axis([real(z)-l real(z)+l imag(z)-l imag(z)+l]) % zoom to given z
end

if 0&verb  % paper figs --------------------------------------------------------
  figure; imagesc(gx,gy,real(u)); caxis(2*[-1 1]); colorbar;
  axis xy equal off; hold on; plot(s.x,'k-','linewidth',0.1);
  bb = [-1.18 -0.85 -.3 0.03]; %bb = [-1.2 -0.68 -.35 .15];
  plot(bb([1 2 2 1 1]), bb([3 3 4 4 3]), 'k-');
  set(gcf, 'paperposition', [0 0 6 5]);
  print -depsc2 wobbly_k250_full.eps

  figure; imagesc(gx,gy,real(u)); caxis(2*[-1 1]);
  axis xy equal off; hold on; plot(s.x,'k-','linewidth',0.1); axis(bb)
  set(gcf, 'paperposition', [0 0 6 5]); axis off
  print -depsc2 wobbly_k250_zoom.eps

  figure; imagesc(gx,gy,log10(abs(uN-ueg)/max(abs(ueg(:))))); axis xy equal;
  colormap(goodbw); colorbar; caxis([-15 0]);
  axis(bb);
  set(gcf, 'paperposition', [0 0 6 5]);
  print -depsc2 wobbly_k250_uNerr.eps
  
  figure; imagesc(gx,gy,log10(abs(u-ueg)/max(abs(ueg(:))))); axis xy equal;
  colormap(goodbw); colorbar; caxis([-15 -11]); axis(bb);
  pthis = -.95-.25i; b = find(abs(z0-pthis)==min(abs(z0-pthis)));
  bs = max(abs(s.Zp(t0(mod(b+[-3:3]-1,B)+1))))*abs(als); % this wei box size 
  nf = find(abs(sf.x - z0(b))<bs*fdist); % dist criterion for local fine pts
  o.nospill = 1;
  for b=1200:1700
    [box] = p.bas{1}.QBXbox(real(t0(b)), [1/(2*B) als*0.4], [], o);
    hold on; plot([box.v; box.v(1)],'-','color',[.5 .5 .5]);
  end
  plot(sf.x(nf), '-', 'linewidth', 2);  % local fine pts
  plot(z0(bthis), '.', 'markersize',10, 'color', [.3 .3 .3]);
  ss = newtzero(s.Zp,0.54); % Schwartz singularities of domain
  hold on; plot(s.Z(ss), '*', 'markersize',4);
  bbz = [-1.155 -1.13 -0.26 -0.235];  % zoom-zoom box
  plot(bbz([1 2 2 1 1]), bbz([3 3 4 4 3]), 'k-');
  plot([0.8*bbz(1)+0.2*bbz(2) -1.17], [bbz(4) -0.1], 'k-');
  axes('position', [0.05 0.53 0.3 0.35]);  % zoom box
  imagesc(gx,gy,log10(abs(u-ueg)/max(abs(ueg(:))))); axis xy equal tight;
  caxis([-15 -11]); axis(bbz); set(gca,'xtick',[], 'ytick',[]);
  for b=1500:1700
    [box] = p.bas{1}.QBXbox(real(t0(b)), [1/(2*B) als*0.4], [], o);
    hold on; plot([box.v; box.v(1)],'-','color',[.5 .5 .5]);
    plot(z0(b), '.', 'markersize',10, 'color', [.3 .3 .3]);
  end
  plot(s.Z(ss), '*', 'markersize',4); plot(v,'k--'); % Schwarz & Gamma_bad
  set(gcf, 'paperposition', [0 0 6 5]);
  print -depsc2 wobbly_k250_qbxerr.eps
end
