%------------------------------------------------------------------------
% NAME:     qp_Se
% 
%           Calculates and returns (no saving) the Se matrix.
% 
% FORMAT:   qp_Se( Q, H, Hd, klevels )
%
% OUT:      -
% IN:       Q          Setting structure.
%           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.
%------------------------------------------------------------------------

% 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 )


%=== Read some vectors to get sizes
%
f  = read_datafile( [Q.SENSOR_DIR,'/',Q.F_SENSOR ], 'VECTOR' );
za = read_datafile( [Q.SENSOR_DIR,'/',Q.ZA_SENSOR ], 'VECTOR' );
%
nf   = length( f);
nza  = length( za );



%=== Measurement thermal noise
%
if Q.MEASNOISE_DO == 2
  %
  out(2,'Doing measurement thermal noise.');
  %
  S = sFromFile( [Q.SENSOR_DIR,'/',Q.MEASNOISE_COVMAT], f );
  %
  [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);
  %
  Se = sparse( row(:), col(:), repmat(s,nza,1) );
  %
else
  %
  Se = sparse( nf*nza, nf*nza );
  %
end



%=== Calibration thermal noise
%
if Q.CALINOISE_DO == 2
  %
  out(2,'Doing calibration thermal noise.');
  %
  S = sFromFile( [Q.SENSOR_DIR,'/',Q.CALINOISE_COVMAT], f );
  %
  [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
%
out(2,'Applying data reduction on noise.');
%
Se = Hd * Se * Hd';



%=== 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] = ...
                                   qp_iter1( Q, H, Hd, tmpdir, klevels, 0, 0 );
  delete_tmp_dir( Q.TMP_AREA, tmpdir );
  %
  Se = Se + Kb * Sb * Kb';
  %
  clear fy Kb 
else
  out(2,'   empty');
end
%
clear Sb


return



%------------------------------------------------------------------------------

%=== This is an old part
if 0

if Q.CALINOISE_DO == 2
  %
  out(2,'Doing calibration thermal noise.');
  %
  Noise  = read_datafile( [Q.SENSOR_DIR,'/',Q.CALINOISE_FILE], 'MATRIX' );
  sigma  = interp1( Noise(:,1), Noise(:,2), f_sensor );
  e      = sigma .* sigma;
  for i = 1:ncorrs
    covar  = ( sigma(1:(nf-i)) .* sigma((1+i):nf) ) * Q.CHCORR_DATA(i);
    eright = [ zeros(i,1); covar ]; 
    eleft  = [ covar; zeros(i,1) ];
    e      = [ eleft, e, eright ];
  end
  d = [];
  ind = -ncorrs:ncorrs;
  nc2 = length( Q.CALINOISE_CORR );
  for i = -nc2:nc2
    d = [ d, i*nf + ind];
  end
  e  = repmat(e,nza,1);
  e0 = e;
  for i = 1:nc2
    e      = [ e0*Q.CALINOISE_CORR(i), e, e0*Q.CALINOISE_CORR(i) ];
  end
  Se = Se + spdiags( e, d, nf*nza, nf*nza );
end

end %=== Old part