% Class 6 Polynomial Interpolation
% Announcements
% Class 5 - video (you watch)
% PSET 2 due tonight at Midnight 
% Makeup office hours for 2 hr after class, 1150
% PSET 3 posted, tricky, start early
%
% Examples of Interpolation with Matlab
% Coded by NFR on 1.31.19
clc
close all
Runge = @(x) 1/(1+25*x^2);
fplot(Runge,[-1 1])
%
% Example 17.11 from the textbook
T = [200 250 300 350 400 450];
rho = [1.708 1.367 1.139 0.967 0.854 0.759];
plot(T,rho,'o')
xlabel('Temperature (K)')
ylabel('Density (kg/m^3')
% First order fit using 300 and 350K data points
coef_1 = polyfit(T(3:4),rho(3:4),1)
Tvec = linspace(200,450,100);
rhovec_1 = polyval(coef_1,Tvec);
hold on
plot(Tvec,rhovec_1,'k:')
% To answer the question of rho at 330K
rho_330K_1 = polyval(coef_1,330)
% Third order interpolation
coef_3 = polyfit(T(2:5),rho(2:5),3);
rhovec_3 = polyval(coef_3,Tvec);
plot(Tvec,rhovec_3,'-r')
rho_330K_3 = polyval(coef_3,330)
% Fifth order interpolation
coef_5 = polyfit(T,rho,5);
rhovec_5 = polyval(coef_5,Tvec);
plot(Tvec,rhovec_5,'--g')
rho_330K_5 = polyval(coef_5,330)
% Inverse Interpolation Problem 
% See the issue of centering/scaling
%
hold off
D = (1920:10:2000);
P = [106.46 123.08 132.12 152.27...
    180.67 205.05 227.23 249.46 281.42];
plot(D,P,'o')
% What is the population in 1963 (normal interplation question)
% NOT inverse
%c_pop = polyfit(D,P,8);
% How to get rid of the error, scale the date vector
Span = 2000-1920; 
Mid = 1960;
Dscal = (D-Mid)/Span;
c_pop = polyfit(Dscal,P,8);
% Now to find the answer w/ the scaled values
Pop_ans = polyval(c_pop,(1963-Mid)/Span)
hold on
plot(1963,Pop_ans,'+')
% Inverse interpolation problem
% What year did the population hit 225?
YearVec = linspace(1920,2000,10000);
Pop_vec = polyval(c_pop,(YearVec-Mid)/Span);
plot(YearVec,Pop_vec,'r.')
% Let's now use logic to find the corresponding year
Rank = abs(Pop_vec-225);
[~,ind] = min(Rank);
Year_225 = YearVec(ind)
plot(Year_225,225,'xk')