%
%
  clear all
%

%
%  Set up the grid and parameter values.  Construct an error with sigma = 2.2.
%
  P = 0:.005:0.8;
  alpha_1 = -389.4;
  alpha_11 = 761.3;
  alpha_111 = 61.5;
  p = 3;
  n = length(P);

  sigma = 2.2;
  var = sigma^2;
  error = sigma*randn(size(P));

%
% Construct the solution, observations, and residuals.
%

  psi = alpha_1*P.^2 + alpha_11*P.^4 + alpha_111*P.^6;
  obs = psi + error;
  R = psi - obs;

%
% Construct the sensitivity matrix, observation error estimator, and covariance estimator.
%

  psi_alpha_1 = P.^2;
  psi_alpha_11 = P.^4;
  psi_alpha_111 = P.^6;

  sigma_est = (1/(n-p))*R*R'

  sensmat = [psi_alpha_1; psi_alpha_11; psi_alpha_111];
  V = sigma_est*inv(sensmat*sensmat');
  rank(V);

%
% Plot the data and residuals.
%

  figure(1)
  plot(P,psi,P,obs,'x','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('Polarization P')
  ylabel('Helmholtz Energy \psi(P)')

  figure(2)
  plot(P,psi - obs,'x',P,0*P,'k','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('Polarization P')
  ylabel('Residuals')