function [spline_param,yint] = naturalcubicsplineINT(x, y, xint)
%
%  function spline_param = naturalcubicsplineINT(x, y, xint)
%
%  Use Algorithm 3.4 on pages 142-143 of Burden and Faires to construct
%  the cubic spline interpolant.
%
%  Input parameters:
%    x     vector of x coordinates of given points.
%    y     vector of y coordinates of given points.
%    xint  the x coordinate at which y is to be interpolated
%
%
%  Output parameters:
%    spline_param  the parameters for the natural cubic spline
%                  in the order
%                  [a_0,b_0,c_0,d_0,....,a_{n-1},b_{n-1},c_{n-1},d_{n-1}]
%    yint          the interpolated value at xint
%
%  Kartik Sivaramakrishnan
%  Department of  Mathematics
%  North Carolina State University
%  September 15, 2007
%

n = length(x); % Number of points

for i=1:n,
  a(i) = y(i);
end 

a = a(:);

for i=1:n-1,
  h(i) = x(i+1)-x(i);
end

h = h(:);
 
% Set up the rhs of the linear equations  
alpha = zeros(n,1);
for i=2:n-1,
   alpha(i) = 3*(a(i+1)-a(i))/h(i) - 3*(a(i)-a(i-1))/h(i-1);
end   
 
% Set up the coefficient matrix for the linear equations
A = zeros(n);
% First set the superdiagonal of coefficient matrix A
for i=2:n,
  if i == 2,
    A(i-1,i) = 0;
  else
    A(i-1,i) = h(i-1);
  end
end

% Next set the subdiagonal of coefficient matrix A
for i=2:n,
  if i == n,
    A(i,i-1) = 0;
  else
    A(i,i-1) = h(i-1);
  end
end  

% Finally set the diagonal entries of A
A(1,1) = 1;
A(n,n) = 1;
for i=2:n-1,
  A(i,i) = 2*(h(i-1)+h(i));
end


% Solve Ac=alpha for c
c = A\alpha;

c = c(:);

for i=n-1:-1:1,
  b(i) = (a(i+1)-a(i))/h(i) - h(i)*(c(i+1)+2*c(i))/3;
  d(i) = (c(i+1)-c(i))/(3*h(i));
end 

b=b(:);
d=d(:);

for i=1:n-1,
  spline_param(4*(i-1)+1) = a(i);
  spline_param(4*(i-1)+2) = b(i);
  spline_param(4*(i-1)+3) = c(i);
  spline_param(4*(i-1)+4) = d(i);
end 

spline_param = spline_param(:);

% Find the interval in which xint lies
if nargout > 1,
  if nargin < 3,
     fprintf('Please enter xint as 3rd input argument');
     return;
  end   
  index = find(x-xint > 0);
  if isempty(index) == 0,
    index = index(1)-1;
  else
    index = n-1;
  end
  val = (xint-x(index));
  yint = spline_param(4*(index-1)+1) + spline_param(4*(index-1)+2)*val + spline_param(4*(index-1)+3)*(val^2) + spline_param(4*(index-1)+4)*(val^3);
end