%------------------------------------------------------------------------
% NAME:     qp_Se
% 
%           Creates the Se matrix.
% 
%           If NARGOUT=0, the Se matrix and its inverse are saved to files. 
%           Else, the Se matrix is returned (not inverted).
% 
% FORMAT:   qp_Se( Q )
%             or
%           Se = qp_Se( Q )
% 
% IN:       Q          Setting structure.
% OUT:      Se         The Se matrix.
%
% TEMPLATE: -
%------------------------------------------------------------------------

% HISTORY: 2001.03.28  Created by Patrick Eriksson.
%          2001.08.25  Modified by Carlos Jimenez


function Se = qp_Se( Q )


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


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

%=== 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 measurement 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
%
load( [Q.OUT,'.hd'], '-mat' );
%
if ~(size(Hd,1)==1 & size(Hd,2)==1 & Hd==1)
  %
  out(2,'Applying data reduction on noise.');
  %
  Se = Hd * Se * Hd';
end
clear Hd



if nargout == 0

  %=== Invert
  Seinv = inv( Se );

  %=== Save
  save( [Q.OUT,'.se'], 'Se' );
  save( [Q.OUT,'.seinv'], 'Seinv' );

end


out(1,-1);
return


%=== This is an old part, not adapted to Qpack
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





%=== Remaining variables are handled by ARTS template
%
if Q.TEMPERATURE==2  | Q.POINTING==2   | Q.CALIBRATION==2 | ...
   Q.CABS_PRIMARY==2 | Q.CABS_IMAGE==2
  %
  tmpdir    = temporary_directory( Q.TMP_AREA );
  %    
  QE.ATMOSPHERE = Q.SETUP_ATM1;
  %
  [cfile,basename] = qsystem( Q, tmpdir, [Q.TEMPLATE_DIR,'/se.tmplt'], QE );
  %
  qsmr_arts( Q, cfile ); 
  %
  Sb = read_artsvar( basename, 'sb' );
  %
  Kb = read_artsvar( basename, 'kb' );
  load( [Q.MODEDIR,'/h'] );
  Kb = H * Kb;
  clear H
  %
  Se = Se + Kb * Sb * Kb';
  %
  delete_tmp_dir( Q.TMP_AREA, tmpdir )
  %
end

end %=== Old part