%
%  lect25.m -- leapfrog one-way wave
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

N =    64; L =   pi;  dx = 2*L/N;
M = 3/8*N; T =  pi/2;  dt =   T/M;

%  set grids & coefficients

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

lambda = dt/dx;

c  = 1.0 + 0.0*sin(xg);
f  = 0.0;

%  row vector of IVs

IV0  = [exp(-4*(xg-pi).^2)];
IV1  = [exp(-4*(xg-1*c*dt-pi).^2)];
%IV1  = [exp(-4*(xg-2*c*dt-pi).^2)];

%  PDE solve with operation count

disp(' ')
disp(' fully-discrete one-way wave')
disp(' ')
figure(1); clf
flops(0);

u = IV0;  
plot(xg,real(u),'r');  hold on;  plot(xg,real(u),'kx');
axis([0 2*pi -0.25 1.25]);

% exact next step
   uold = u;
   u    = IV1;

for k=2:M
   unew  = uold - lambda*c.*([u(2:end) u(1)] - [u(end) u(1:end-1)]) + f;

   uold = u;  u = unew;

% interval plotting
%   if (mod(k,M/8)==0);  plot(xg,real(u),'k');  end
end

plot(xg,real(u),'b');  hold on;  plot(xg,real(u),'kx');

opcount = flops;  
disp(' ');  disp([' flops = ' num2str(opcount)]);  disp(' ')
xlabel('\bf x-axis');  ylabel('\bf u-axis')
title('\bf one-way wave (c=1,f=0)')