function output = csapi(x,y,xx)
%CSAPI Cubic spline interpolant with not-a-knot end condition.
%
%   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
%   data values may be vector-valued, in which case Y(:,j) is the j-th
%   value.
%
%   The alternative call
%
%   PP  = CSAPI(X,Y)
%
%   returns the ppform of the cubic spline interpolant instead, for
%   later use with FNVAL, FNDER, FNPLT, etc.
%
%   For example, 
%
%      values = csapi([-1:5]*(pi/2),[-1 0 1 0 -1 0 1], linspace(0,2*pi));
%
%   gives a surprisingly good fine sequence of values for the sine over its
%   period.
%
%   It is also possible to interpolate to data on a rectangular grid and
%   obtain values of the interpolant on a rectangular grid or at scattered
%   points, by the calls
%
%   VALUES = CSAPI({X1, ...,Xm},Y,XX)
%   or
%   PP = CSAPI({X1, ...,Xm},Y)
%
%   in which Y is expected to have size [d,length(X1),...,.length(Xm)]
%   (or [length(X1),...,.length(Xm)] if the function is to be scalar-valued),
%   and XX is either a list of m-vectors XX(:,j) or else a cell-array 
%   {XX1,...,XXm} specifying the m-dimensional grid at which to evaluate
%   the interpolant.
%
%   For example, the statements
%
%      x = -1:.2:1; y=-1:.25:1; [xx, yy] = ndgrid(x,y); 
%      z = sin(10*(xx.^2+yy.^2)); pp = csapi({x,y},z);
%      fnplt(pp)
%
%   produce the picture of an interpolant to a bivariate function. 
%   Use of MESHGRID instead of NDGRID here would produce an error.
%
%   See also CSAPE, SPAPI, SPLINE, NDGRID.

%   Carl de Boor Feb.25, 1989
%   cb : 22 feb 94 (use sparse matrices, handle vector-valued ordinates)
%   cb : 9 may '95 (use .' instead of ')
%   cb : 26oct97 (also handle gridded data)
%   cb : 28nov98 (handle sparse Y, mention NDGRID in help)
%   cb : 06feb00 (site/value instead of abscissa/ordinate)
%   cb : 19jul00 (correct misprint in help)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.11 $

if iscell(x)     % we are to handle gridded data
   
   m = length(x);
   sizey = size(y);
   switch length(sizey)
     case m  % grid values of a scalar-valued function
        if issparse(y), y = full(y); end 
        sizey = [1 sizey]; y = reshape(y, sizey); 
     case m+1
     otherwise
        error(['If X is a cell-array of length m, then Y must be an ', ...
               'm- or (m+1)-dimensional array.'])
   end
   
   v = y; sizev = sizey;
   for i=m:-1:1   % carry out coordinatewise interpolation
      [b,v,l,k] = ppbrk(csapi1(x{i}, reshape(v,prod(sizev(1:m)),sizev(m+1)) ));
      breaks{i} = b;
      sizev(m+1) = l*k; v = reshape(v,sizev);
      if m>1
         v = permute(v,[1,m+1,2:m]); sizev(2:m+1) = sizev([m+1,2:m]);
      end
   end
   % At this point, V contains the tensor-product pp coefficients;
   % It remains to make up the formal description:
   pp = ppmak(breaks, v);

else             % we have univariate data
   pp = csapi1(x,y);
end

if nargin==2
   output = pp;
else
   output = ppual(pp,xx);
end

function pp = csapi1(x,y)
%CSAPI1 Univariate cubic spline interpolant with not-a-knot end condition.

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

n = length(x); [xi,ind] = sort(x); xi = xi(:);
pp = [];
if n<2, error('There should be at least two data points.'), end
if all(diff(xi))==0, error('The data sites should be distinct.'), end

[yd,yn] = size(y); % if y happens to be a column matrix, change it to
                   % the expected row matrix.
if yn==1, yn = yd; yd = 1; y=reshape(y,1,yn); end

if yn~=length(x)
   error('Number of data sites and data values should match.')
end

yi=y(:,ind).'; dd = ones(1,yd);
dx=diff(xi); divdif = diff(yi)./dx(:,dd);

if (n==2), % the interpolant is a straight line
   pp=ppmak(xi.',[divdif.' yi(1,:).'],yd);
elseif (n==3), % the interpolant is a parabola
   yi(2:3,:)=divdif;
   yi(3,:)=diff(divdif)*(1/(xi(3)-xi(1)));
   yi(2,:)=yi(2,:)-yi(3,:)*dx(1);
   pp = ppmak( [xi(1),xi(3)], yi([3 2 1],:).' );
else % set up the sparse linear system for solving for the slopes at  xi .
   c = spdiags([ [dx(2:n-1);0;0] ...
                 [0;2*(dx(2:n-1)+dx(1:n-2));0] ...
                 [0;0;dx(1:n-2)] ], [-1 0 1], n, n);
     % the first two and last two equations are special:
   xi31 = xi(3)-xi(1); xin = xi(n)-xi(n-2);
   c(1,1:2) = [dx(2) xi31]; c(n,n-1:n) = [xin dx(n-2)];
   b = zeros(n,yd);
   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,dd).*divdif(1:n-2,:)+dx(1:n-2,dd).*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 ...   (protect current spparms settings)
   mmdflag = spparms('autommd');
   spparms('autommd',0);    % suppress pivoting
   s=c\b;
   spparms('autommd',mmdflag);

     %                          ... and convert to pp form

   c4 = (s(1:n-1,:)+s(2:n,:)-2*divdif(1:n-1,:))./dx(:,dd);
   c3 = (divdif(1:n-1,:)-s(1:n-1,:))./dx(:,dd) - c4;
   pp = ppmak(xi.', ...
     reshape([(c4./dx(:,dd)).' c3.' s(1:n-1,:).' yi(1:n-1,:).'],(n-1)*yd,4),yd);
end
