% Lab 3 Divided Difference Interpolation
clear; clc;
theta=(pi/180)*[0 10 20 30 40 50 60 70 80 90];
% theta=[0 10 20 30 40 50 60 70 80 90];
ft=[0 17.3 34.2 50 64.3 76.6 86.6 93.9 98.4 100];
% Calculation of divided differences
for i=1:1:9
    fd1(i,1)=(ft(1,i)-ft(1,i+1))/(theta(i)-theta(i+1));
end
% 2nd differences
for i=1:1:8
    fd2(i,1)=(fd1(i,1)-fd1(i+1,1))/(theta(i)-theta(i+2));
end
% 3rd differences
for i=1:1:7
    fd3(i,1)=(fd2(i,1)-fd2(i+1,1))/(theta(i)-theta(i+3));
end
% 4th differences
for i=1:1:6
    fd4(i,1)=(fd3(i,1)-fd3(i+1,1))/(theta(i)-theta(i+4));
end
% 5th differences
for i=1:1:5
    fd5(i,1)=(fd4(i,1)-fd4(i+1,1))/(theta(i)-theta(i+5));
end
% 6th differences
for i=1:1:4
    fd6(i,1)=(fd5(i,1)-fd5(i+1,1))/(theta(i)-theta(i+6));
end
% 7th differences
for i=1:1:3
    fd7(i,1)=(fd6(i,1)-fd6(i+1,1))/(theta(i)-theta(i+7));
end
% 8th differences
for i=1:1:2
    fd8(i,1)=(fd7(i,1)-fd7(i+1,1))/(theta(i)-theta(i+8));
end
% 9th differences
for i=1:1:1
    fd9(i,1)=(fd8(i,1)-fd8(i+1,1))/(theta(i)-theta(i+9));
end
%------------------------
x=0:0.1:pi/2;
f=ft(1,1)+(x-theta(1))*fd1(1,1)+(x-theta(1)).*(x-theta(2))*fd2(1,1)+ ...
    (x-theta(1)).*(x-theta(2)).*(x-theta(3))*fd3(1,1)+ ...
    (x-theta(1)).*(x-theta(2)).*(x-theta(3)).*(x-theta(4))*fd4(1,1)+ ...
    (x-theta(1)).*(x-theta(2)).*(x-theta(3)).*(x-theta(4)).*(x-theta(5))*fd5(1,1)+ ...
    (x-theta(1)).*(x-theta(2)).*(x-theta(3)).*(x-theta(4)).*(x-theta(5)).*(x-theta(6))*fd6(1,1)+ ...
    (x-theta(1)).*(x-theta(2)).*(x-theta(3)).*(x-theta(4)).*(x-theta(5)).*(x-theta(6)).*(x-theta(7))*fd7(1,1)+ ...
    (x-theta(1)).*(x-theta(2)).*(x-theta(3)).*(x-theta(4)).*(x-theta(5)).*(x-theta(6)).*(x-theta(7)).*(x-theta(8))*fd8(1,1)+ ...
    (x-theta(1)).*(x-theta(2)).*(x-theta(3)).*(x-theta(4)).*(x-theta(5)).*(x-theta(6)).*(x-theta(7)).*(x-theta(8)).*(x-theta(9))*fd9(1,1);
polar(x,f); grid on; hold on
polar(theta,ft,'X');