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

  NN = [ 50 100 200 400 800 ] ;

  for k=1:4
    
    N = NN(k) ;% 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 = Q0(Xc) ;
  
    % initial time and time step
    InitialTime = 0 ;
    FinalTime   = 20 ;
    StepTime    = 1 ;
    CFLratio    = 0.95 ;
  
    Qesatta = Q0(Xc - 20) ;

    [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                      Xnodes,       ...
                                      InitialTime,  ...
                                      FinalTime,    ...
                                      StepTime,     ...
                                      CFLratio,     ...
                                      @maxSpeed,    ...
                                      @doStepADER2 ) ;
                                  
    % plot the solution at final time step
    [XX,YY] = buildXYForPlot( Xnodes, Qframes{end} ) ;
    [X1,Y1] = buildXYForPlot( Xnodes, Qesatta ) ;
    plot( XX, YY, X1, Y1) ;
  
    err(k) =  max(abs(YY-Y1)) ;
    disp( sprintf( 'Errore %g', err(k) ) ) ;
  end
  
  disp( sprintf( 'Order %g', log(err(1)/err(2))/log(2) ) ) ;  
  disp( sprintf( 'Order %g', log(err(2)/err(3))/log(2) ) ) ;
  disp( sprintf( 'Order %g', log(err(3)/err(4))/log(2) ) ) ;

end

%
%
%
function Q = Q0( Xc ) 
  Q = 1+sin(Xc) ;
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
%
%  Given average states compute unlimited slopes
%
%  INPUTS:
%
%     Q = average values
%
%  OUTPUTS:
%
%     DELTA = slopes
%
%     The slopes are built by finite difference of sucessive averages.
%     In this way slopes associated to interfaces are used to copute cell slopes.
%     The slopes associated to a cell is the minimum slope between the two interface slopes.
%
function DELTA = SlopeNoLimiter( Q )
  % compute interface slopes, SLOPE(k) is the slopes of interface between cell k and k+1
  DELTA = (Q(:,3:end) - Q(:,1:end-2))./2 ;
end

%
% given the vector (matrix) a and b return
% res(k) = sign(a)*min(abs(a(k)),abs(b(k))) if a(k)*b(k)>0
% res(k) = 0 otherwise
%
function res = minmod( a, b )
  res = zeros( size(a) ) ;
  I1 = find( a > 0 & b > 0 ) ;
  I2 = find( a < 0 & b < 0 ) ;
  I3 = find( a.*b <= 0 ) ;
  res(I1) = min(a(I1),b(I1)) ;
  res(I2) = max(a(I2),b(I2)) ;
  res(I3) = 0 ;
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, @SlopeNoLimiter, @Leading0, @Leading1, @Source ) ;

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

end