function [ t, y, schemeData, stats ] = ...
                            odeCFL1(schemeFunc, tspan, y0, options, schemeData)
% odeCFL1: integrate a CFL constrained ODE (eg a PDE by method of lines).
%
% [ t, y, schemeData, stats ] = ...
%                           odeCFL1(schemeFunc, tspan, y0, options, schemeData)
%
% Integrates a system forward in time by CFL constrained timesteps
%   using a first order forward Euler scheme
%   (which happens to be the first order TVD RK scheme).
%
% parameters:
%   schemeFunc	 Function handle to a CFL constrained ODE system
%                  (typically an approximation to an HJ term, see below).
%   tspan        Range of time over which to integrate (see below).
%   y0           Initial condition vector 
%                  (typically the data array in vector form).
%   options      An option structure generated by odeCFLset 
%                  (use [] as a placeholder if necessary).
%   schemeData   Structure passed through to schemeFunc.
%
%
%   t            Output time(s) (see below).
%   y            Output state (see below).
%   schemeData   Output version of schemeData (see below).
%   stats        Structure containing integration statistics (see below).
%
% A CFL constrained ODE system is described by a function with prototype
%
%        [ ydot, stepBound ] = schemeFunc(t, y, schemeData)
%
%   where t is the current time, y the current state vector and
%   schemeData is passed directly through.  The output stepBound
%   is the maximum allowed time step that will be taken by this function
%   (typically the option parameter factorCFL will choose a smaller step size).
%
% The time interval tspan may be given as 
%   1) A two entry vector [ t0 tf ], in which case the output will 
%      be scalar t = tf and a row vector y = y(tf).
%   2) A vector with three or more entries, in which case the output will
%      be column vector t = tspan and each row of y will be the solution
%      at one of the times in tspan.  Unlike Matlab's ode suite routines,
%      this version just repeatedly calls version (1), so it is not
%      particularly efficient.
%
% Note that using this routine for integrating HJ PDEs will usually
%   require that the data array be turned into a vector before the call
%   and reshaped into an array after the call.  Option (2) for tspan should 
%   not be used in this case because of the excessive memory requirements
%   for storing solutions at multiple timesteps.
%
% The output version of schemeData will normally be identical to the input
%   version, and therefore can be ignored.  However, if a PostTimestep 
%   routine is used (see odeCFLset) then schemeData may be modified during 
%   integration, and the version of schemeData at tf is returned in this
%   output argument.
%
% The output parameter stats is a structure containing the following fields
%   (if tspan contains only two entries):
%
%   .numSteps    The number of full timesteps taken to reach the final time.
%   .numSubsteps The number of substeps taken for each timestep
%                  (for odeCFL1 this will always be 1).
%   .cpuTime     The computational time taken to reach the final time
%                  (as measured by Matlab's cputime command).
%   .stepTimes   The timestamps of the intermediate steps
%                  (a vector of length numSteps).
%   .stepBounds  The CFL bounds on step lengths for each of the 
%                  intermediate steps (a vector of length numSteps).
%
% If tspan contains three or more entries, then stats is a cell vector
%   of the same length as tspan.  Then stats{i} contains
%
%   i = 1        The empty vector [].
%   otherwise    The structure described above for the integration
%                  which finished at final time tspan(i).

% Copyright 2004 Ian M. Mitchell (mitchell@cs.ubc.ca).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Ian Mitchell, 5/14/03.
% Calling parameters modified to more closely match Matlab's ODE suite
%   Ian Mitchell, 2/6/04.
% Added stats output parameter.  Ian Mitchell 8/18/04.

  %---------------------------------------------------------------------------
  % How close (relative) do we need to be to the final time?
  small = 100 * eps;
  
  %---------------------------------------------------------------------------
  % Make sure we have the default options settings
  if((nargin < 4) | isempty(options))
    options = odeCFLset;
  end
  
  %---------------------------------------------------------------------------
  numT = length(tspan);
  
  %---------------------------------------------------------------------------
  % If we were asked to integrate forward to a final time.
  if(numT == 2)

    t = tspan(1);
    y = y0;
    steps = 0;
    startTime = cputime;
    
    while(tspan(2) - t > small * abs(tspan(2)))

      % Approximate the derivative and CFL restriction.
      [ ydot, stepBound ] = feval(schemeFunc, t, y, schemeData);

      % CFL bound on timestep, but not beyond the final time.
      deltaT = min([ options.factorCFL * stepBound; ...
                     tspan(2) - t; options.maxStep ]);
      y = y + deltaT * ydot;
      t = t + deltaT;

      steps = steps + 1;

      if(nargout > 3)
        % The caller wants integration statistics.

        if(steps == 1)
          % Preallocate space for the statistics structure.
          numStepsEstimate = 2 * ceil((tspan(2) - tspan(1)) / deltaT);
          stats.stepTimes = zeros(numStepsEstimate, 1);
          stats.stepBounds = zeros(numStepsEstimate, 1);
        end

        % Store this intermediate step's statistics.
        stats.stepTimes(steps) = t;
        stats.stepBounds(steps) = stepBound;
      end

      % If there is a post-timestep routine, call it.
      if(~isempty(options.postTimestep))
        [ y schemeData ] = feval(options.postTimestep, t, y, schemeData);
      end

      % If we are in single step mode, then do not repeat.
      if(strcmp(options.singleStep, 'on'))
        break;
      end
    
    end
    
    endTime = cputime;

    if(nargout > 3)
      % Complete the integration statistics.
      stats.numSteps = steps;
      stats.cpuTime = endTime - startTime;

      % This statistic depends only on the integration scheme.
      stats.numSubsteps = 1;

      % Truncate the likely overallocated intermediate step stat vectors.
      stats.stepTimes = stats.stepTimes(1:steps);
      stats.stepBounds = stats.stepBounds(1:steps);
    end

    if(strcmp(options.stats, 'on'))
      fprintf('\t%d steps in %g seconds from %g to %g\n', ...
              steps, endTime - startTime, tspan(1), t);
    end

  %---------------------------------------------------------------------------
  % If we were asked for the solution at multiple timesteps,
  %   just use recursion over each pair of timesteps.
  elseif(numT > 2)
    t = reshape(tspan, numT, 1);
    y = zeros(numT, length(y0));
    y(1,:) = y0';
    if(nargout > 3)
      stats = cell(numT, 1);
      stats{1} = [];
    end
    for i = 2 : numT
      [ t(i), yout, schemeData, stats{i} ] = ...
         odeCFL1(schemeFunc, [ t(i-1), t(i) ], y(i-1,:)', schemeData, options);
      y(i,:) = yout';
    end
    
  %---------------------------------------------------------------------------
  else
    error('tspan must contain at least two entries');
  end
