% Class 6 - Interpolation (polynomials)
% Coded by NFR on 9.12.2019
%
% Problem 17.11 as an example
close all
clc
T = (200:50:450);
D = [1.708 1.367 1.139 0.967 0.854 0.759];
% GOAL: Density at 330K
% Using a first order polynomial
% Two points are needed for 1st order
% Pick the two closest to the desired value
% Points 3 and 4
p = polyfit(T(3:4),D(3:4),1); %<-- n = order of polynomial
% p(1) = slope, p(2) = intercept
% Use polyval to find solution
Sol1 = polyval(p,330)
% Third order fit
% 4 points
p = polyfit(T(2:5),D(2:5),3);
Sol2 = polyval(p,330)
% 5th Order fit
p = polyfit(T,D,1);
Sol3 = polyval(p,330)
% Plot all the points
plot(T,D,':o',330,Sol1,'r+',330,Sol2,'bo',330,Sol3,'k+')
% One more example, to show wrinkle of 
% inverse interpolation
% (this will also show poorly conditioned error)
Y_vec = (1920:10:2000);
P_vec = [106.46 123.08 132.12 152.27 180.67...
    205.05 227.23 249.46 281.42];
plot(Y_vec,P_vec,'o')
xlabel('Year')
ylabel('Population (Millions)')
% Show plot of switching axes, 
% why you DON'T do this for inverse interpolation
figure
plot(P_vec,Y_vec,'o')
ylabel('Year')
xlabel('Population (Millions)')
% You can't use polyfit/polyval in reverse
% (usually) because the y axis data is not 
% at even intervals
%
% Regular interpolation
% What is the pop in 1936?
% Test an 8th order polynomial
p = polyfit(Y_vec(1:4),P_vec(1:4),3);
Sol = polyval(p,1936);
plot(Y_vec,P_vec,'o',1936,Sol,'r+')
xlabel('Year')
ylabel('Population')
% Inverse interpolation
% Conditional logic method to solve
% Can use root finding (later)
% In what year did the pop hit 190 M?
p = polyfit(Y_vec(4:7),P_vec(4:7),3);
x_vec = linspace(1950,1980,1000);
sol_vec = polyval(p,x_vec);
[~,ind] = min(abs(sol_vec-190));
Year = x_vec(ind);
fprintf('The year is %.2f\n',Year)