%Demo for water wave
g=9.8;H=0.1; 
L=1000; %shallow water or long wave
%L=10; % intemediate depth
%L=0.1; % deep water

nx=10000;
x=((1:nx)-nx/2-1)/nx*L;
eta0=exp(-(x/0.01/L).^2);
%eta0=zeros(1,nx);eta0(find(abs(x/L)<0.1))=1;%a top-hat initial condition

u0=0*x;
etaf=fft(eta0);

%wavenumber
wn=[(0:nx/2) (-nx/2+1:1:-1)]/L*2*pi;
%frequency from the dispersion
omega=sqrt(g*abs(wn).*tanh(abs(wn)*H));

%plot(abs(wn)*H,omega./abs(wn)/sqrt(g*H))

nt=30;
dt=L/4/nt;

for it=1:nt 
  eta=real(ifft(etaf.*(exp(-1i*omega*dt*it)+exp(1i*omega*dt*it))/2));
  %add up the two eigenmodes
  tmp=eta;tmp(:)=0;
  tmp=[0 g*wn(2:end)./omega(2:end).*etaf(2:end).*exp(-1i*omega(2:end)*dt*it)-g*wn(2:end)./omega(2:end).*etaf(2:end).*exp(1i*omega(2:end)*dt*it)]/2;  
  u=real(ifft(tmp));
  
  set(gca,'NextPlot','replacechildren')
  subplot(2,1,1)
  hold off
  plot(x/L,eta0,'r')
  hold on
  plot(x/L,eta)
  xlabel('x/L');ylabel('eta');
  axis([-0.5 0.5 0 1])
  subplot(2,1,2)
  plot(x/L,u)
  xlabel('x/L');ylabel('u');

  h=gcf;
  M(:,it)=getframe(h);
end
%movie(h,M,1)


