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

  global A_jacobian ;
  
  A_jacobian = [ 0 1 ; 1 0 ] ;

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

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

  [ 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(1,size(Q,2)) ;
end

function F = Flux( Q )
  F = JFlux( Q ) * Q ;
end

function A = JFlux( Q )
  global A_jacobian ;
  A = A_jacobian ;
end

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

function D0 = Leading0( QL, QR )
  global A_jacobian ;
  D0 = LinearRiemannProblem( QL, QR, A_jacobian ) ;
end

function D1 = Leading1( QL, QR, D0 )
  global A_jacobian ;
  D1 = LinearRiemannProblem( QL, QR, A_jacobian ) ;
end
%
% Advancing with simple Godunov 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 = Leading0(QL, QR) ;
  F     = Flux( QSTAR ) ;

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

end