% script file lin_adv.m
% linear advection equation with
%        upward differences
%        downward differences
%        Lax-Wendroff scheme
%        Lax-Friedrich scheme
% and periodic initial data
%
% last changes: 02.06.03 by Vita
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clear all; close all;

EXAMPLE = 'B'; % 'A' or 'B' for different initial data
METHOD = 'downward'; %'upward','downward','LW' or 'LF'
ANIMATE = 1; % 1 or 0
             % 1: shows each time step
	     % 0: only 4 moments plotted

Tmax = 10; 


a = 1; % advection speed
n = 101; % #(grid points) in x-direction

CFL = 1; 


% initialize 
t = 0;
u0 = zeros(n,1);
if isequal(EXAMPLE,'A')
  ax = 0; bx = 1;
  x = linspace(ax,bx,n);
  
  left = find(x<0.5);
  right = find(x>=0.5);
  u0(left) = x(left);
  u0(right) = 1-x(right);
elseif isequal(EXAMPLE,'B')
  ax = -5; bx = 2;
  x = linspace(ax,bx,n);
  
  left = find(x<0);
  right = find(x>=0);
  u0(left) = 1;
  u0(right) = -1;
else
  error('Only A and B examples implemented!')
end

% spatial step
h = x(2)-x(1);

% choose k (time step) according to CFL
% CFL = a*lambda = a*k/h
% => k = CFL*h/a
k = CFL*h/abs(a);
lambda = k/h;

u = u0;
u_new = zeros(n,1);
figure(1); 

% PLOTTING TIMES

if ANIMATE==0
  time_steps = floor(Tmax/k); % how many steps will be done
  plot_time_step = floor(time_steps/4);
  plot_times_ind = [0:plot_time_step:time_steps]; % indices
  plot_times = k*plot_times_ind; % values
  plot_ctr = 0;
  
  subplot(4,1,1)
  title(['Method: ',METHOD,' Example: ',EXAMPLE,' CFL = ',num2str(CFL)])
  hold on
end

while t<Tmax
  t = t+k;
  
  switch METHOD
   case 'upward'
    u_new = u - a*lambda*([u(2:n);u(1)] - u);
   case 'downward'
    u_new = u - a*lambda*(u - [u(n);u(1:n-1)]);
   case 'LW'
    u_new = u - a*lambda/2*([u(2:n);u(1)]-[u(n);u(1:n-1)]) + ...
	    a^2*lambda^2/2*([u(2:n);u(1)]-2*u+[u(n);u(1:n-1)]);
   case 'LF'
    u_new = ([u(2:n);u(1)]+[u(n);u(1:n-1)])/2 - ...
	    a*lambda/2*([u(2:n);u(1)]-[u(n);u(1:n-1)]);
  otherwise
    error('only upward, downward, LW and LF implemented!')
  end
  
  u = u_new;

  if ANIMATE==1
    plot(x,u,'b.-')
    drawnow
  elseif ANIMATE==0
    
    % VISUALIZE
    if min(abs(t-plot_times))<=1e-12
      plot_ctr = plot_ctr+1;
      subplot(4,1,plot_ctr);
      plot(x,u,'b.-',x,u0,'r')
      legend(sprintf('t=%2.2f',t),0)
      drawnow;
    end
  else
    warning('no visualization!')
  end
  
end



  
  
  
  