% Unit 4 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  polynomial interpolation using un-evenly distributed data
%   plus fit of combined data sample to non-linear function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, close all, clc

% the data points
%-----------------
x=[1.1     1.2     2.4     3.1     4.1     5.5     6.1     7.0     8.4     8.7];
f=[1.65    1.86    3.45    3.55    2.60    1.73    1.19    0.88    0.31    0.27];
N=length(x);

% 1. plot data points
% (0.5 point)
%---------------------
plot(x,f,'o','MarkerSize',10,'LineWidth',3)

% 2. calculate f(x) at mid-points using polynomial interpolation
% (5 points)
%----------------------------------------------------------------
il = 2:N-2;
ir = 3:N-1;
xmid = x(il)+(x(ir)-x(il))./2;
fmid = p4interpol(x,f,xmid);

% 3. plot new mid-point data
% (0.5 point)
%----------------------------
hold on
plot(xmid,fmid,'rs','MarkerSize',10,'LineWidth',3)

% 4. combined sample of data points
% (1 point)
%-----------------------------------
xcom = [x xmid];
fcom = [f fmid];

% 5. obtain the amplitude of best-fit Planck curve
% (3 points)
%--------------------------------------------------
% Planck curve (without the amplitude)
g=@(x)(x.^3./(exp(x)-1));

% generate the column(!) matrix X containing dfda at x()
dfda = [g(xcom)];           % this is a line vector
X    = transpose([dfda]);   % this now is a column matrix

% calculate the optimisation matrix
Y    = inv(transpose(X)*X)*transpose(X);

% multiply optimisation matrix with f() (as a column vector!)
C = Y*fcom'

% plot best-fit Planck curve
m=linspace(min(xcom),max(xcom),100);
plot(m,C.*g(m)), hold on

legend('data points','polynomial interpolation to mid-points','best-fit Planck curve to combined data points','Location','NorthEast')
xlabel('x','FontSize',20)
ylabel('f(x)','FontSize',20)

