% naive sum for RPWs
clear
n = 7000;
ppw = 3.7;
real = 1; coordinner = 1;       % whether to inner loop over coords or blobs
expt = 'r';      % 'r': RPWs, 'b':test bessel J_0

ns = ((1:n)-n/2-1);   % offset integer lattice
dx = 2*pi/ppw; x = dx*ns;            % 1d position grid (same for both coords)
dk = ppw/n;    k =  dk*ns;           % 1d wavenumber grid
%N = 2*ceil(pi/dk); dt = 2*pi/N; t = ((1:N)-1/2)*dt;  % # samples around S^1
N = ceil(1.5*ceil(pi/dk)); if mod(N,2)==1, N=N+1; end  % make N even
%N = ceil(0.5*ceil(pi/dk)); if mod(N,2)==1, N=N+1; end  % make N even
dt = 2*pi/N; t = ((1:N)-1/2)*dt;  % # samples around S^1
fprintf('%d blobs\n', N);
% real case: cos coeffs then sin coeffs. complex case, c-numbers around S^1:
if expt=='b'
  if real, a = zeros(1,N); a(1:N/2)=1/(2*pi); else, a = ones(1,N) / (2*pi); end
else
  a = randn(1,N); if ~real, a = a + 1i*randn(1,N); end
end
  
f = zeros(n,n,'single'); x = single(x); a = single(a);
tic;% profile clear; profile on;
if coordinner % ........ inner coord loop (more memory needed)
[xx yy] = meshgrid(x);
if real
  for i=1:N/2, if floor(20*i/N)>floor(20*(i-1)/N),fprintf('i=%d\n',i); end
    fprintf('i=%d / %d\n', i, N/2)
  d = cos(t(i))*xx + sin(t(i))*yy;
  f = f + a(i) * cos(d) + a(i+N/2) * sin(d);
    end
    f = 2*f;
else
  for i=1:N, if floor(10*i/N)>floor(10*(i-1)/N),fprintf('i=%d\n',i); end
      f = f + a(i) * exp(1i*(cos(t(i))*xx + sin(t(i))*yy));
    end
end
else  % ....... inner blobs loop (less memory-intensive)
  kx = cos(t)'; ky = sin(t)';  % col vecs of wavevectors
  if real
  for i=1:numel(x)
    i, %if floor(10*i/n)>floor(10*(i-1)/n),fprintf('x=%g\n',x(i)); end
    for j=1:numel(x)
      d = kx(1:N/2)*x(i)+ky(1:N/2)*x(j);
      f(i,j) = a * [cos(d); sin(d)];  % real case
    end, end
    f = 2*f;
  else, error('complex coordinner=0 not implemented!');
  end
end
f = dt * f;        % quadrature approximation
toc,%profile off; profile viewer

%clear xx yy d; save f_7000.mat
load f_7000.mat

v = sum(f(:).^2)/n/n;
f = f/sqrt(v);                     % fix the std deviation to be 1
figure; if real, imagesc(x,x,abs(f)'); else, imagesc(x,x,abs(real(f))'); end
set(gca,'ydir','normal');axis equal;colormap(jet(256));caxis([0 3]);
if expt=='b', figure; [xx yy] = meshgrid(x);
  imagesc(x,x,log10(abs(besselj(0,sqrt(xx.^2+yy.^2)) - f))');
  set(gca, 'ydir', 'normal'); axis equal; caxis([-10 0]); colorbar;
  xlabel('x_1'); ylabel('x_2'); axis tight;
end

% movie of zooming out    ... ?
%fc = abs(f)>3.0;
%fcd = downsample(fc,8);

u=5.5;             % single downsampled plot
fc = abs(f)>u;
%figure; imagesc(fc');set(gca,'ydir','normal');axis equal;colormap(1-gray);
fcd = downsample(fc,8); % downsample by factor 8
figure; imagesc(fcd');set(gca,'ydir','normal');axis equal;colormap(jet(256));
%colormap(1-gray);
title('fcd'); caxis(caxis/2); m = size(fcd,1);

% ...................movie of resampled extreme value sets ..................
figure; set(gcf, 'position', [0 0 m m]); set(gca, 'position', [0 0 1 1]);
clear Mov, Mov2
cs = 0:0.04:6;
for i=1:numel(cs)
  fc = abs(f)>cs(i); fcd = downsample(fc,8);
  imagesc(fcd');axis off; set(gca,'ydir','normal');
  caxis(caxis*0.8); colormap(gray(256)); % saturate slightly
  text(0.85*m,0.95*m,sprintf('u=%.2f',cs(i)), ...
       'fontsize',36,'color',[1 0.5 0]);
  drawnow; Mov(i) = getframe(gcf);
  colormap(jet(256)); drawnow; Mov2(i) = getframe(gcf); 
end
nam = 'clev7000ppw3.7resamp8'; movie2avi(Mov,sprintf('%s.avi',nam)); % AVI (uncompr)
system(sprintf('mencoder %s.avi -o %s.mpg -ovc lavc', nam, nam)); % make mpeg
nam = 'clev7000ppw3.7resamp8c'; movie2avi(Mov2,sprintf('%s.avi',nam)); % AVI (uncompr)
system(sprintf('mencoder %s.avi -o %s.mpg -ovc lavc', nam, nam)); % make mpeg