function yi = interp1cubic(x,y,xi)
% INTERP1CUBIC 1-D interpolation with cubic interpolant
%    YI = INTERP1(X,Y,XI) interpolates to find YI, the values of the
%    underlying function Y at the points in the array XI. X and Y must
%    be vectors of length N.
%    The interpolant is a fouth order accurate approximation to smooth
%    data, composed of local cubic polynomials using 4 data points each.

% 02/2015 by Benjamin Seibold
%            http://www.math.temple.edu/~seibold/

n = length(x); % number of data points
xc = [-inf,x(3:n-2),inf]; % cell boundaries (n-3 cells)
yi = xi;
for j = 1:n-3 % index of cell
    ind = xc(j)<xi&xi<=xc(j+1); % evaluation points in that cell
    xic = xi(ind);
    yi(ind) = ...
        y(j)*((xic-x(j+1)).*(xic-x(j+2)).*(xic-x(j+3)))./...
        ((x(j)-x(j+1)).*(x(j)-x(j+2)).*(x(j)-x(j+3)))+...
        y(j+1)*((xic-x(j)).*(xic-x(j+2)).*(xic-x(j+3)))./...
        ((x(j+1)-x(j)).*(x(j+1)-x(j+2)).*(x(j+1)-x(j+3)))+...
        y(j+2)*((xic-x(j)).*(xic-x(j+1)).*(xic-x(j+3)))./...
        ((x(j+2)-x(j)).*(x(j+2)-x(j+1)).*(x(j+2)-x(j+3)))+...
        y(j+3)*((xic-x(j)).*(xic-x(j+1)).*(xic-x(j+2)))./...
        ((x(j+3)-x(j)).*(x(j+3)-x(j+1)).*(x(j+3)-x(j+2)));
end
