function output = spline(x,y,xx)

 %SPLINE Cubic spline data interpolation.

 %   PP = SPLINE(X,Y) provides the piecewise polynomial form of the

 %   cubic spline interpolant to the data values Y at the data sites X,

 %   for use with the evaluator PPVAL and the spline utility UNMKPP.

 %   X must be a vector.

 %   If Y is a vector, then Y(j) is taken as the value to be matched at X(j),

 %   hence Y must be of the same length as X  -- see below for an exception

 %   to this.

 %   If Y is a matrix or ND array, then Y(:,...,:,j) is taken as the value to

 %   be matched at X(j),  hence the last dimension of Y must equal length(X) --

 %   see below for an exception to this.

 %

 %   YY = SPLINE(X,Y,XX) is the same as  YY = PPVAL(SPLINE(X,Y),XX), thus

 %   providing, in YY, the values of the interpolant at XX.  For information

 %   regarding the size of YY see PPVAL.

 %

 %   Ordinarily, the not-a-knot end conditions are used. However, if Y contains

 %   two more values than X has entries, then the first and last value in Y are

 %   used as the endslopes for the cubic spline.  If Y is a vector, this

 %   means:

 %       f(X) = Y(:end-),  Df(min(X))=Y(),    Df(max(X))=Y(end).

 %   If Y is a matrix or N-D array with SIZE(Y,N) equal to LENGTH(X)+, then

 %   f(X(j)) matches the value Y(:,...,:,j+) for j=:LENGTH(X), then

 %   Df(min(X)) matches Y(:,:,...:,) and Df(max(X)) matches Y(:,:,...:,end).

 %

 %   Example:

 %   This generates a sine-like spline curve and samples it over a finer mesh:

 %       x = :;  y = sin(x);

 %       xx = :.:;

 %       yy = spline(x,y,xx);

 %       plot(x,y,'o',xx,yy)

 %

 %   Example:

 %   This illustrates the use of clamped or complete spline interpolation where

 %   end slopes are prescribed. In this example, zero slopes at the ends of an

 %   interpolant to the values of a certain distribution are enforced:

 %      x = -:; y = [ . 1.12 2.36 2.36 1.46 . . ];

 %      cs = spline(x,[ y ]);

 %      xx = linspace(-,,);

 %      plot(x,y,'o',xx,ppval(cs,xx),'-');

 %

 %   Class support for inputs x, y, xx:

 %      float: double, single

 %

 %   See also INTERP1, PCHIP, PPVAL, MKPP, UNMKPP.

 %   Carl de Boor --

 %   Copyright - The MathWorks, Inc.

 %   $Revision: 5.18.4.6 $  $Date: // :: $

 % Check that data are acceptable and, if not, try to adjust them appropriately

 [x,y,sizey,endslopes] = mychckxy(x,y);

 n = length(x); yd = prod(sizey);

 % Generate the cubic spline interpolant in ppform

 dd = ones(yd,); dx = diff(x); divdif = diff(y,[],)./dx(dd,:);

 if n==

    if isempty(endslopes) % the interpolant is a straight line

       pp=mkpp(x,[divdif y(:,)],sizey);

    else         % the interpolant is the cubic Hermite polynomial

       pp = pwch(x,y,endslopes,dx,divdif); pp.dim = sizey;

    end

 elseif n==&&isempty(endslopes) % the interpolant is a parabola

    y(:,:)=divdif;

    y(:,)=diff(divdif')'/(x()-x());

    y(:,)=y(:,)-y(:,)*dx();

    pp = mkpp(x([,]),y(:,[  ]),sizey);

 else % set up the sparse, tridiagonal, linear system b = ?*c for the slopes

    b=zeros(yd,n);

    b(:,:n-)=*(dx(dd,:n-).*divdif(:,:n-)+dx(dd,:n-).*divdif(:,:n-));

    if isempty(endslopes)

       x31=x()-x();xn=x(n)-x(n-);

       b(:,)=((dx()+*x31)*dx()*divdif(:,)+dx()^*divdif(:,))/x31;

       b(:,n)=...

       (dx(n-)^*divdif(:,n-)+(*xn+dx(n-))*dx(n-)*divdif(:,n-))/xn;

    else

       x31 = ; xn = ; b(:,[ n]) = dx(dd,[ n-]).*endslopes;

    end

    dxt = dx(:);

    c = spdiags([ [x31;dxt(:n-);] ...

         [dxt();*(dxt(:n-)+dxt(:n-));dxt(n-)] ...

         [;dxt(:n-);xn] ],[-  ],n,n);

    % sparse linear equation solution for the slopes

    mmdflag = spparms('autommd');

    spparms('autommd',);

    s=b/c;

    spparms('autommd',mmdflag);

    % construct piecewise cubic Hermite interpolant

    % to values and computed slopes

    pp = pwch(x,y,s,dx,divdif); pp.dim = sizey;

 end

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

 function [x,y,sizey,endslopes] = mychckxy(x,y)

 %CHCKXY check and adjust input for SPLINE and PCHIP

 %   [X,Y,SIZEY] = CHCKXY(X,Y) checks the data sites X and corresponding data

 %   values Y, making certain that there are exactly as many sites as values,

 %   that no two data sites are the same, removing any data points that involve

 %   NaNs, reordering the sites if necessary to ensure that X is a strictly

 %   increasing row vector and reordering the data values correspondingly,

 %   and reshaping Y if necessary to make sure that it is a matrix, with Y(:,j)

 %   the data value corresponding to the data site X(j), and with SIZEY the

 %   actual dimensions of the given values.

 %   This call to CHCKXY is suitable for PCHIP.

 %

 %   [X,Y,SIZEY,ENDSLOPES] = CHCKXY(X,Y) also considers the possibility that

 %   there are two more data values than there are data sites.

 %   If there are, then the first and the last data value are removed from Y

 %   and returned separately as ENDSLOPES. Otherwise, an empty ENDSLOPES is

 %   returned.  This call to CHCKXY is suitable for SPLINE.

 %

 %   See also PCHIP, SPLINE.

 %   Copyright - The MathWorks, Inc.

 % make sure X is a vector:

 if length(find(size(x)>))>

   error(message('MATLAB:chckxy:XNotVector'))

 end

 % ensure X is real

 if any(~isreal(x))

   error(message('MATLAB:chckxy:XComplex'))

 end

 % deal with NaN's among the sites:

 nanx = find(isnan(x));

 if ~isempty(nanx)

    x(nanx) = [];

    warning(message('MATLAB:chckxy:nan'))

 end

 n=length(x);

 if n<

   error(message('MATLAB:chckxy:NotEnoughPts'))

 end

 % re-sort, if needed, to ensure strictly increasing site sequence:

 x=x(:).';

 dx = diff(x);

 if any(dx<), [x,ind] = sort(x); dx = diff(x); else ind=:n; end

 if ~all(dx), error(message('MATLAB:chckxy:RepeatedSites')), end

 % if Y is ND, reshape it to a matrix by combining all dimensions but the last:

 sizey = size(y);

 while length(sizey)>&&sizey(end)==, sizey(end) = []; end

 yn = sizey(end);

 sizey(end)=[];

 yd = prod(sizey);

 if length(sizey)>

    y = reshape(y,yd,yn);

 else

    % if Y happens to be a column matrix, change it to the expected row matrix.

    if yn==

        yn = yd;

        y = reshape(y,,yn);

        yd = ;

        sizey = yd;

    end

 end

 % determine whether not-a-knot or clamped end conditions are to be used:

 nstart = n+length(nanx);

 if yn==nstart

    endslopes = [];

 elseif nargout==&&yn==nstart+

    endslopes = y(:,[ n+]); y(:,[ n+])=[];

    if any(isnan(endslopes))

       error(message('MATLAB:chckxy:EndslopeNaN'))

    end

    if any(isinf(endslopes))

        error(message('MATLAB:chckxy:EndslopeInf'))

    end

 else

    error(message('MATLAB:chckxy:NumSitesMismatchValues',nstart, yn))

 end

 % deal with NaN's among the values:

 if ~isempty(nanx)

     y(:,nanx) = [];

 end

 y=y(:,ind);

 nany = find(sum(isnan(y),));

 if ~isempty(nany)

    y(:,nany) = []; x(nany) = [];

    warning(message('MATLAB:chckxy:IgnoreNaN'))

    n = length(x);

    if n<

      error(message('MATLAB:chckxy:NotEnoughPts'))

    end

 end

function [x,y,sizey,endslopes] = mychckxy(x,y)

%CHCKXY check and adjust input for SPLINE and PCHIP

%   [X,Y,SIZEY] = CHCKXY(X,Y) checks the data sites X and corresponding data

%   values Y, making certain that there are exactly as many sites as values,

%   that no two data sites are the same, removing any data points that involve

%   NaNs, reordering the sites if necessary to ensure that X is a strictly

%   increasing row vector and reordering the data values correspondingly,

%   and reshaping Y if necessary to make sure that it is a matrix, with Y(:,j)

%   the data value corresponding to the data site X(j), and with SIZEY the

%   actual dimensions of the given values.

%   This call to CHCKXY is suitable for PCHIP.

%

%   [X,Y,SIZEY,ENDSLOPES] = CHCKXY(X,Y) also considers the possibility that

%   there are two more data values than there are data sites.

%   If there are, then the first and the last data value are removed from Y

%   and returned separately as ENDSLOPES. Otherwise, an empty ENDSLOPES is

%   returned.  This call to CHCKXY is suitable for SPLINE.

%

%   See also PCHIP, SPLINE.

%   Copyright 1984-2011 The MathWorks, Inc.

% make sure X is a vector:

if length(find(size(x)>1))>1

  error(message('MATLAB:chckxy:XNotVector'))

end

% ensure X is real

if any(~isreal(x))

  error(message('MATLAB:chckxy:XComplex'))

end

% deal with NaN's among the sites:

nanx = find(isnan(x));

if ~isempty(nanx)

   x(nanx) = [];

   warning(message('MATLAB:chckxy:nan'))

end

n=length(x);

if n<2

  error(message('MATLAB:chckxy:NotEnoughPts'))

end

% re-sort, if needed, to ensure strictly increasing site sequence:

x=x(:).';

dx = diff(x);

if any(dx<0), [x,ind] = sort(x); dx = diff(x); else ind=1:n; end

if ~all(dx), error(message('MATLAB:chckxy:RepeatedSites')), end

% if Y is ND, reshape it to a matrix by combining all dimensions but the last:

sizey = size(y);

while length(sizey)>2&&sizey(end)==1, sizey(end) = []; end

yn = sizey(end);

sizey(end)=[];

yd = prod(sizey);

if length(sizey)>1

   y = reshape(y,yd,yn);

else

   % if Y happens to be a column matrix, change it to the expected row matrix.

   if yn==1

       yn = yd;

       y = reshape(y,1,yn);

       yd = 1;

       sizey = yd;

   end

end

% determine whether not-a-knot or clamped end conditions are to be used:

nstart = n+length(nanx);

if yn==nstart

   endslopes = [];

elseif nargout==4&&yn==nstart+2

   endslopes = y(:,[1 n+2]); y(:,[1 n+2])=[];

   if any(isnan(endslopes))

      error(message('MATLAB:chckxy:EndslopeNaN'))

   end

   if any(isinf(endslopes))

       error(message('MATLAB:chckxy:EndslopeInf'))

   end

else

   error(message('MATLAB:chckxy:NumSitesMismatchValues',nstart, yn))

end

% deal with NaN's among the values:

if ~isempty(nanx)

    y(:,nanx) = [];

end

y=y(:,ind);

nany = find(sum(isnan(y),1));

if ~isempty(nany)

   y(:,nany) = []; x(nany) = [];

   warning(message('MATLAB:chckxy:IgnoreNaN'))

   n = length(x);

   if n<2

     error(message('MATLAB:chckxy:NotEnoughPts'))

   end

end