% cbm, 20060303

% warning: might (ha!) contain some bugs

% if anyone can see how to get rid of the fftshift inside the surf
% plot, I'd love to know:  FIXED thanks to David Alexander.

Nx=100;  Ny=100;
Lx=1;  Ly=1;

dx=Lx/Nx;  dy=Ly/Ny;
x=[-Lx/2:dx:Lx/2-dx];
y=[-Lx/2:dy:Ly/2-dy];
[xg,yg]=meshgrid(x,y);

dk=2*pi/Lx;  dj=2*pi/Ly;
kk=[-Nx/2:Nx/2-1]*dk;
ll=[-Ny/2:Ny/2-1]*dj;
[kg,lg]=meshgrid(kk,ll);

% initial conditions: a Gaussian blob
% there must be a bug here because this has height one, but the
% scale of the solutions is more like 5000.  There must be a 1/dx
% or something hanging around
q=exp(-(xg.^2 + yg.^2)/(1/1000));

surf(q);
pause(.1);

qh=fftshift(fft2(fftshift(q)));
%qh=fft2(q);

surf(real(qh));
title('real');
pause(.1);
surf(imag(qh));
title('imag');
pause(.1);

qh(:,1) = 0;
qh(1,:) = 0;
a = real(qh(:));

% it would be too slow to use all the modes, so we make an index
% for those modes larger than something
Indices = find(abs(a) > 0.25);
% You could try this if your computer is fast
%Indices = 1:N

disp('how many modes?:');
size(Indices)
disp('percentage:');
length(Indices)/length(a)

ks = sqrt(kg(:).^2+lg(:).^2);
N = length(ks);

kset = kg(:);
lset = lg(:);

g = 9.81;
h = 1;

dt = .01;
Nt = 150;
for jt = 1:Nt;
  t = dt*jt
  eta = zeros(size(xg));
  %for j = 1:N
  % iterate of the set of Indices
  for jj = 1:length(Indices)
    j = Indices(jj);
    omega = sqrt(g*ks(j)*tanh(ks(j)*h));
    k = kset(j);
    l = lset(j);
    eta = eta + a(j)*cos(k*xg + l*yg - omega*t);
  end

  figure(1); clf;
  H = surf(x,y,(eta));
  %H = surf(x,y,fftshift(eta));
  title(['t = ' num2str(t)]);
  xlabel('x'); ylabel('y');
  axis([-.5 .5 -.5 .5 -5000 5000]);
  shading interp;
  fname = sprintf('frame%.3d', jt);
  print('-dpng', '-r75', fname);
  pause(0);
end
