% Engg. Maths for Adv. Studies - Autumn 2019 IIT Dharwad
% Module Numerical Methods Assignment 01
% Sample Solution using Octave - Question 1
% Author - Samarth Raut
% Last modified - 17th Nov. 2019
%**********************************************************

% Code for Cubic Spline demonstration using 8 datapoints
% dataTable=[ SrNo | x | y] 

% Assumption - 8 equispaced datapoints

clear;
% Sinc function
dataTable = ...
[0  0.0     0.12507;
 1  2.0    -0.12828;
 2  4.0    -0.10022;
 3  6.0     0.66500;
 4  8.0     0.66500;
 5 10.0    -0.10022;
 6 12.0    -0.12828;
 7 14.0     0.12507]

% sin function
%dataTable = ...
%[ 1  0.00000   0.00000
%  2 0.44880   0.43388
%  3 0.89760   0.78183
%  4 1.34640   0.97493
%  5 1.79520   0.97493
%  6 2.24399   0.78183
%  7 2.69279   0.43388
%  8 3.14159   0.00000]
 
 t=dataTable(:,2);
 v=dataTable(:,3);
 
 h=t(2)-t(1)  % Assumption - equispaced points on x axis
 
 A=[1.0     0     0     0     0     0     0     0;
    h/6   2*h/3 h/6   0     0     0     0     0;
    0     h/6   2*h/3 h/6   0     0     0     0;
    0     0     h/6   2*h/3 h/6   0     0     0;
    0     0     0     h/6   2*h/3 h/6   0     0;
    0     0     0     0     h/6   2*h/3 h/6   0;
    0     0     0     0     0     h/6   2*h/3 h/6;
    0     0     0     0     0     0     0     1.0; ]
  
f=zeros(size(v,1),1);
%Free Run out conditions
f(1,1)=0;
f(size(v,1),1)=0;

for i=2:(size(v,1)-1)
  f(i)=(v(i-1)-2*v(i)+v(i+1))/h;
end
f
g2t=A\f

%Evaluation of the function value using Cubic Spline for given x
j=1;
span=t(size(t,1),1)-t(1);
nFrag=30;
for k=1:nFrag;
  x=t(1) + k/nFrag*span;

  i=floor((x-t(1))/h)+1;
  if i>=size(g2t,1)
    i=size(g2t,1)-1;
  end
  
y = g2t(i)/6*( (t(i+1)-x)^3/h - h*(t(i+1)-x) ) ...
    + g2t(i+1)/6*( (x - t(i))^3/h - h*(x-t(i)) ) ...
    + v(i)/h*(t(i+1)-x) + v(i+1)/h*(x-t(i));

X1(j,1)=x;
Y1(j,1)=y;
j=j+1;
end

[X1, Y1];
%%plot(t,v, 's', X1,Y1, '*', X1, sin(X1-7.0)./(X1-7.0) )
%% Analytic does not overlap  as input data for Assignment 01 was scaled in addition to translation...Oops!:)
plot(t,v, 's', X1,Y1, '*--')
legend("input data", "Cubic Spline approx", " Analytic Sinc Fn")
%plot(t,v, 'sf', X1,Y1, '*', X1, sin(X1) )
%legend("input data", "Cubic Spline approx", " Analytic Sin Fn")
xlabel("time")
ylabel("vel")