clear all;
close all;

x = 1 : 5; % weekdays (1 - for Monday, 2 - for Tuesday, 3 - for Wednesday, 4 - for Thursday, 5 - for Friday)
y = [7607.40,7544.60,7561.20,7526.40,7515.30]; % TSE 300 index
z = [1848.20,1836.50,1860.20,1843.40,1805.20]; % NASDAQ index

% this script computes cubic runout splines through given data points (x,y) and (x,z)
n = length(y) - 1; % n = 4 for five data points
h = x(2) - x(1); % h = 1 for given data
A = 4*diag(ones(n-1,1)) + diag(ones(n-2,1),1) + diag(ones(n-2,1),-1); % coefficient matrix
A(1,1) = 6; A(n-1,n-1) = 6; A(1,2) = 0; A(n-1,n-2) = 0;
by = h^(-2)*6*(y(1:n-1) - 2*y(2:n) + y(3:n+1)); % right-hand-side
bz = h^(-2)*6*(z(1:n-1) - 2*z(2:n) + z(3:n+1));
my = A\by'; my = [2*my(1)-my(2); my; 2*my(n-1)-my(n-2)]'; % vector of curvatures
mz = A\bz'; mz = [2*mz(1)-mz(2); mz; 2*mz(n-1)-mz(n-2)]';

a = y(1:n); b = (y(2:n+1)-y(1:n))/h - h*(my(2:n+1)+2*my(1:n))/6;
c = 0.5*my(1:n); d = (my(2:n+1)-my(1:n))/(6*h); % coefficients of the cubic spline interpolants
xInt = 0:0.01:5;
for j = 1 : length(xInt)
	if xInt(j) <= x(1)
        iInt(j) = 1;
    elseif xInt(j) ~= x(n+1)
        iInt(j) = sum(x <= xInt(j));
	else
		iInt(j) = n;
	end
end
yInt = a(iInt) + b(iInt).*(xInt-x(iInt)) + c(iInt).*(xInt - x(iInt)).^2 + d(iInt).*(xInt - x(iInt)).^3;
figure(1); plot(x,y,'b*',xInt,yInt,'g'); 

% middle of Wednesday
xx = 2.5; ii = 2; 
yy = a(ii) + b(ii)*(xx-x(ii)) + c(ii)*(xx - x(ii))^2 + d(ii)*(xx - x(ii))^3;

a = z(1:n); b = (z(2:n+1)-z(1:n))/h - h*(mz(2:n+1)+2*mz(1:n))/6;
c = 0.5*mz(1:n); d = (mz(2:n+1)-mz(1:n))/(6*h); % coefficients of the cubic spline interpolants
zInt = a(iInt) + b(iInt).*(xInt-x(iInt)) + c(iInt).*(xInt - x(iInt)).^2 + d(iInt).*(xInt - x(iInt)).^3;
figure(2); plot(x,z,'b*',xInt,zInt,'g'); 
zz = a(ii) + b(ii)*(xx-x(ii)) + c(ii)*(xx - x(ii))^2 + d(ii)*(xx - x(ii))^3;
fprintf('In the middle of Wednesday: TSE 300 = %10.2f; NASDAQ = %10.2f\n',yy,zz);

yyS = spline(x,y,xx); zzS = spline(x,z,xx);
fprintf('Indices with MATLAB "spline": TSE 300 = %10.2f; NASDAQ = %10.2f\n\n',yy,zz);
% notice that the values of MATLAB splines are the same as the cubic runout spline

ySpline = spline(x,y,xInt); 
yNorm = norm(ySpline-yInt);
zSpline = spline(x,z,xInt);
zNorm = norm(zSpline-zInt);
fprintf('Norms of the difference between MATLAB "spline" and cubic runout spline:\n TSE 300 = %20.15f; NASDAQ = %20.15f',yNorm,zNorm);