function z=hcubspl(knots,data,derivatives,points);
%            function z=hcubspl(knots,data,derivatives,points);
%            HERMITE INTERPOLATION BY PIECEWISE CUBICS
%% This file provides a function that will interpolate function values
%% and derivative values at a set of points by piecewise cubic polynomials.
%% This is somethimes called Hermite Cubic Spline Interpolation
%%
%% INPUT:
%%   knots - points at which the values are to be interpolated
%%   data - function values to be interpolated
%%   derivatives - derivative values to be interpolated
%%   points - the points at which the interpolating function is to be evaluated.
%%   The input should be all row vectors OR all column vectors
%% 
%%  OUTPUT:
%%    z - the values of the interpolating function at points

% Compute the lengths of intervals
n=length(knots)-1;
h=knots(2:n+1)-knots(1:n);
% Compute the divided differences g[x_j,x_{j+1}]
DDg=(data(2:n+1)-data(1:n))./h;
% Compute the values of the coefficients of the cubic polynomials
% on each interval
c0=data(1:n);
c1=derivatives(1:n);
c3=(derivatives(1:n)+derivatives(2:n+1)-2*DDg)./(h.^2);
c2=((DDg-derivatives(1:n))./h)-(c3.*h);

%% Evaluate the spline at given values
%% Here we determine which interval the point is in and evaluate the
%%spline using the coefficients for that interval
z=zeros(1,length(points));
for j=0:n-1;
 Ij=find(points>=knots(j+1)-eps&points<=knots(j+2)+eps);
  if isempty(Ij)==0
  z(Ij)=(((points(Ij)-knots(j+1))*c3(j+1)+c2(j+1)).*(points(Ij)-knots(j+1))+...
            c1(j+1)).*(points(Ij)-knots(j+1))+c0(j+1);
  end;
end;
