%
% Compute the clamped cubic spline.
%
% function [a, b, c, d] = ccspline( x, f, fpa, fpb )
%
% input:
%        x:     vector containing the 
%               interpolation points
%
%        f:     vector containing the 
%               function values to be interpolated 
%
%        fpa:   value of f prime at left boundary
%
%        fpb:   value of f prime at right boundary
%
% output:
%
%        a:     vector containing the 
%               coefficients a(i)
%
%        b:     vector containing the 
%               coefficients b(i)
%
%        c:     vector containing the 
%               coefficients c(i)
%
%        d:     vector containing the 
%               coefficients d(i)
%
%        iflag  error flag
%               iflag = 0  spline is computed
%               iflag = 1  number of interpolation points and function
%                          values are not the same

%

function [a, b, c, d, iflag] = ccspline( x, f, fpa, fpb );

iflag = 0;

n = size(x(:),1);
if ( n ~= size(f(:),1) )
   iflag = 1;
   return
end

% allocate arrays
a = zeros(n,1);
b = zeros(n,1);
c = zeros(n,1);
d = zeros(n,1);
e = zeros(n,1);
h = zeros(n,1);

% compute interval length h(i) = x(i+1) - x(i);
h(1:n-1) = x(2:n) - x(1:n-1);

% Compute the   n by n  tridiagonal matrix  T
d(1)     = 2*h(1);
e(1)     = h(1);
d(2:n-1) = 2*( h(1:n-2)+h(2:n-1) );
e(2:n-1) = h(2:n-1);
d(n)     = 2*h(n-1);


% right hand side of the system (stored in c)
t1   =  3*(f(2) - f(1)) / h(1);
c(1) = t1 - 3*fpa;
for i = 2:n-1
    t2   =  3*(f(i+1) - f(i)) / h(i);
    c(i) = t2 - t1;
    t1   = t2;
end
c(n) = 3*fpb - t1;

%  solve  T c = c
[d, e, ierr] = tridiagLDL( d, e );
[c, ierr] = tridiagLDLsl( d, e, c );


%  compute a(i), d(i)
for i = 1:n-1
    a(i) = f(i);
    d(i) = (c(i+1) - c(i)) / (3*h(i));
end
a(n) = f(n);

%  compute b(i)
for i = 1:n-1
    b(i) = (a(i+1) - a(i))/h(i) - h(i)*(2*c(i) + c(i+1))/3;
end