%------------------------------------------------------------------------
% NAME:     qpoem_invchar
%
%           Performs a characterisation of the general retrieval performance
%           assuming linearity.
%
%           If the function has output arguments, the main characterisation
%           quantities are returned. If there are no return arguments,
%           figures with the following plots are produced:
%
%           Measurement response 
%             The measurement response is the sum of the averaging kernels.
%             The sum is calculated only for the retrieval identity
%             sub-matrix. 
%
%           Total and observation error
%    
%           The FWHM of the averaging kernels
%             FWHM = full width at half maximum
%
%           The averaging kernels
%
%           Correlation FWHM
%             The distance over which the error correlation is >= 0.5.
%
%           Max/min correlation to other retrieval identities.
%             Correlation to all other altitudes etc. is considered.
%
% FORMAT:   qpoem_invchar( Q [,mr_limit] )
%
%              or
%
%           [kx_names,kx_index,kx_aux,PTZ,measres,s_tot,s_obs,vertres,A,C] =
%                                                qpoem_invchar( Q [,mr_limit] )
%
% OUT:      -
% IN:       Q          Setting structure.
% OPTIONAL: mr_limit   Do not include altitudes with a measurement response
%                      below this limit. Default is 0.3.
%                      To set MR_LIMIT=0, can result in crash.
%------------------------------------------------------------------------

% HISTORY: 2001.08.23  Created by Patrick Eriksson


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


if ~exist('mr_limit'),   mr_limit = 0.3; end



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


%=== Run ARTS to get Kx etc.
%
[fy,H,Hd,Kx,kx_names,kx_index,kx_aux] = qpoem_iter1( Q, tmpdir, 0 );
clear H Hd
%
np  = size( Kx, 2 );
nrq = size( kx_names, 1 );



%=== Load some other quantities from Q.OUT
%
load( [Q.OUT,'.sxinv'], '-mat' );
load( [Q.OUT,'.seinv'], '-mat' );
%
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(Sxinv,1)
  error('Sizes of Kx and Sx do not match (recalculation needed?).')
end 



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



%=== Calculate main quantities
KtSeinv = Kx' * Seinv;
Dy      = inv( Sxinv + KtSeinv*Kx ) * KtSeinv;
A	= Dy * Kx;
AI	= A - eye(size(A));

%
clear KtSeinv Sxinv Seinv
%
load( [Q.OUT,'.se'], '-mat' );
load( [Q.OUT,'.sx'], '-mat' );
%
C       = Dy*Se*Dy';
s_obs	= sqrt( diag( C ));
C       = C + AI*Sx*AI';
s_tot   = sqrt( diag( C ));
C       = C ./ (s_tot*s_tot');   
%
clear AI Sx Se



%=== Calculate measurement response
%
z       = zeros(np,1);
measres = zeros(np,1);
vertres = zeros(np,1);
%
for i = 1:nrq
  ind = kx_index(i,1) : kx_index(i,2);
  measres(ind) = sum( A(ind,ind)' )';
end


%=== Read ptz
%
PTZ = read_datafile( Q.APRIORI_PTZ, 'MATRIX', 1 );



%=== Return if output arguments ===           %=== 
%                                             %===  
if nargout                                    %=== Return ?
  return                                      %===
end                                           %=== 



%=== Set general axes settings
%
A1 = axes_def; A2 = axes_def; A3 = axes_def;
A4 = axes_def; A5 = axes_def; A6 = axes_def;
%
A1.POSITION = [0.08 0.68 0.30 0.26];
A2.POSITION = [0.40 0.68 0.58 0.26];
A3.POSITION = [0.08 0.37 0.30 0.26];
A4.POSITION = [0.40 0.37 0.58 0.26];
A5.POSITION = [0.08 0.06 0.30 0.26];
A6.POSITION = [0.40 0.06 0.58 0.26];
%
A1.XGRID    = 'on';
A2.XGRID    = 'on';
A3.XGRID    = 'on';
A5.XGRID    = 'on';
A6.XGRID    = 'on';
%
A2.YTICKLABEL = [];
A4.YTICKLABEL = [];
A6.YTICKLABEL = [];



%=== First figure - 1
if isempty(get(0,'Children'))     % To handle the case when no window
  fig = 0;                        % is opened
else
  fig = gcf - 1;
end



%=== Loop the retrieval identities ============================================
%
for i = 1:nrq


  %= Names of present retrieval identity
  [group,name] = qp_kinfo( kx_names{i} );
  %= Keep just main name
  idash   = min(find( name == '-' ));
  if ~isempty( idash )
    name    = name( 1: idash-1 ); 
  end

  %= Index in x of present retrieval identity
  ind = kx_index(i,1):kx_index(i,2);

  %= Get index for points with measurement response above limit
  ip = find( measres(ind) >= mr_limit );


  %= Do something only if there are points > mr_limit
  if isempty(ip)
    fprintf('\nRetrieval identity: %s\n', name );
    fprintf('----------------------------------\n');
    fprintf('All data below limit on measurement response\n\n');
  else

    %= Create figure
    fig = fig + 1;
    figure(fig);
    fig_size(18,24,'cm',0.8);


    %=== Scalar identity ======================================================
    if length( ind ) == 1

      row    = 0.016;
      fsize  = 10;
      y0     = 0.97;
      xoff   = 0.02;
      textw  = 0.25;

      set(gcf,'Unit','Cent');
      set(gca,'Position',[-0.01 -0.01 1 1]);
      set(gca,'Box','Off');
      set(gca,'TickLength',[0 0])
    
      irow = 1;
      text2( xoff, y0-irow*row, sprintf('%s',name),fsize*1.4,1);
  
      irow = irow+2;
      text2(xoff,y0-irow*row,'Measurement response',fsize,1);
      text2(xoff+textw,y0-irow*row,sprintf(': %.2f',measres(ind)),fsize,0);

      irow = irow+2;
      text2(xoff,y0-irow*row,'Total error',fsize,1);
      text2(xoff+textw,y0-irow*row,sprintf(': %.2e',s_tot(ind)),fsize,0);
      irow = irow+1;
      text2(xoff,y0-irow*row,'Observation error',fsize,1);
      text2(xoff+textw,y0-irow*row,sprintf(': %.2e',s_obs(ind)),fsize,0);

      irow = irow+2;
      text2(xoff,y0-irow*row,'Correlation with',fsize,1);
      irow = irow+2;
      for j = 1:nrq
        if j~=i
          nr    = j-(j>i);
          ind2 = kx_index(j,1):kx_index(j,2);
          [gg,nn] = qp_kinfo( kx_names{j} ); 
          if mod(nr,2)== 1
            x = xoff;
          else
            x = 0.5;
          end
          text2(x,y0-irow*row,sprintf('%s',nn),fsize,0);
          if length(ind2) > 2
            text2(x+textw,y0-irow*row,...
                          sprintf(': %.2f <-> %.2f',min(C(ind2,ind)),...
                                                    max(C(ind2,ind))),fsize,0);
          else
            text2(x+textw,y0-irow*row,sprintf(': %.2f',C(ind2,ind)),fsize,0);  
          end
          if mod(nr,2)== 0
            irow = irow + 1;
          end
        end
      end


    %=== Vector identity ======================================================
    else

      ip = ind(min(ip)):ind(max(ip));

      %= Some plotting quantities
      %
      if strcmp( group, 'Species' )
        %
        z(ind)       = interpp( PTZ(:,1), PTZ(:,3), kx_aux(ind,1) );
        %
        yfac = 1000; ytext = 'Altitude [km]'; yunit = 'km';
        xfac = 100;  xunit = '%';
        %
      elseif strcmp( group, 'Temperature' )
	%
        z(ind)       = interpp( PTZ(:,1), PTZ(:,3), kx_aux(ind,1) );
        %
        yfac = 1000; ytext = 'Altitude [km]'; yunit = 'km';
        xfac = 1;    xunit = 'K';
        %    
      elseif strcmp( group, 'Continuum' )
	%
        z(ind)       = interpp( PTZ(:,1), PTZ(:,3), kx_aux(ind,1) );
        %
        yfac = 1000; ytext = 'Altitude [km]'; yunit = 'km';
        xfac = 1;    xunit = '1/m';
        %    
      else
	error(sprintf( 'Unknown retrieval identity: %s', name ));
      end


      %=== Figure 1 ===
      %
      A1.XLABEL   = 'Measurement response [%]';
      A1.YLABEL   = ytext;
      subplot(3,2,1)
      plot( measres(ip)*100, z(ip)/yfac );
      axis([mr_limit*100 119 min(z(ip))/yfac max(z(ip))/yfac]);
      ax_set( A1 );
      

      %=== Figure 2 ===
      %
      A2.XLABEL   = ['Error [',xunit,']'];
      subplot(3,2,2)
      pline( s_tot(ip)*xfac, z(ip)/yfac, 1 ); 
      hold on
      pline( s_obs(ip)*xfac, z(ip)/yfac, 2 );
      axis([0 max(s_tot(ip)*xfac)*1.2 min(z(ip))/yfac max(z(ip))/yfac])
      sleg = str2mat('Total','Observation');
      plineleg(sleg,0.65,0.7,0.12,10,0.07);
      ax_set( A2 );
      

      %=== Figure 3 ===
      %
      vertres(ind) = fwhm( z(ind), A(ind,ind) );
      %
      A3.XLABEL   = ['FWHM [',yunit,']'];
      A3.YLABEL   = ytext;
      subplot(3,2,3)
      plot( vertres(ip)/yfac, z(ip)/yfac )
      axis([0 max(vertres(ip)/yfac)*1.2 min(z(ip))/yfac max(z(ip))/yfac])
      ax_set( A3 );
      

      %=== Figure 4 ===
      %
      A4.XLABEL   = 'Averaging kernels [-]';
      subplot(3,2,4)
      plot( A(ip,ip), z(ip)/yfac )
      axis([-0.15 1.05 min(z(ip))/yfac max(z(ip))/yfac]);
      hold on
      for j = ip
        plot( [-2 2], [z(j) z(j)]/yfac, ':' )
      end
      hold off
      ax_set( A4 );
      

      %=== Figure 5 ===
      %
      A5.XLABEL   = ['Correlation FWHM [',yunit,']'];
      A5.YLABEL   = ytext;
      subplot(3,2,5)
      cl = fwhm( z(ind), abs(C(ind,ind)) );
      plot( cl/yfac, z(ind)/yfac )
      axis([0 max(cl/yfac)*1.2 min(z(ip))/yfac max(z(ip))/yfac])
      ax_set( A5 );
      

      %=== Figure 6 ===
      %
      A6.XLABEL   = 'Max/min correlation [-]';
      subplot(3,2,6)
      dtext = (max(z(ip))-min(z(ip)))/yfac/10;
      sleg = [];
      for j = 1:nrq
        if j~=i
          nr    = j-(j>i);
          ind2 = kx_index(j,1):kx_index(j,2);
          if length(ind2) > 2
            pline( max(C(ind2,ind)), z(ind)/yfac, nr );hold on
            pline( min(C(ind2,ind)), z(ind)/yfac, nr );
          else
            pline( C(ind2,ind), z(ind)/yfac, nr );hold on
          end
          [gg,nn] = qp_kinfo( kx_names{j} ); 
          %= Keep just main name
          idash   = min(find(nn=='-'));
          if ~isempty( idash )
            nn    = nn( 1: idash-1 ); 
          end
          sleg    = strvcat( sleg, nn );
        end
      end
      axis([-1 1 min(z(ip))/yfac max(z(ip))/yfac]);
      plineleg(sleg,0.05,0.93,0.12,10,0.06);
      hold off
      ax_set( A6 );
      
      
      suptitle(name,0.97,14);

    end
  end
end






%==============================================================================
%===
%=== Sub-functions
%===
%==============================================================================

function A = axes_def() 
  %
  A.BOX           = 'on';
  %
  A.FONTSIZE      = 10;
  A.FONTWEIGHT    = 'normal';
  %
  A.XLABEL_SIZE   = 10;
  A.XLABEL_WEIGHT = 'bold';
  %
  A.YLABEL_SIZE   = 10;
  A.YLABEL_WEIGHT = 'bold';
  %
return