function [varargout] = SlopeLimiterModel( what, varargin )

  global SlopeLimiter_parameters ;

  switch ( what )

    case 'setup'
      expected( what, 1, nargin, 0, nargout ) ;
      SlopeLimiter_parameters . limiter = lower(  varargin{1} ) ;

    case 'eval_slope'
      expected( what, 1, nargin, 1, nargout ) ;
      varargout{1}  = SlopeLimiter( varargin{1} ) ;

    otherwise
      error(['SlopeLimiterModel bad argument: ' what]) ;
  end
end

function expected( what, ein, nin, eout, nout )
  if ein ~= nin-1 
    error( ['in SlopeLimiterModel call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
  end
  if eout ~= nout 
    error( ['in SlopeLimiterModel call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;  
  end
end

function DELTA = SlopeLimiter( US )
  global SlopeLimiter_parameters ;
  limiter = SlopeLimiter_parameters . limiter ;
  [r,c] = size(US) ;

  for k=1:r
    UL = US(r,1:end-2) ;
    UC = US(r,2:end-1) ;
    UR = US(r,3:end-0) ;

    P    = UR - UC ;
    Q    = UC - UL ;
    EPSI = 1E-10 * ones(size(Q)) ;

    % Compute RATIO=DUPW/DLOC of upwind to local changes 
    % some workaround to avoid division by zero are added
    R = sign(Q).*P./max(abs(Q),EPSI) ;
    S = sign(P).*Q./max(abs(P),EPSI) ;
    DELTA(r,:) = min( Q .* StandardLimiter( limiter, R ), ...
                      P .* StandardLimiter( limiter, S )  ) ;
  end
end
