function [x, t, u] = advectloop(tmax, level, lambda, al, x0l, deltal)
% advcectloop Solves the advection equation using an O(h^2) FDA.
%
% Generates finite-difference solution of advection equation using O(h^2) 
% difference techniques (except for right boundary, where O(h) 
% method is used) as discussed in class.
%
% This implementations uses a loop over the spatial index rather
% that whole-array operations.   See advect.m for the whole-array
% version.
%
% Input arguments
%
%      tmax:   (real scalar) Final solution time.
%      level:  (integer scalar) Discretization level.
%      lambda: (real scalar) Courant number: delt / delx.
%      al, x01, deltal: (real scalars) Parameters defining initial 
%        gaussian profile.
%
% Output arguments
%
%      t:      (real vector) Vector of length nt (where nt is defined
%              implicitly via level, tmax and lambda) containing 
%              discrete times (time mesh).
%      x:      (real vector) Vector of length nx = 2^level + 1
%              containing discrete space coordinates (space mesh).
%              angular displacement at discrete times t(n).
%      u:      (real 2D array) 2D Array of size nt x nx containing
%              full solution of finite difference equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % Enables/disables local tracing and frequency of tracing in 
   % main evolution loop
   trace = 0;
   tracefreq = 50;

   if trace
      fprintf('advectloop: Argument dump follows\n');
      tmax, level, lambda, al, x0l, deltal, ar, x0r, deltar
   end

   % Define 0-length defaults for output vectors/matrices
   x = zeros(0,0);
   t = zeros(0,0);
   u = zeros(0,0);

   % Check arguments for validity
   err = 0;
   if tmax <= 0 
      fprintf('tmax must be > 0\n'); err = 1;
   end
   if level ~= fix(level)
      fprintf('level is not an integer\n'); err = 1;
   else
      if level < 3 
         fprintf('level must be >= 3\n'); err = 1;
      end
   end
   if lambda <= 0
      fprintf('lambda must be >= 0\n'); err = 1;
   end
   if (x0l < 0.0) | (x0l > 1.0) 
      fprintf('x0l must be in range 0.0 - 1.0\n'); err = 1;
   end
   % Bail out if there were any errors in the arguments
   if err
      return
   end

   % Define x and t vectors and related mesh parameters
   nx = 2^level + 1;
   x = linspace(0.0, 1.0, nx);
   delx = x(2) - x(1);

   delt = lambda * delx;
   % Use round to guard against possible round-off effects
   nt = round(tmax / delt + 1);
   t = linspace(0.0, tmax, nt);
   delt = t(2) - t(1);

   if trace
      fprintf('advectloop: Dump of mesh parameters follows:\n');
      nx, delx, nt, delt1, delt
   end

   % Initialize solution array.  This is important for efficiency:
   % the code would execute much more slowly if rows were added to
   % u "dynamically", i.e. as the solution progressed.
   u = zeros(nt, nx);

   % Initialize first and second time steps using the function
   % advectl which returns the exact solution evaluated on the 
   % mesh x at time t (0 in the first call, delt in the second).
   u(1,:) = advectl(x, 0,    al, x0l, deltal);
   u(2,:) = advectl(x, delt, al, x0l, deltal);

   % Advance the discrete solution the required number of time steps
   % using the difference equations.
   if trace
      fprintf('advectloop: Evolution beginning\n');
   end
   for n = 2 : nt-1
      if rem(n, tracefreq) == 0 
         fprintf('advectloop: Step %d of %d\n', n, nt);
      end

      % Interior equations.
      for j = 2 : nx - 1
         u(n+1, j) = u(n-1, j) + lambda *  (u(n, j+1)  - u(n, j-1));
      end

      % Right boundary.
      u(n+1, nx) = 0;

      % Left boundary (first order approximation of advection equation).
      u(n+1, 1) = (u(n-1, 1) + lambda * u(n+1,2) / 2) / (1 + lambda / 2);
   end
end
