function s = fspline(xd,yd,dyd,x)

%  clamped (or fixed) cubic spline at equally spaced data points
%  xd, yd : data vectors (they must have the same length)
%  dyd : a 2-vector containing the values of first derivative at two endpoints
%  x : points where s(x) is evaluated

% some additional information is at the end of the file

nd=length(xd);
ndd=length(yd);
% errors
if nd < 3
	error('You need at least three data points');
end
if nd ~= ndd 
	error('xd and yd must have the same length');
end


a=xd(1);
b=xd(nd);
np=length(x);
h=xd(2)-xd(1);

d=4*ones(nd,1);
dT=ones(nd,1);
dT(1)=2;
dB=ones(nd,1);
dB(nd)=2;
for j=1:nd
	z(j)=6*yd(j);
end;
z(1)=z(1)+2*h*dyd(1);
z(nd)=z(nd)-2*h*dyd(2);
w=tri(d,dB,dT,z); 
ww=[w(2)-2*h*dyd(1) w w(nd-1)+2*h*dyd(2)];

% evaluate interpolation fn
for ix=1:np
	sum=0;
	for k=1:nd+2
		xk=a-h+h*(k-1);
		xx=(x(ix)-xk)/h;
		sum=sum+ww(k)*bbspline(xx);
	end;
	s(ix)=sum;
end;


function y=bbspline(x)
% Calculate the value of a cubic B-spline at point x
x=abs(x) ;
if x>2,
  y=0 ;
else
  if x>1,
    y=(2-x)^3/6 ;
  else
    y=2/3-x^2*(1-x/2) ;
  end ;
end ;


% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end


%  the procedure uses cubic B-splines

%  Example:
%  xd = [ 1 2 3 4 5];  yd = [-1 2 -1 0 -1];
%  dyd = [-3 2];
%  x = 1:0.05:5;  y = fspline(xd,yd,dyd,x);
%  plot(xd,yd,'o',x,y,'r')
%  fspline(xd,yd,dyd,3.1)

%  It has been tested on MATLAB, version R2010b and version R2012a
%  version: 1.0
%  March 5, 2013















