%
%  lect17.m -- semi-discrete one-way wave
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

movie = 1;

%  parameters

global c dx

c = 1.0;  
M =   16;  T = 2*pi;  dt = T/M;
N =   80;  L =   pi;  dx = 2*L/N;

%  set grids & IVs

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

IVs  = [sin(2*xg + cos(xg))]';
%IVs  = [sign(sin(xg))]';
Xode = 'X_1way';  method = 'ode45';

Tacc = 1;

options = odeset('reltol',Tacc*1e-3,'abstol',Tacc*1e-6,'stats','on');

disp(' ')
disp([method ' for ' Xode])
disp(' ')

%  ODE solve with operation count

flops(0);
[t y] = feval(method,Xode,tg,IVs,options);
disp(' ')
opcount = flops;  disp(['flops = ' num2str(opcount)])

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  solution plot

if (movie~=1) then

figure(1);clf

yMin = min(min(y));  yMax = max(max(y));  

subplot(2,1,1);  title('\bf one-way wave:  upwind difference')
%subplot(2,1,2);  title('\bf one-way wave:  centred difference')

xlabel(['\bf x-axis:  N = ' num2str(N)]);  ylabel('\bf u-axis');
axis([0 2*L-dx yMin yMax])

for k=1:size(t,1)
%  plot(xg,y(k,:),'kx')
  hold on
  plot(xg,y(k,:),'b--')
  pause(0)
end

plot(xg,y(1,:),'r')
plot(xg,y(end,:),'r')

error = sqrt(sum((y(1,:)-y(end,:)).^2)*dx)

else 

%
%  movie version
%

figure(1);clf
yMin = min(min(y));  yMax = max(max(y));  

subplot(2,1,1);

for k=1:size(t,1)
  plot(xg,y(1,:),'r')
  hold on
%  plot(xg,y(k,:),'kx')
  plot(xg,y(k,:),'b--')
  hold off

  title('\bf one-way wave:  the movie')
  xlabel(['\bf x-axis:  N = ' num2str(N)]);  ylabel('\bf u-axis');
  text(pi/10,-0.5,['time = ' num2str(t(k))]); 
  axis([0 2*L-dx yMin yMax])
  pause(0.2)
end

end