% Quadratic interpolation demonstration
% DRB 06-Dec-2010

lam = [10 11 12];           % x-axis values
y = [6 10 9];               % y-axis values

% matlab solution using polyfit
p = polyfit(lam,y,2);       % quadratic fit

% matlab solution using matrix inverse
Q = [lam(1)^2 lam(1) 1 ; lam(2)^2 lam(2) 1 ; lam(3)^2 lam(3) 1];
Y = y(:);
q = inv(Q)*Y;

% efficient solution using translated coordinates
phi = lam-lam(1);
z = y - y(1);
cbar = 0;
bbar = (z(3)-phi(3)^2/phi(2)^2*z(2))/(phi(3) - phi(3)^2/phi(2));
abar = (z(2)-phi(2)*bbar)/phi(2)^2;
a = abar;
b = bbar-2*lam(1)*abar;
c = cbar +abar*lam(1)^2 - bbar*lam(1)+y(1);

% maximum
lamstar = -b/(2*a);

lamc = (min(lam)-1):0.01:(max(lam)+1);
plot(lam,y,'x',...
     lamc,p(1)*lamc.^2+p(2)*lamc+p(3),'--',...
     lamc,q(1)*lamc.^2+q(2)*lamc+q(3),'--',...
     lamc,a*lamc.^2+b*lamc+c,'--',...
     lamstar,a*lamstar.^2+b*lamstar+c,'rp','Linewidth',2);
grid on
xlabel('\lambda');
ylabel('y=|A(\lambda)|^2');