function pp = NKCSpp(xx,ff)
% usage: pp = NKCSpp(x,f)
% description: this function returns the not-a-knot cubic spline
% with nodes at the entries of the row vector of abscissas x and
% corresponding ordinates in the row vector f, excluding the 
% second and penultimate entries of each.

% local variables:
% x,f: vectors with not-a-knot entries removed
% nn,m,n: size of node and value arrays
% j: index variable
% h: stepsize vector
% d: divided difference vector
% Coef: coefficient matrix for the system
% Rhs: right hand side vector
% h0, h1, hnm2, hnm1: extra step sizes
% M: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

nn = length(xx);
m= length(ff);
if (m ~= nn | nn<4)  % error check for size and order
  error('Size mismatch of abscissa and function ordinate or size problem');
end;
x = [xx(1), xx(3:nn-2), xx(nn)];
f = [ff(1), ff(3:nn-2), ff(nn)];
n = nn - 2; % adjust size to reflect size of x
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
d = (f(2:n) - f(1:n-1))./h;
hh = xx(2:nn) - xx(1:nn-1);
% now set up the system as in ClassNotes, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row 
Coef(1, 1:2) = [hh(1)+2*hh(2), 2*hh(1)+hh(2)];
Rhs(1) = 6*(ff(1)/hh(1) - h(1)*ff(2)/(hh(1)*hh(2)) + ff(3)/hh(2));
for j = 2:(n-1) % interior rows
  Coef(j, (j-1):(j+1)) = [h(j-1)/6, (h(j-1)+h(j))/3, h(j)/6];
  Rhs(j) = d(j) - d(j-1);
end;
% last row 
Coef(n, (n-1):n) = [hh(nn-2)+2*hh(nn-1), 2*hh(nn-2)+hh(nn-1)]; 
Rhs(n) = 6*(ff(nn-2)/hh(nn-2) - h(n-1)*ff(nn-1)/(hh(nn-2)*hh(nn-1)) + ff(nn)/hh(nn-1));
% solve for M
M = Coef\Rhs;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai*(x-xi)^3+bi*(x-xi)^2+ci*(x-xi)+di 
% we deduce from si''(x)=6*ai*(x-xi)+2*bi that
% bi = Mi/2 and ai = (M(i+1)-Mi)/(6*hi)
coefs(:,2) = M(1:n-1)/2;
coefs(:,1) = (M(2:n) - M(1:n-1))./(6*h');
% evaluation gives di = yi
coefs(:,4) = f(1:n-1)';
% evaluation of s' gives formula for ci
coefs(:,3) = d - ((M(2:n) + 2*M(1:n-1))'/6).*h;
pp = PP(4,n,x,coefs);
