%
%  leap-frog, central differences for the linearized shallow water
%  equations (scaled)
%     u_t +  h_x = 0, periodic on [0,2L]
%     h_t +  u_x = 0, periodic on [0,2L]
%
%  u and h unknowns are STAGGERED
%
%  G = g H/U^2 from the original scaled problem and in most cases U = 1
%  except we can set U=0 to turn this term off
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  close all;  clear;
%  parameters

global dx dt
%
        disp(' ');
        disp('Leap Frog Shallow Water, Staggered')
        M =  64;  L = pi;  dx = 2*L/M;
        N =  128;  T = 2*pi;  dt = T/N; % MC
%        N =  14.25*128;  T = 14.25*2*pi;  dt = T/N; % BEN	
        nplot = 8;
%        nplot = 1;	
k = [-M/2:M/2-1];
        
%  set grids & coefficients

xg = 0:dx:2*L-dx;  xh = xg + 0.5*dx;
tg = 0:dt:T;
lambda = dt/dx;
CFL = lambda

u0  = [exp(-4*(xg-L).^2)];
u0k = fftshift(fft(u0));  k = [-M/2:M/2-1];  u0k = abs(uk);
figure (2);
semilogy(k,u0k,'rx');  hold on
h0  = zeros(size(xg));
IVu = u0;  IVh = h0;
%
%  Forward Euler, central differences first time step
u1 = u0 - lambda*(h0 - [h0(M) h0(1:M-1)]);
h1 = h0 - lambda*([u0(2:M) u0(1)] - u0);
%
%  PDE solve with operation count

disp(' ')
disp(' leap-frog Shallow Water')
disp(' ')
figure(1);  clf 
plot(xg,u0,'r');  hold on;  plot(xh,h0,'g');
axis([0 2*pi -1 1.25]);

%
%  Time stepping
for n=2:N
   t = n*dt;
% centered diff in space
u2 = u0 - 2*lambda*(h1 - [h1(M) h1(1:M-1)]);
h2 = h0 - 2*lambda*([u1(2:M) u1(1)] - u1);
% interval plotting
   if (mod(n,nplot)==0);  
      hold off
      plot(xg,u2,'r');  hold on;
      plot(xh,h2,'g');  hold off;
      axis([0 2*pi -1 1.2])
   end
   u0 = u1;
   u1 = u2;
   h0 = h1;
   h1 = h2;
end

disp('We are at the end');
hold on
plot(xg,u2,'r'); plot(xh,h2,'g')
plot(xg,IVu,'r:');  plot(xh,IVh,'g:');

xlabel('\bf x-axis');  ylabel('\bf u-axis')
title('\bf one-way wave (c=1,f=0)')

uk = fftshift(fft(u2));  uk = abs(uk);
figure(2)
hold on
semilogy(k,uk,'b+')