%
%  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]
%  The u and h unknowns are UNSTAGGERED
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  close all;  clear;
%  parameters

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

xg = 0:dx:2*L-dx;  x = [xg 2*L];
tg = 0:dt:T;
lambda = dt/dx;
CFL = lambda

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

disp(' ')
disp(' leap-frog Shallow Water')
disp(' ')
figure(1);  clf 
plot(xg,u0,'r');  hold on;  plot(xg,h0,'g');
axis([0 2*pi -1 1.25]);
pause
%
%  Time stepping
for n=2:N
   t = n*dt;
% centered diff in space
u2 = u0 - lambda*([h1(2:M) h1(1)] - [h1(M) h1(1:M-1)]);
h2 = h0 - lambda*([u1(2:M) u1(1)] - [u1(M) u1(1:M-1)]);
% interval plotting
   if (mod(n,nplot)==0);  
      hold off
      plot(xg,u2,'r');  hold on;
      plot(xg,h2,'g');  hold off;
      axis([0 2*pi -1 1.2])
      pause;
   end
   u0 = u1;
   u1 = u2;
   h0 = h1;
   h1 = h2;
end
disp('We are at the end');
hold on
plot(xg,u2,'r'); plot(xg,h2,'g')
plot(xg,IVu,'r:');  plot(xg,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)
semilogy(k,uk)