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


function Se = qp_Se( Q )



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



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


Se = 0;


%=== Common part for measurement and calibration noise
if Q.CHCORR_ON
  ncorrs = length( Q.CHCORR_DATA );
else
  ncorrs = 0;
end




%=== Measurement thermal noise
%
if 1
  %
  out(2,'Doing measurement thermal noise.');
  %
  Noise  = read_datafile( [Q.SENSOR_DIR,'/',Q.MEASNOISE_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
  Se       = spdiags( repmat(e,nza,1), -ncorrs:ncorrs, nf*nza, nf*nza );
end



%=== Calibration thermal noise
%
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



%=== Include data reduction
%
if ~( (size(Se,1)==1 & size(Se,2)==1) )
  %
  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
end


%=== This is an old part, not adapted to Qpack
if 0

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


if nargout == 0

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

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

end


out(1,-1);