% load grid
[h hm]=loadHeightGridVectors('c:\McArtim\grids\height_grid.txt');

% a priori state vector and covariance matrix
p=loadLIProfile('c:\McArtim\medium\NO2_prf.txt', hm);
n=20;		% use only first n components of the state vector
xa=n(1:n);
[Sa Sai]=buildDiagonalCovarianceMatrix(0.5*xa);	% allow 50% error
alpha=1;		% use a priori
	
% regularisation matrix
R=buildRegularisationMatrix(h(1:n));
beta=1e-18;	% reg. parameter

% measurement vector and covariance matrix
Y=load('c:\measurement\NO2\DSCDs.txt');
y=Y(:,1);	% assume first column to have the SCDs
ye=Y(:,2);	% second column: SCDs errors
[Se Sei]=buildDiagonalCovarianceMatrix(ye);

% load K matrix from previous McArtim run
[M Me]=loadMcArtimResultMatrix('c:\McArtim\BoxSensitivities.txt', 4, n);
[K Ke]=buildDifferentialJacobiMatrix(M, Me);

% conversion km -> cm
K*=1e5;
Ke*=1e5;

% inversion
B=K'*Sei;
C=inverse(B*K+alpha*Sai+beta*R);	% <- this is the algebraic inversion

% state retrieval
x=C*(B*y+alpha*Sai*xa);
xe=C*(B*ye+alpha*Sai*xa);

% averaging kernel
AK=C*B*K;
writeMatrixToFile (AK, 'AK.txt');

% construct a nice result file
Profile=[hm, x, xe, diag(AK)];
writeMatrixToFile (Profile, 'profile.txt');