function output = csaps(x,y,p,xx)
% CSAPS	Cubic smoothing spline.
%
%      values = csaps(x,y,p,xx)
%
% returns the values at  xx  of the cubic smoothing spline for the given data
%  (x,y)  and depending on the smoothing parameter  p  from  [0,1] . For  p=0,
% this is the least-squares straight line fit to the data, while, on the other
% extreme, i.e., for  p=1 , this is the `natural' cubic spline interpolant.
%
%      The alternative call
%
%      pp = csaps(x,y,p)
%
% returns the pp-form of the cubic smoothing spline instead, for later use 
% with  fnval,  fnder, etc.

% C. de Boor / latest change: Sep.2, 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(:);
output=[];
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 smoothing spline is the straight line interpolant
      pp=ppmak(xi',[diff(yi)./diff(xi) yi(1)],1);
   else % set up the linear system for solving for the 2nd derivatives at  xi .
        % this is taken from (XIV.6)ff of the `Practical Guide to Splines'
        % with the diagonal matrix  D = eye(n,n)
      dx=diff(xi);divdif=diff(yi)./dx;
      odx=ones(n-1,1)./dx;
     if (n<50),
        diag([dx(1:n-2);0],-1) + 2*diag([0;dx(2:n-1)+dx(1:n-2);0]) ...
                      + diag([0;dx(2:n-1)],1);
        R=ans(2:n-1,2:n-1);
        diag([odx(1:n-2);0],-1) - diag([0;odx(2:n-1)+odx(1:n-2);0])...
                      + diag([0;odx(2:n-1)],1);
        Qt=ans(2:n-1,:);
        % solve for the 2nd derivatives 
        u=(6*(1-p)*Qt*Qt'+p*R)\diff(divdif);
     else, % use explicit Gauss elimination to handle the larger banded
           % matrices. The number 50 was picked without much thought.
           % For very large  n  , iteration might be a good alternative.
        % generate the five diagonals of Qt*Qt'
        P=zeros(n-2,6);
        oddx=odx(1:n-2)+odx(2:n-1);
        odx.*odx;
        P(:,3)=ans(1:n-2)+oddx.*oddx+ans(2:n-1);
        P(:,4)=-[(oddx(1:n-3)+oddx(2:n-2)).*odx(2:n-2);0];
        P(:,5)=[odx(2:n-3).*odx(3:n-2);0;0];
        % combine with the three diagonals of R to form (6*(1-p)*Qt*Qt'+p*R)
        P=(6*(1-p))*P;
        P(:,3)=P(:,3)+(2*p)*(dx(2:n-1)+dx(1:n-2));
        P(:,4)=P(:,4)+p*dx(2:n-1);
        P(:,2)=[0;P(1:n-3,4)];
        P(:,1)=[0;0;P(1:n-4,5)];
        % solve for the 2nd derivatives 
        % u=(6*(1-p)*Qt*Qt'+p*R)\diff(divdif);
        % by Gauss elimination without pivoting
        P(:,6)=diff(divdif);
        for j=1:n-4;
           P(j,4:6)=(1/P(j,3))*P(j,4:6);
           P(j+1,[3 4 6])=P(j+1,[3 4 6])-P(j+1,2)*P(j,[4:6]);
           P(j+2,[2 3 6])=P(j+2,[2 3 6])-P(j+2,1)*P(j,[4:6]);
        end
        j=n-3;
           P(j,4:6)=(1/P(j,3))*P(j,4:6);
           P(j+1,[3 4 6])=P(j+1,[3 4 6])-P(j+1,2)*P(j,[4:6]);
        j=n-2;
           P(j,6)=P(j,6)/P(j,3);
        j=n-3;
           P(j,6)=P(j,6)-P(j,4)*P(j+1,6);
        for j=n-4:-1:1;
           P(j,6)=P(j,6)-P(j,[4 5])*P(j+[1 2],6);
        end
        u=P(:,6);
     end
        % ... and convert to pp form
      % Qt'*u=diff([0;diff([0;u;0])./dx;0])
      yi = yi - 6*(1-p)*diff([0;diff([0;u;0])./dx;0]);c3 = [0;p*u;0];
      c2=diff(yi)./dx-dx.*(2*c3(1:n-1)+c3(2:n));
      pp=ppmak(xi',[diff(c3)./dx,3*c3(1:n-1),c2,yi(1:n-1)],1);
   end
   if nargin==3,
      output=pp;
   else
      output=ppual(pp,xx);
   end
end
 
