function [c,xi] = spline_interpolation_cubic(t,f,bc,df)
% SPLINE_INTERPOLATION_CUBIC: cubic spline interpolation
% [c,xi] = spline_interpolation_cubic(t,f,bc,df)
% interpolation of the data (t(j),f(j)) by a spline 
% with coefficients c and knots xi 
% boundary conditions bc: 
% 'natural', 'not-a-knot', or 'derivative' (slopes df(1:2) at the first and last point) 

% add 3 knots to left and right
M = [4 -3; 3 -2; 2 -1];
xi = [M*t(1:2); t; M(3:-1:1,:)*t([end, end-1])];

% form knot differences tau(k+j)-tau(k), j=1,2,3
h1 = diff(xi);
h2 = h1(1:end-1)+h1(2:end);
h3 = h2(1:end-1)+h1(3:end);

% assemble matrix by diagonals of B-spline values
dl = (h1(4:end-2).^2)./(h2(3:end-2).*h3(2:end-2));
dr = (h1(3:end-3).^2)./(h2(3:end-2).*h3(3:end-1));
d = 1-dl-dr;
A = diag([dl; 0],-1) + diag([0; d; 0]) + diag([0; dr],1);

% modify first and last row according to boundary conditions
% define right hand side

switch bc

case 'natural' 
   f = [0; f; 0];
   A(1,1:3) = [1/h3(2), -1/h3(2)-1/h3(3), 1/h3(3)]/h2(3);
   A(end,end-2:end) = [1/h3(end-2), ... 
      -1/h3(end-2)-1/h3(end-1), 1/h3(end-1)]/h2(end-2);

case 'not-a-knot'
   f = [0; f; 0];
   p1 = 1/(h1(4)*h2(3)*h3(3)) + 1/(h1(4)*h2(4)*h3(3)) + 1/(h1(5)*h2(4)*h3(3)); 
   p2 = 1/(h1(4)*h2(4)*h3(4)) + 1/(h1(5)*h2(4)*h3(4)) + 1/(h1(5)*h2(5)*h3(4));
   A(1,1) = -1/(h1(4)*h2(3)*h3(2));
   A(1,5) = -1/(h1(5)*h2(5)*h3(5)); 
   A(1,2) = -A(1,1)+p1; 
   A(1,3) = -p1-p2; 
   A(1,4) = -A(1,5)+p2;
   d = size(f,1)-5;
   p1 = 1/(h1(d+4)*h2(d+3)*h3(d+3)) + 1/(h1(d+4)*h2(d+4)*h3(d+3)) + ...
	1/(h1(d+5)*h2(d+4)*h3(d+3)); 
   p2 = 1/(h1(d+4)*h2(d+4)*h3(d+4)) + 1/(h1(d+5)*h2(d+4)*h3(d+4)) + ...
        1/(h1(d+5)*h2(d+5)*h3(d+4));
   A(end,end-4) = -1/(h1(d+4)*h2(d+3)*h3(d+2));
   A(end,end) = -1/(h1(d+5)*h2(d+5)*h3(d+5)); 
   A(end,end-3) = -A(end,end-4)+p1; 
   A(end,end-2) = -p1-p2; 
   A(end,end-1) = -A(end,end)+p2;

case 'derivative'
   f = [df(1); f; df(2)]; 
   A(1,1:3) = 3*[-h1(4)/h3(2), -h1(3)/h3(3)+h1(4)/h3(2), h1(3)/h3(3)]/h2(3); 
   A(end,end-2:end) = 3*[-h1(end-2)/h3(end-2), ...
      h1(end-2)/h3(end-2)-h1(end-3)/h3(end-1), ...
      h1(end-3)/h3(end-1)]/h2(end-2);
end

c = A\f;

