%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Fit parabola to observations.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Given: (1) p data points (xi_i,eta_i).
%
%         (2) Background parabola
%             eta = b0 + b1*xi + b2*xi^2
%
%  Required: best fit parabola
%            eta = a0 + a1*xi + a2*xi^2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Data
xi  = [ 1    2     3     4     5    ]';
eta = [ 4.08 4.95  5.01  5.97  6.04 ]';

stem(xi,eta,'ro')
axis([0,6,3,8])
hold on
pause

y = eta;

H = [   1   xi(1)  xi(1)^2 ;
        1   xi(2)  xi(2)^2 ;
        1   xi(3)  xi(3)^2 ;
        1   xi(4)  xi(4)^2 ;
        1   xi(5)  xi(5)^2   ] ;

%  Background field
b0 = 2.5; b1 = 1.5; b2 = -0.2;
xb = [  b0  
        b1  
        b2 ];

xipara = 0:0.1:6;
etapara = b0 + b1*xipara + b2*xipara.^2;
plot(xipara,etapara,'m')
pause

%% Variances.
sigmab = 1.0;
sigmao = 1.0;

%  Error covariance matrices
sigbsq = sigmab^2;
B = [  sigbsq     0      0    
          0    sigbsq    0
          0       0   sigbsq  ];

sigosq = sigmao^2;
R = [  sigosq     0      0      0      0 
          0    sigosq    0      0      0
          0       0   sigosq    0      0
          0       0      0   sigosq    0
          0       0      0      0   sigosq ];


%% Solve the OI equations.
W = B*H'*inv(R+H*B*H');
xa = xb + W*(y-H*xb);

a0 = xa(1);
a1 = xa(2);
a2 = xa(3);

xipara = 0:0.1:6;
etapara = a0 + a1*xipara + a2*xipara.^2;
plot(xipara,etapara,'b')
hold off

return
