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

global REPORT_LEVEL
r_l = Q.QP_LEVEL;


%=== Read some vectors to get sizes
%
if ~Q.ANTENNA_ON & ~Q.DSB_ON & ~Q.BACKEND_ON
  f  = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.F_MONO ], 'VECTOR' );
  za = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.ZA_PENCIL ], 'VECTOR' );
else
  f  = read_datafile( [Q.SENSOR_DIR,'/',Q.F_SENSOR ], 'VECTOR' );
  za = read_datafile( [Q.SENSOR_DIR,'/',Q.ZA_SENSOR ], 'VECTOR' );
end
%
nf        = length( f);
nza       = length( za );


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


%=== Measurement thermal noise

if Q.MEASNOISE_DO

  out(2,'Doing measurement thermal noise.');

  % -- getting covariance matrix
  % -- NOTE: Doing each za at one time to avoid a covariance
  %          matrix of dimension f_y x f_y, even play with sparse
  %          matrices is very slow

  REPORT_LEVEL = 0;
  se = sFromFile([Q.SENSOR_DIR,'/',Q.MEASNOISE_COVMAT],f);
  REPORT_LEVEL = r_l;

  % -- loading data reduction  
  load( [Q.OUT,'.hd'], '-mat' );

  n_red = size(Hd,1);

  if Hd==1
    % -- no reduction

    [dummy,p]= spdiags(se);
    if p==0  % is diagonal

      % -- save as diagonal matrix
      B= repmat(diag(se),nf*nza,1);
      Se = spdiags(B,0,nf*nza,nf*nza);

    else

      % -- fill the elements of a sparse_zero
      %    with possible run out of memory if
      %    sizes too large
      Se = sparse(zeros(nf*nza,nf*nza));
      for i= 1:nza
        ind          = (1+((nza-1)*nf)):(nf*nza);
        Se(ind,ind)  = se;
      end 
      % -- make full if more than 2/3 of values are used
      if length(find(H)) > size(H,1)*size(H,2)*0.67
        Se = full(Se);
      end
    end

  else

    % -- reduction
    %    doing Hd * Se * Hd' in pieces, it seems
    %    faster than using a big sparse matrix and
    %    filling the non-zero elements for big
    %    sizes, right?
    %    It is done by constructing Hd*Se by filling
    %    Ha taking advantage of the structure of Se
    %    because no correlation in za is allowed here. 

    n_red = size(Hd,1);
    Ha    = zeros(n_red,nf*nza);
    for i = 1:nza
      ind = (1+((i-1)*nf)):(nf*i);
      Ha(:,ind)  = Hd(:,ind) * se;
    end
    Se    = Ha * Hd';

  end

end


clear Hd

%=== Calibration thermal noise
%
if Q.CALINOISE_DO == 2
  %
  % -- old part, no updated as calibartion still not clear
  %
  out(2,'Doing calibration thermal noise not implemeneted');   
  %
end

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



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


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