function yi = interpBelow(x, y, xi, dimensions)  
%   YI = interpBelow(X, Y, XI, dimensions) interpolates to find YI,  
%   the values of the underlying function Y at the points in the array XI. 

%   Copyright 2008-2019 The MathWorks, Inc.
   
    n = length(x);
    ds = dimensions;
    
    if isscalar(ds)
        yMat = y;        
    else       
        prodDs = prod(ds);
        yMat = reshape(y,[prodDs n]);
    end

    if n<1
        error(message('Sldv:SldvDataInterp:NotEnoughPts'))
    end

    if ~isvector(x)
        error(message('Sldv:SldvDataInterp:Xvector'));
    end
        
    % Prefer column vectors for x
    xCol = x(:);

    xiCol = xi(:);

    if isscalar(ds)
        siz_yi = [length(xi) ds];
    else
        siz_yi = [ds length(xi)];
    end

    if ~isreal(x)
        error(message('Sldv:SldvDataInterp:ComplexX'))
    end

    if  ~isreal(xi)
        error(message('Sldv:SldvDataInterp:ComplexInterpPts'))
    end

    if (any(isnan(xCol)))
        error(message('Sldv:SldvDataInterp:NaNinX'));
    end


    % Start the algorithm
    % We now have column vector xCol, column vector or 2D matrix yMat and
    % column vector xiCol.

    h = diff(xCol);
    eqsp = (norm(diff(h),Inf) <= eps(norm(xCol,Inf)));
    if any(~isfinite(xCol))
        eqsp = 0; % if an INF in x, x is not equally spaced
    end
    if eqsp
        h = (xCol(n)-xCol(1))/(n-1);
    end


    if any(h < 0)
        [xCol,p] = sort(xCol);
        yMat = yMat(p,:);
        if eqsp
            h = -h;
        else
            h = diff(xCol);
        end
    end
    if any(h == 0)
        error(message('Sldv:SldvDataInterp:RepeatedValuesX'));
    end

    numelXi = length(xiCol);

    if ~eqsp && any(diff(xiCol) < 0)
        [xiCol,p] = sort(xiCol);
    else
        p = 1:numelXi;
    end

    % Find indices of subintervals, x(k) <= u < x(k+1),
    % or u < x(1) or u >= x(m-1).
    if isempty(xiCol)
        k = xiCol;
    else
        [~,k] = histc(xiCol,xCol);
        k(xiCol<xCol(1) | ~isfinite(xiCol)) = 1;
        k(xiCol>=xCol(n)) = n;
    end

    yiMat(:,p) = yMat(:,k);
    if isscalar(ds)
        yi = reshape(yiMat.',siz_yi);        
    else
        yi = reshape(yiMat,siz_yi);
    end

end