%
% File written by Annalinda Arroyo on 10/20/05
% Modified 10/25/05 by Peter Blomgren
% Modified 10/08/08 by Joe Mahaffy
%
% To call/use this function type: splineinterpolate

set(0,'defaultlinelinewidth',2)
set(0,'defaultaxesfontsize',16)
set(0,'defaultaxesfontweight','bold')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Given a set of points (saved under the name duckdata) this program
% will generate the coefficients of the spline that passes through the
% nodes (x) (to find one of the coefficients requires the trisolve
% function taken from the class web cite). Finally, the program will
% evaluate the spline at arbitary pionts (y) (this is done using the
% evalspline function also taken from the class web site), and then
% graphs both the original points and the polynomial that follows
% between them allowing you to see the fit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
load lynxdata
f = yd';                       % reset f as a column vector
x = xd';                       % reset x as a column vector
n = length(x);                % set n equal to the number of x-values

a = f;                        % set a equal to the column vector 'f'

h = diff(x);

% Matrix 'T' is the compact form of an nxn matrix 'A' where T(:,1) is
% the sub-diag of A, T(:,2) is the diag of A, and T(:,3) is the
% super-diag of A

T      = zeros(n,3);
T(1,2) = 1;                   % set both entries T(1,2)=T(n,2)=1
T(n,2) = 1; 

rhs    = zeros(n,1);

% Build matries 'h', 'T', and 'rhs' using the column matrix of the
% x-values %
T(2:(n-1),1) = h(1:length(h)-1);
T(2:(n-1),3) = h(2:length(h));
T(2:(n-1),2) = 2 * ( T(2:(n-1),1) + T(2:(n-1),3) );

da = diff(a);

rhs(2:(n-1)) = 3 ./ h(2:(n-1)) .* da(2:(n-1)) - 3 ./ h(1:(n-2)) .* da(1:(n-2));

c = trisolve(T,rhs);          % function that solves for c where T*c=rsh

% Build matries 'b' and 'd' using 'h' from pervious for loop, and
% 'c' from the trisolve function %
for i = 2:n
    b(i-1,1) = (1/h(i-1,1))*(a(i,1)-a(i-1,1)) - ...
	(h(i-1,1)/3)*(2*c(i-1,1)+c(i,1));
    d(i-1,1) = (c(i,1)-c(i-1,1)) / (3*h(i-1,1));
end

y = min(x):0.01:max(x);       % set y to arbitary pts to evaluate
                              % the spline at 
y = y';                       % reset y as a column vector 

S = evalspline(a,b,c,d,x,y);  % function that uses the coeff
                              % defining the spline to evaluate the
                              % spline at each pt in 'y' 

% graph pts (x,f) and the spline that fits these pts

plot(y,S,'b-',x,f,'mo','markersize',7)

grid on