%------------------------------------------------------------------------
% NAME:     qp_Se
% 
%           Calculates the Se matrix.
%
%           If there are no output arguments, Se and its inverse are
%           stored to files.
%
%           If H and Hd not are given, those variables are loaded from
%           files. 
% 
% FORMAT:   Se = qp_Se( Q [, H, Hd, klevels ] )
%
% OUT:      -
% IN:       Q          Setting structure.
% OPTIONAL: H          The sensor/data-reduction transfer matrix.
%           Hd         The data-reduction transfer matrix.
%           klevels    Vector with the do-levels to consider.
%                      See further qp_iter1.m.
%                      Default is 2.
%------------------------------------------------------------------------

% HISTORY: 2001.03.28  Created by Patrick Eriksson.
%          2001.08.25  Modified by Carlos Jimenez
%          2001.10.13  Included variables beside thermal noise (PE)


function Se = qp_Se( Q, H, Hd, klevels )


out(1,1);
out(1,'Setting up Se.');


%=== Set input arguments if not given
%
if nargin < 2
  load( [Q.OUT,'.h'], '-mat' );
end
if nargin < 3
  load( [Q.OUT,'.hd'], '-mat' );
end
if nargin < 4
  klevels = 2;
end


%=== Read some vectors to get sizes
%
load( [Q.OUT,'.f_sensor'], '-mat');
load( [Q.OUT,'.za_sensor'], '-mat');
%
nf   = length( f_sensor );
nza  = length( za_sensor );



%=== Measurement thermal noise
%
if Q.MEASNOISE_DO == 2
  %
  out(2,'Doing measurement thermal noise.');
  %
  S = sFromFile( fullfile( Q.SENSOR_DIR, Q.MEASNOISE_COVMAT ), f_sensor );
  %
  [row,col,s] = find( S ); 
  %
  n = length( row );
  %
  row = ones(n,1)*(0:nf:(nza-1)*nf) + repmat(row,1,nza);
  col = ones(n,1)*(0:nf:(nza-1)*nf) + repmat(col,1,nza);
  %
  if isfield( Q, 'MEASNOISE_SCFAC') & ~isempty( Q.MEASNOISE_SCFAC )
    if length(Q.MEASNOISE_SCFAC) ~= nza
      error( ...
     'The length of MEASNOISE_SCFAC is not equal to number of zenith angles.');
    end
    w  = ones(nf,1) * vec2row(Q.MEASNOISE_SCFAC);
    w  = (w(:)).^2;
    Se = sparse( row(:), col(:), repmat(s,nza,1).* w );
  else
    Se = sparse( row(:), col(:), repmat(s,nza,1) );
  end
  %
else
  %
  Se = 0;
  %
end


%=== Calibration thermal noise
%
if Q.CALINOISE_DO == 2
  %
  out(2,'Doing calibration thermal noise.');
  %
  S = sFromFile( fullfile( Q.SENSOR_DIR, Q.CALINOISE_COVMAT ), f_sensor );
  %
  [row,col,s] = find( S ); 
  %
  n = length( row );
  %
  row = ones(n,1)*(0:nf:(nza-1)*nf) + repmat(row,1,nza);
  row = repmat( row(:), 1, nza );
  col = repmat( col, 1, nza ); 
  col = ones(n*nza,1)*(0:nf:(nza-1)*nf) + repmat(col(:),1,nza);
  %
  Se = Se + sparse( row(:), col(:), repmat(s,nza*nza,1) );
  %
end



%=== Include data reduction
%
if size(Se,1) > 1
  out(2,'Applying data reduction on noise.');
  Se = Hd * Se * Hd';
end



%=== Variables beside thermal noise
%
out(2,'Creating Sb:');
Sb = qp_Sxb( Q, klevels );
%
if ~isempty(Sb)
  %
  out(2,'Calculating Kb ...');
  %
  tmpdir = temporary_directory( Q.TMP_AREA );
  [fy, Kb, kb_names, kb_index, kb_aux, VMRs] = ...
                                   qp_iter1( Q, H, Hd, tmpdir, klevels, 0, 0 );
  delete_tmp_dir( Q.TMP_AREA, tmpdir );
  %
  Se = Se + Kb * Sb * Kb';
  %
else
  out(2,'   empty');
end




%=== Save if no output argument.
%    Both Se and the inverse are saved.
%
if ~nargout 
  %
  save( [Q.OUT,'.se'], 'Se' );
  %
  Seinv = Se \ speye(size(Se));
  save( [Q.OUT,'.seinv'], 'Seinv' );
  %
end


out(1,-1);