%  Test Advancing routine by solving the problem
%
%   Q(x,t) + Q(x,t)  = 0    -1 <= x <= 1
%    t        x
%
%            /
%            |  0 if x < -1/2
%   Q(x,0) = |  1 if x > -1/2 and x < 1/2
%            |  0 if x > 1/2
%             \
%               +---+
%               |   |
%      ---------+   +----------
%
%  Cyclic boundary conditions
%
%  Last revision: 28 May 2009
%
function TestAder2LinearTransport

  N     = 100 ; % number of cells
  N4    = N/4 ;

  % setup intial conditions
  Qinit = [ zeros(1,N4) ones(1,N-2*N4) zeros(1,N4) ] ;
  
  % setup the boundary of the cells
  Xnodes = -1:2/N:1 ;
  
  % initial time and time step
  InitialTime = 0 ;
  FinalTime   = 20 ;
  StepTime    = 0.1 ;
  CFLratio    = 0.95 ;

  [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                    Xnodes,       ...
                                    InitialTime,  ...
                                    FinalTime,    ...
                                    StepTime,     ...
                                    CFLratio,     ...
                                    @maxSpeed,    ...
                                    @doStepADER2 ) ;

  [ Qframes1, Tframes1 ] = Advancing( Qinit,        ...
                                      Xnodes,       ...
                                      InitialTime,  ...
                                      FinalTime,    ...
                                      StepTime,     ...
                                      CFLratio,     ...
                                      @maxSpeed,    ...
                                      @doStepGodunov ) ;

                                % plot the solution at final time step
  [XX,YY] = buildXYForPlot( Xnodes, Qframes{end} ) ;
  [X1,Y1] = buildXYForPlot( Xnodes, Qframes1{end} ) ;
  plot( XX, YY, X1, Y1 ) ;
                             
end
%
% Return a vector with the maximum speed of each cell
%
function SPEED = maxSpeed( Q )
  SPEED = ones(size(Q)) ;
end

%
% Flux for linear transport
%
function F = Flux( Q )
  F = Q ;
end

%
% Jacobians of the flux (vectorial form)
%
function A = JFlux( Q )
  A = ones(size(Q)) ;
end
%
% Solve the Riemann problem for the linear transport equation exactly.
%
function USTAR = Riemann(UL, UR)
  USTAR = UL ;
end

function S = Source( Q )
  S = zeros(size(Q)) ;
end

function D0 = Leading0( QL, QR )
  D0 = Riemann(QL,QR) ;
end

function D1 = Leading1( DL, DR, D0 )
  D1 = Riemann(DL, DR) ;
end
%
% Advancing with ADER2 numerical scheme
%
function Qnew = doStepADER2( Time, DT, Xnodes, Q )

  % extrapolation
  Qe = [ Q(:,end-1) Q(:,end) Q Q(:,1) Q(:,2) ] ;

  [F,S] = ADER2Flux( Qe, Xnodes, DT, @Flux, @WENOSlopeLimiter, @Leading0, @Leading1, @Source ) ;

  DX   = Xnodes(2)-Xnodes(1) ;
  Qnew = Q - (DT/DX) .* ( F(:,2:end) - F(:,1:end-1) ) + DT * S ;

end

%
% Advancing with simple Godunov numerical scheme
%
function Qnew = doStepGodunov( Time, DT, Xnodes, Q )
  % extrapolation (with only one ghost cell by side)
  Qe = [ Q(:,end)  Q Q(:,1) ] ;
  % reconstruction no needed

  % compute numerical flux
  QL    = Qe(:,1:end-1) ; % cell on the left
  QR    = Qe(:,2:end)   ; % cell on the right
  QSTAR = Riemann(QL, QR) ;
  F     = Flux(QSTAR) ;
  S     = Source(Q) ;

  DX   = Xnodes(2) - Xnodes(1) ; % cells are assumed all of equal size
  Qnew = Q - (DT/DX) .* ( F(2:end) - F(1:end-1) ) + DT * S ;

end

