% numerical rank of J_l(xy) kernel on [0,1] x [0,2pi.W]. Barnett 5/22/19
clear; close all

W = 20;   % overall dimless freq param = K/2pi
m = 2e2; x = (0.5:m-0.5)/m; dx = 1/m;   % discr grid
y = 2*pi*W*x;
l = 5;            % order
ker = besselj(l,y'*x);   % kernel func eval on pts
figure(1); imagesc(x,y,ker); axis xy; xlabel x; ylabel y; colorbar
title(sprintf('J_%d(xy) over [0,1]x[0,K],  K=%.3g',l,2*pi*W))
discrop = sqrt(2*pi*W)*dx*ker;      % prefac is for quadrature in integral op
fprintf('2-norm of operator = %.3g\n',norm(discrop))
% explore rank vs log(1/eps) ...
for digits=2:2:10
  eps = 10.^-digits;
  fprintf('eps=%.3g\t eps-rank=%d\n',eps,rank(discrop,eps))
end

pause;      % l animation, fixed eps ...
eps = 1e-2;
figure(2);
for l=0:8*W         % orders
  ker = besselj(l,y'*x);
  imagesc(x,y,ker); axis xy; xlabel x; ylabel y; caxis(0.5*[-1 1]); colorbar;
  discrop = sqrt(2*pi*W)*dx*ker;      % prefac is for quadrature in integral op
  s=sprintf('K=%.3g\t l=%d\t eps-rank=%d',2*pi*W,l,rank(discrop,eps));
  title(s); disp(s); drawnow;
end
