%  Given average states compute WENO limited slopes
%
%  Name of program: WENOSlopeLimiter.m
%
%  Developed by: Enrico Bertolazzi
%                Department of Mechanical and Structural Engineering
%                enrico.bertolazzi@ing.unitn.it
%                http://www.ing.unitn.it/~bertolaz/
%
%  This code is inspired on the library
%    NUMERICA
%    A Library of Source Codes for Teaching,
%    Research and Applications,
%    by E. F. Toro
%    Published by NUMERITEK LTD,
%    Website: www.numeritek.coM
%
%  The algorithm are described in the book:
%
%    E. F. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics"
%                 Springer-Verlag, 1997, Second Edition, 1999
%
%  Last revision: 27 May 2009
%
%  The scheme can be visualized in the following way
% 
%   +===+  The cell where the average value Qi is "stored".
%
%   +...+  Ghost cell. These cells are used to extrapolate values for the
%          construction of the bounary conditions. Two ghost cell for boundary
%          are enough for all the considered fluxes.
%
%   +xxx+  Uncomputed values. These cells do not cotains valid values.
%          Are used to avoid shifting of indices or complex indicization.
%
%   Averages states
%
%   +...+...+===+===+===+===+===    ===+===+...+...+
%     Q1  Q2  Q3                                QN+4
%
%   Limited slopes, the slopes S1 and SN+4 are not computed (indeed it is not possible)
%
%   +xxx+===+===+===+===+===+===    ===+===+===+xxx+
%         S1  S2  S3                        SN+1
%
%  INPUTS:
%
%     Q = average values
%
%  OUTPUTS:
%
%     SLOPE = limited 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 = WENOSlopeLimiter( Q )

  % compute interface slopes, SLOPE(k) is the slopes of interface between cell k and k+1
  SLOPE = Q(:,2:end) - Q(:,1:end-1) ;

  % left slopes
  LSLOPE = SLOPE(:,1:end-1);

  % right slopes
  RSLOPE = SLOPE(:,2:end);

  % find minimum slopes
  DELTA = minmod( LSLOPE, RSLOPE ) ;

  % is equivalent to the loop
  %
  % for k=1:lenght(SLOPE)
  %   if abs(LSLOPE(k)) < abs(RSLOPE(k))
  %     SLOPE(k) = LSLOPE(k);
  %   else
  %     SLOPE(k) = RSLOPE(k);
  %   end
  % end

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
