function output = csapi(x,y,xx)
% CSAPI	Cubic spline interpolation.
%
%      values = csapi(x,y,xx)
%
% returns the values at  xx  of the cubic spline interpolant to the
% given data  (x,y)  (using the not-a-knot end condition).
%      The alternative call
%
%      pp = csapi(x,y)
%
% returns the pp-form of the cubic spline interpolant instead, for later
% use with  fnval,  fnder, etc.

% C. de Boor / latest change: Feb.25, 1989
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.

%     Generate the cubic spline interpolant in pp form, depending on
% the number of data points.

n=length(x);[xi,ind]=sort(x);xi=xi(:);
if n<2,
   error('There should be at least two data points!')
elseif all(diff(xi))==0,
   error('The data abscissae should be distinct!')
elseif n~=length(y),
   error('Abscissa and ordinate vector should be of the same length!')
else   
   yi=y(ind);yi=yi(:);
   if (n==2), % the interpolant is a straight line
      pp=ppmak(xi',[diff(yi)./diff(xi) yi(1)],1);
   elseif (n==3), % the interpolant is a parabola
      yi(2:3)=diff(yi)./diff(xi);
      yi(3)=(yi(3)-yi(2))/(xi(3)-xi(1));
      yi(2)=yi(2)-yi(3)*(xi(2)-xi(1));
      pp = ppmak([xi(1),xi(3)],yi([3 2 1])',1);
   else % set up the linear system for solving for the slopes at  xi .
      dx=diff(xi);divdif=diff(yi)./dx;xi31=xi(3)-xi(1);xin=xi(n)-xi(n-2);
      c = diag([dx(2:n-1);0],-1)+2*diag([0;dx(2:n-1)+dx(1:n-2);0]) ...
                    + diag([0;dx(1:n-2)],1);
      c(1,1:2)=[dx(2) xi31];
      c(n,n-1:n) = [xin dx(n-2)];
      b=zeros(n,1);
      b(1)=((dx(1)+2*xi31)*dx(2)*divdif(1)+dx(1)^2*divdif(2))/xi31;
      b(2:n-1)=3*(dx(2:n-1).*divdif(1:n-2)+dx(1:n-2).*divdif(2:n-1));
      b(n)=...
      (dx(n-1)^2*divdif(n-2)+(2*xin+dx(n-1))*dx(n-2)*divdif(n-1))/xin;

        % solve for the slopes and convert to pp form
      s=c\b;
      c4=(s(1:n-1)+s(2:n)-2*divdif(1:n-1))./dx;
      c3=(divdif(1:n-1)-s(1:n-1))./dx - c4;
      pp=ppmak(xi',[c4./dx c3 s(1:n-1) yi(1:n-1)],1);
   end
   if nargin==2,
      output=pp;
   else
      output=ppual(pp,xx);
   end
end
