%  Test Advancing routine by solving the problem
%
%   Q (x,t) + Q(x,t) Q (x,t) = 0    -pi <= x <= pi
%    t                x
%
%   Q(x,0) = 1+sin(x)
%
%  Cyclic boundary conditions
%
%  Last revision: 28 May 2009
%
function TestAder2Burger

  N = 100 ; % number of cells

  % setup the boundary of the cells
  DX     = 2*pi/N ;
  Xnodes = [-pi:DX:pi] ;

  % setup intial conditions
  Xc    = (Xnodes(1:end-1)+Xnodes(2:end))/2 ; % center of the cells
  Qinit = 1+sin(Xc) ;

  % initial time and time step
  InitialTime = 0 ;
  FinalTime   = 10 ;
  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 = abs(Q) ;
end

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

function F = Flux( Q )
  F = (Q .* Q) ./ 2 ;
end
%
% Solve the Riemann problem for the inviscid Burgers equation exactly.
%
function USTAR = Riemann(UL, UR)
  % UL         Left data state
  % UR         Right data state
  % S          Shock speed
  % USTAR      Sampled state
  ISW = find ( UL <  UR ) ;
  IRW = find ( UL >= UR ) ;
  % Solution is a shock wave
    % Compute shock speed S
    S(ISW) = 0.5*(UL(ISW) + UR(ISW)) ;
    % Sample the state along the t-axis
    IGE = find( S(ISW) >= 0 ) ;
    ILT = find( S(ISW) <  0 ) ;
    USTAR(IGE) = UL(IGE) ;
    USTAR(ILT) = UR(ILT) ;
  % Solution is a rarefaction wave.
    % There are 3 cases:
    I1 = find( UL >= 0 ) ;
    I2 = find( UR <= 0 ) ;
    I3 = find( UL < 0 & UR > 0 ) ;
    USTAR(I1) = UL(I1) ; % Right supersonic rarefaction
    USTAR(I2) = UR(I2) ; % Left supersonic rarefaction
    USTAR(I3) = 0 ;      % Transonic rarefaction
end

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

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

function D1 = Leading1( QL, QR, D0 )
  D1 = zeros(size(QL)) ;
  I1 = find( D0 >  0 ) ;
  I2 = find( D0 <= 0 ) ;
  D1(I1) = QL(I1) ;
  D1(I2) = QR(I2) ;
  %for k=1:size(QL,2)
  %  D1(k) = LinearRiemannProblem( QL(k), QR(k), A(k) ) ;
  %end
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 = Riemann(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

