%------------------------------------------------------------------------
% NAME:     qpcls_invchar
%
%           Calculates or plots the error due to each error source with
%           DO-level 1 or 2.
%
%           If the function has output arguments, the errors are returned. 
%           If there are no return arguments, figures are produced by 
%           QP_ERRORPLOT.
%
% FORMAT:   qpcls_error( Q, mr_limit )
%
%              or
%
%           [PTZ,E,Et] = qpcls_errorplot( Q ,mr_limit )
%
% OUT:      PTZ        Matrix with pressure, temperature and altitude.
%           E          Error for each source.
%           Et         Name on the error sources as cell array.
% 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 [PTZ,E,Et] = qpcls_errorplot( 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,'PLOTTING INDIVIDUAL ERRORS: CONSTRAINED LEAST SQUARES');
out(1,0);
out(1,'Calculating Kx ...');

%=== 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 ] = ...
                                         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_errorplot
%
[PTZ,E,Et] = ...
      qp_errorplot( Q, Dy, Kx, kx_names, kx_index, kx_aux, ~nargout, mr_limit);



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