function y = splineA(xd,yd,x,Ends,Ders)

% Nonuniform Cubic Spline Interpolation
% Includes any combination of first or second derivative end conditions
% If you want the formulas for the resulting cubics, use splineB.m

% Usage (examples at end of file)
% For a natural cubic spline (this is the default version)
%   y = splineA(xd,yd,x)
% If one or both ends involve a derivative condition, then
%   y = splineA(xd,yd,x,Ends,Ders)

% Output:
% If n = length of x, then y is a n-vector containing the values of the
% cubic spline at the respective x locations.

% Required arguments:
%   xd, yd: data (row) vectors with xd(i)<xd(i+1)
%   x: point or points (row vector) where y(x) is to be evaluated

% Optional arguments (only needed if a derivative end condition is used):
%   Ends = [ a  b ]: this specifies which end conditions to use
%       a = 1: the first derivative is specified at left end
%       a = 2: the second derivative is specified at left end
%       b = 1: the first derivative is specified at right end
%       b = 2: the second derivative is specified at right end
%   Ders = [ c  d ]: this specifies the derivative values
%       c = value of 1st, or 2nd, derivative at left end
%       d = value of 1st, or 2nd, derivative at right end
%       This row vector is required if Ends is included in calling command.

% check for errors
n=length(xd);
nd=length(yd);
if nargin < 2 || nargin==4
    error('You are missing one of the required arguments');
end
if n < 2
    error('You need at least three data points')
end
if n ~= nd
    error('xd and yd must have the same length')
end
if n ~= length(unique(xd))
    error('the xd values can not repeat')
end

% set the defaults
if nargin == 3
    Ends=[2 2];
    Ders=[0 0];
end

% the code begins
% extend nodes for B_0(x) and B_{n+1}(x)
h=diff(xd); h1=h(1); hn=h(n-1);
h=[h1 h1 h1 h hn hn hn];
Xd=[xd(1)-3*h1 xd(1)-2*h1 xd(1)-h1 xd xd(n)+hn xd(n)+2*hn xd(n)+3*hn];

% compute the coefficients for the B_i's
[B,G]=coeffs(h);

% interpolate data
if Ends(1)==2
    alpha0=G(4,1)/G(2,1);
    beta0=G(6,1)/G(2,1);
    gamma0=Ders(1)/G(2,1);
elseif  Ends(1)==1
    alpha0=G(3,1)/G(1,1);
    beta0=G(5,1)/G(1,1);
    gamma0=Ders(1)/G(1,1);
else
    error('Ends values are incorrect')
end
B0=Bi(xd(1),0,B,Xd);
Aa(1)=2/3-B0*alpha0;
Ac(1)=Bi(xd(1),2,B,Xd)-B0*beta0;
Ab(1)=0;
yd(1)=yd(1)-B0*gamma0;
for ii=2:n-1
    Aa(ii)=2/3;
    Ac(ii)=Bi(xd(ii),ii+1,B,Xd);
    Ab(ii)=Bi(xd(ii),ii-1,B,Xd);
end
if Ends(2)==2
    alphan=G(2,2)/G(6,2);
    betan=G(4,2)/G(6,2);
    gamman=Ders(1)/G(6,2);
elseif  Ends(2)==1
    alphan=G(1,2)/G(5,2);
    betan=G(3,2)/G(5,2);
    gamman=Ders(2)/G(5,2);
else
    error('End values are incorrect')
end
Bn1=Bi(xd(n),n+1,B,Xd);
Aa(n)=2/3-Bn1*betan;
Ac(n)=0;
Ab(n)=Bi(xd(n),n-1,B,Xd)-Bn1*alphan;
yd(n)=yd(n)-gamman*Bn1;
w=tri(Aa,Ab,Ac,yd);
a0=-alpha0*w(1)-beta0*w(2)+gamma0;
anp1=-alphan*w(n-1)-betan*w(n)+gamman;
w=[a0 w anp1];

% evaluate interpolation function
nx=length(x);
y=zeros(1,nx);
for ix=1:nx
    sum=0;
    for k=1:n+2
        sum=sum+w(k)*Bi(x(ix),k-1,B,Xd);
    end
    y(ix)=sum;
end


% local functions used by algorithm
% 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

function [B,G] =coeffs(h)
% compute the values of the coefficients of the the B_i's
% B = (n+2)x16 [16=(4 subintervals for B_i)*(4 coeffs a,b,c,d)]
% G = 6x2
%   1st col=B_0'(x1),B0''(x1),B_1'(x1),B1''(x1),B_2'(x1),B2''(x1)
%   2nd col=B_n-1'(xn),B_n-1''(xn),B_n-1'(xn),B_n''(xn),B_n+1'(xn),B_n+1''(xn)
% N=n+5; number of B_i's=n+2
N=length(h);
B=zeros(N-3,16);
G=zeros(6,2);
% i = j-1
for j=1:N-3
    hm2=h(j); hm1=h(j+1); h1=h(j+2); h2=h(j+3);
    
    % find d_{i-2} and d_{i+2}
    d2=solver(hm2,hm1,h1,h2);
    
    % interval x_{i-2}<x<x_{i-1}
    B(j,4)=d2(1);
    
    % interval x_{i-1}<x<x_{i}
    B(j,5)=d2(1)*hm2^3;
    B(j,6)=3*d2(1)*hm2^2;
    B(j,7)=3*d2(1)*hm2;
    Tm3=(hm2+hm1)^3-hm1^3;
    B(j,8)=(2/3-Tm3*d2(1))/hm1^3;
    
    % interval x_{i}<x<x_{i+1}
    B(j,9)=-d2(2)*h2^3;
    B(j,10)=3*d2(2)*h2^2;
    B(j,11)=-3*d2(2)*h2;
    T3=(h2+h1)^3-h1^3;
    B(j,12)=(-2/3-T3*d2(2))/h1^3;
    
    % interval x_{i+1}<x<x_{i+2}
    B(j,16)=d2(2);
    
    % B_0
    if j==1
        G(1,1)=3*d2(2)*h2^2;
        G(2,1)=-6*d2(2)*h2;
    end
    
    % B_1
    if j==2
        G(3,1)=3*d2(2)*h2*(h2+2*h1)+3*B(j,12)*h1^2;
        G(4,1)=-6*(d2(2)*h2+B(j,12)*h1);
    end
    
    % B_2
    if j==3
        G(5,1)=3*d2(1)*hm2^2;
        G(6,1)=6*d2(1)*hm2;
    end
    
    % B_{n-1}
    if j==N-5
        G(1,2)=3*d2(2)*h2^2;
        G(2,2)=-6*d2(2)*h2;
    end
    
    % B_n
    if j==N-4
        G(3,2)=3*d2(2)*h2*(h2+2*h1)+3*B(j,12)*h1^2;
        G(4,2)=-6*(d2(2)*h2+B(j,12)*h1);
    end
    
    % B_{n+1}
    if j==N-3
        G(5,2)=3*d2(1)*hm2^2;
        G(6,2)=6*d2(1)*hm2;
    end
    
end

% find d_{i-2} and d_{i+2}
function d2=solver(hm2,hm1,h1,h2)
a=h1*hm2*(hm2+hm1)^2;
b=-hm1*h2*(h2+h1)^2;
c=h1^2*hm2*(hm2+hm1)*(hm2+2*hm1);
d=hm1^2*h2*(h2+h1)*(h2+2*h1);
A=[[a b];[c d]];
bb=[2*(hm1+h1)/3; -2*(hm1^2-h1^2)/3];
d2=A\bb;

% evaluate B_i(x)
function g=Bi(x,i,B,Xd)
j=i+1;
if x<Xd(j) || x>Xd(j+4)
    g=0;
elseif x<Xd(j+1)
    g=B(j,4)*(x-Xd(j))^3;
elseif x<Xd(j+2)
    g=B(j,5)+B(j,6)*(x-Xd(j+1))+B(j,7)*(x-Xd(j+1))^2+B(j,8)*(x-Xd(j+1))^3;
elseif x<Xd(j+3)
    g=B(j,9)+B(j,10)*(x-Xd(j+3))+B(j,11)*(x-Xd(j+3))^2+B(j,12)*(x-Xd(j+3))^3;
else
    g=B(j,16)*(x-Xd(j+4))^3;
end


% % Examples
%
% % generate some random data
% n = 8;
% xd = rand(1,n)+0:n;
% yd = randi(n,1,n);
% x = linspace(xd(1),xd(n),100);
%
% % natural spline with plot
% y = splineA(xd,yd,x);
% plot(x,y,'r','LineWidth',1.6)
% hold on
% plot(xd,yd,'ob','MarkerSize',8)
%
% % clamped spline
% y = splineA(xd,yd,x,[1 1],[0.2 -1]);
%
% % combination examples
% y = splineA(xd,yd,x,[1 2],[-0.2 0.6]);
% y = splineA(xd,yd,x,[2 1],[-0.2 0.6]);
%
% % check on a particular value
% y0 = splineA(xd,yd,pi)
% y0 = splineA(xd,yd,7.1,[2 1],[-0.2 0.6])

% Comments:
% 1) The nodes xd do not need to be uniformly spaced, but it is required
% that they are sequential, so xd(i)<xd(i+1).
% 2) This program replaces, and extends, two earlier spline programs that I
% wrote: nspline.m and fspline.m
% 3) The code is reasonably fast, and as example, with length(xd)=200
% and length(x)=1000, the time is about 8 msec.
% 4) The algorithm uses a B-spline construction, and it is described in the
% text "Introduction to Scientific Computing and Data Analysis" (Holmes, 2016)

%  Versions:
%  version 1.0: May 3, 2018
%  version 1.01: May 6, 2018 (made a few minor changes in the comments; the 
%  code was not changed)



















