%  General Advancing procedure used to test numerical solver for hyperbolic equations
%
%  Name of program: Advancing
%
%  Developed by: Enrico Bertolazzi
%                Department of Mechanical and Structural Engineering
%                enrico.bertolazzi@ing.unitn.it
%                http://www.ing.unitn.it/~bertolaz/
%
%  Last revision: 2 June 2009
%
%  INPUTS:
%
%     Qinit       = initial cell average values
%     Xnodes      = boundary of the cells
%     InitialTime = starting simulation time
%     FinalTime   = end of simulation time
%     StepTime    = the time difference at which the solution is saved
%     CFLratio    = the Courant–Friedrichs–Lewy ratio at which the solver can advance
%
%     maxSpeed    = function with one vector argument which compute maximum speed signal for each cell.
%                   Usage: feval( maxSpeed, Q ).  
%
%     doStep      = function which perform step of the numerical scheme
%                   Usage: feval( doStep, Time, DT, DX, XC, Q ).
%                         Time = Actual time
%                         DT   = advancing time step
%                         DX   = cells size
%                         XC   = cells center
%                         Q    = cells actual state
%
%  OUTPUTS:
%
%     Tframes     = vector with the time values of saved frame
%     Qframes     = cell array with the solution frames computed
%                   Qframes{1}   = contains the solution at time Tframes(1)
%                   Qframes{2}   = contains the solution at time Tframes(2)
%
%                   ...
%
%                   Qframes{end} = contains the solution at time Tframes(end)
%
function [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                           Xnodes,       ...
                                           InitialTime,  ...
                                           FinalTime,    ...
                                           StepTime,     ...
                                           CFLratio,     ...
                                           maxSpeed,     ...
                                           doStep )

  % Initial conditions are copied in the vector Q
  Q = Qinit ;

  % Compute cell size DX, the size of the cells are assumed all equals.
  DX = Xnodes(2) - Xnodes(1) ;

  % save initial frame
  Qframes{1} = Q ;
  Tframes(1) = InitialTime ;

  kframe   = 1 ; % number of saved framed 
  Time     = InitialTime ; % initial time
  NextTime = InitialTime + StepTime ; % Time to save next frame

  % Time marching procedure, a maximum of 10000 steps are allowed
  for NITER=1:10000

    % Courant-Friedrichs-Lewy (CFL) condition imposed
    DT = CFLratio * DX / max( feval( maxSpeed, Q ) ) ;
        
    % reduce DT if necessary for the last step
    if Time+DT >= FinalTime
      DT = FinalTime - Time ;
      % perform last step
      Q = feval( doStep, Time, DT, Xnodes, Q ) ;
      % save last frame
      kframe          = kframe+1 ;
      Qframes{kframe} = Q ;
      Tframes(kframe) = FinalTime ;
      % exit loop
      break ;
    end

    % adjust time step if necessary to match the time grid
    DT1 = NextTime - Time ;
    while DT >= DT1 % if DT contains a grid step perform smaller step

      % perform step
      Qnew = feval( doStep, Time, DT1, Xnodes, Q ) ;

      % Solution frame is saved
      kframe          = kframe+1 ;
      Qframes{kframe} = Qnew ;
      Tframes(kframe) = NextTime ;

      % check if we have reached the final mesh point
      if FinalTime < NextTime 
        break ;
      end ;

      % update next time frame
      NextTime = NextTime + StepTime ;
      DT1      = NextTime - Time ;
    end

    % perform step
    Q = feval( doStep, Time, DT, Xnodes, Q ) ;

    % update time
    Time = Time + DT ;
  end

end
