function pp = CCSpp(x,f,fpl,fpr)
% usage: pp = CCSpp(x,f,fpl,fpr)
% description:  this function returns the clamped cubic spline
% structure with knots at the entries of the row vector of abscissas
% x and corresponding ordinates in the row vector f. The scalars
% fpl and fpr are values of f'(x) at left and right endpoints of x.

% local variables:
% n,m: 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
% M: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

n = length(x);
m= length(f);
if (m ~= n) % error check for size and order
  error('Size mismatch of abscissa and function ordinate');
end;
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
d = (f(2:n) - f(1:n-1))./h;
% now set up the system as in ClassroomNotes, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row
Coef(1, 1:2) = [h(1)/3, h(1)/6];
Rhs(1) = d(1) - fpl;
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) = [h(n-1)/6, h(n-1)/3]; 
Rhs(n) = fpr - d(n-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);

