function [a,b,c,d,sp]=nma_spline(x,y)
% computes the cubic splines between any number of points.

%function [a,b,c,d]=nma_spline(x,y)
%
% computes the cubic splines between any number of points.
% a x^3 + b x^2 + c x + d
%
% INPUT:
%   x: an array of the x-coordinates of the points.
%   y: an array of the y-coordinates of the points.
%
% OUTPUT:
%   a: an array of the 'a' coefficients of the splines.
%   b: an array of the 'a' coefficients of the splines.
%   c: an array of the 'a' coefficients of the splines.
%   d: an array of the 'a' coefficients of the splines.
%   sp: a cell array. each cell contains a string that
%       represents a spline.
%
% EXAMPLE:
%
% x =[ 1     2     3];
% y =[ 3     4     5];
% [a,b,c,d,sp]=nma_spline(x,y);
% sp{1}
% ans =
% +0.000000*(x-1.000000)^3 +0.000000*(x-1.000000)^2 +1.000000*(x-1.000000) +3.000000

% copyright: Nasser M. Abbasi

n=length(x);
interval = zeros(n-1,1);
for i=1:n-1
    interval(i)=x(i+1)-x(i);
end

j=0;
C = zeros(n-2,1);
for i=2:n-1
    j=j+1;
    C(j)=6*( (y(i+1)-y(i))/interval(i)- (y(i)-y(i-1))/interval(i-1));
end

j=0;
A=zeros(n-2,n-2);
for i=2:n-1
    j=j+1;
    A(j,j)=2*( interval(j) + interval(j+1) );
    if(i<n-1)
        A(j+1,j)=interval(j+1);
        A(j,j+1)=interval(j+1);
    end
end

%now solve As=b
S1=A\C';

S(1)=0;
S(2:n-1)=S1;
S(n)=0;

a = zeros(n-1,1);
b = zeros(n-1,1);
c = zeros(n-1,1);
d = zeros(n-1,1);
sp = cell(n-1,1);

for i=1:n-1
    a(i)= (S(i+1)-S(i))/(6*interval(i));
    b(i)= S(i)/2;
    c(i)= (y(i+1)-y(i))/interval(i)-( 2*interval(i)*S(i) + interval(i)*S(i+1) )/6;
    d(i)= y(i);
    
    sp{i}=sprintf('%+f*(x-%f)^3 %+f*(x-%f)^2 %+f*(x-%f) %+f',a(i),x(i),b(i),x(i),c(i),x(i),d(i));
end