%
% Compute Ader Flux
%
function [F,S] = ADER2Flux( Q,            ...
                            Xnodes,       ...
                            DT,           ...
                            Flux,         ...
                            SlopeLimiter, ...
                            Leading0,     ...
                            Leading1,     ...
                            Source )

  nr = size(Q,1) ;
  DX = Xnodes(2)-Xnodes(1) ;

  % Reconstruct phase
    % construct limited slopes (given N+2 averages return N slopes)
    DQ = feval( SlopeLimiter, Q ) ;

    % extrapolate from the left
    QL = Q(:,2:end-2) + DQ(:,1:end-1)./2 ;

    % extrapolate from the right
    QR = Q(:,3:end-1) - DQ(:,2:end)./2 ;

    % compute slopes
    DELTA = DQ ./ DX ;

    % left slopes
    DELTAL = DELTA(:,1:end-1) ;

    % right slopes
    DELTAR = DELTA(:,2:end) ;

  % End of Reconstruct phase

  % Compute leading term
    D0 = Leading0( QL, QR ) ;
    D1 = Leading1( DELTAL, DELTAR, D0 ) ;
  % End Compute leading term

  % Evaluate ADER2 state
    D0L = D0 - D1 ./2 ;
    D0R = D0 + D1 ./2 ;
    %
    % In order to eliminate the loop 
    %
    %Qader = D0 ;
    %for k=1:size(Qader,2)
    %  A0         = feval( JFlux, D0(:,k) ) ;
    %  Qader(:,k) = Qader(:,k) + (DT/2) * (feval(Source, D0(:,k)) - A0 * D1(:,k)) ;
    %end
    %
    % The second order approxination A0 * D1 ~ Flux( D0+D1/2 ) - Flux(D0-D1/2)
    %
    Qader = D0 + (DT/2) * ( feval( Source, D0 ) - feval( Flux, D0R ) + feval( Flux, D0L ) ) ;
    F     = feval( Flux, Qader) ;
  % End Evaluate numerical flux

  % Evaluate source
    Qs = Q(:,3:end-2) ; % consider only "true" cells
    DD = DELTA(:,2:end-1) ./2 ; % consider only half of the slope of "true" cells
    %
    % In order to eliminate the loop 
    %
    % for k=1:size(Qs,2)
    %   Qs(:,k) = Qs(:,k) + (DT/2) * ( feval( Source, Qs(:,k) ) - feval( JFlux, Qs(:,k) )*DELTA(:,k+1) ) ;
    % end
    %
    % The second order approxination JFlux(Qs) DELTA = Flux( Qs+DELTA/2 ) - Flux(Qs-DELTA/2)
    %
    DL = Qs - DD ;
    DR = Qs + DD ;  
    
    Qs = Qs + (DT/2) * ( feval( Source, Qs ) - feval( Flux, DR ) +  feval( Flux, DL ) ) ;
    S  = feval( Source, Qs ) ;
  % End Evaluate source

end
