%------------------------------------------------------------------------
% NAME:     qpcls_invchar
%
%           Performs a characterisation of the general retrieval performance
%           assuming linearity for CLS inversions.
%
%           If the function has output arguments, the main characterisation
%           quantities are returned. If there are no return arguments,
%           figures are produced by QP_INVCHAR.
%
% FORMAT:   qpcls_invchar( Q, mr_limit )
%
%              or
%
%           [kx_names,kx_index,kx_aux,PTZ,measres,s_tot,s_obs,A,C] =
%                                                qpcls_invchar( Q ,mr_limit )
%
% OUT:      kx_names   As the corresponding ARTS variable.
%           kx_index   Start and stop index for each retrieval identity.
%                      A 2 column matrix.
%           kx_aux     As the corresponding ARTS variable.
%           PTZ        Matrix with pressure, temperature and altitude.
%           measres    Measurement responses as a vector.
%           s_tot      Total retrieval error (sum of smoothing and observation
%                      errors).
%           s_obs      Observation error.
%           vertres    Vertical resolutions.
%           A          Averaging kernel matrix.
%           C          Error correlation matrix.
% IN:       Q          Setting structure.
%
% OPTIONAL  mr_limit   Do not include altitudes with a measurement response
%                      below this limit (only valid for plots). Default is 0.3.
%                      To set MR_LIMIT=0, can result in crash.
%------------------------------------------------------------------------

% HISTORY: 2001.10.14  Created by Patrick Eriksson


function [kx_names,kx_index,kx_aux,PTZ,measres,s_tot,s_obs,A,C] = qpcls_invchar( Q, mr_limit )

%=== default value
if ~exist('mr_limit')
  mr_limit = 0.3;
end

%=== Create temporary dir
%
tmpdir    = temporary_directory( Q.TMP_AREA );

out(1,1);
out(1,'ERROR CHARACTERISATION: CONSTRAINED LEAST SQUARES');
out(1,0);
out(1,'Calculating Kx, Kb...');

%=== Load H matrices
%
load( [Q.OUT,'.h'], '-mat' );
load( [Q.OUT,'.hd'], '-mat' );

 
%=== Use QP_ITER1 to get Kx etc.
%
[ fy, Kx, kx_names, kx_index, kx_aux, kx_aux_abs,Kb, kb_names, cfile, basename] = ...
                                         qp_iter1( Q, H, Hd, tmpdir, 3, 0, 0 );
clear H Hd


out(1,'Calculating Dy ...');


%=== Load Seinv and Sxinv from Q.OUT
%
load( [Q.OUT,'.seinv'], '-mat' );
if strcmp( lower( Q.RETRIEVAL_METHOD ), 'oem' )
  load( [Q.OUT,'.sxinv'], '-mat' );
  RM = Sxinv;
  clear Sxinv
end
%
if length(fy) ~= size(Seinv,1)
  error('Length of F(y) and size of Se do not match (recalculation needed?).')
end 
%
if size(Kx,2) ~= size(RM,1)
  error('Sizes of Kx and regularization matrix do not match (recalculation needed?).')
end 



%=== Delete the temporary directory
delete_tmp_dir( Q.TMP_AREA, tmpdir );



%=== Calculate Dy
KtSeinv = Kx' * Seinv;
Dy      = inv( RM + KtSeinv*Kx ) * KtSeinv;


%=== Clear varaibles not needed
%
clear Seinv RM KtSeinv


%=== Use qp_invchar
%
[PTZ, measres,s_tot,s_obs,A,C] = ...
        qp_invchar( Q, Dy, Kx, kx_names, kx_index, kx_aux, ~nargout, mr_limit);

out(1,-1);
if out(1)
  fprintf('\n');
end

% save all the variables 
%strsave=[' save ', Q.OUT, '_retr_Char.mat Dy Kx kx_names kx_index kx_aux kx_aux_abs PTZ measres s_tot s_obs A C']; 
%eval(strsave);