%------------------------------------------------------------------------
% NAME:     qp_Hd
%
%           Sets up Hd and includes this into H. 
%
%           See the README file for options.
%
% FORMAT:   [H,f_y,za_y,Hd] = qp_Hd(Q,H,Hd,f_y,za_y)
%
% OUT:      H
%           f_y
%           za_y
%           Hd
% IN:       Q              The setting structure.
%           H
%           Hd
%           f_y
%           za_y
%------------------------------------------------------------------------

% HISTORY: 2001.03.29  Created by Patrick Eriksson.


function [H,f_y,za_y,Hd] = qp_Hd(Q,H,Hd,f_y,za_y)



%=== Average of subsequent spectra
if Q.BINVIEW_ON
  %
  [H,f_y,za_y,Hd] = hBinView( H, f_y, za_y, Hd, Q.BINVIEW_DATA );
  %
end



%=== Averaging of neighbouring channels (binning)
if Q.BINNING_ON
  %
  [H,f_y,za_y,Hd] = hBinFile( Hd, f_y, za_y, Q.BINNING_FILE );
  %
end



%=== Data reduction by K
if Q.KRED_ON

  error('K-based data reduction not yet implemented')

  % The code below is old

  QE.ATMOSPHERE = Q.SETUP_ATM1;
  QE.DO_NONLIN  = 0;

  tmpdir    = temporary_directory( Q.TMP_AREA );

  [cfile,basename] = qsystem( Q, tmpdir, [Q.TEMPLATE_DIR,'/iter1.tmplt'], ...
    								    QE );
  %=== Run ARTS and load Kx
  qsmr_arts( Q, cfile ); 
  Kx = read_artsvar( basename, 'kx' );
  Kx = H*Kx;

  %=== Scale KX by elements of SX
  Dx = qsmr_Sx( Q );
  np = size( Dx, 1 );
  Kx = Kx * sqrt(Dx);

  %=== Calculate the reduction transfer matrix
  [U,S,V]    = svd( Kx, 0 );
  clear S V Kx
  if Q.DATARED_N > size(U,2)
    error(sprintf( ...
   '%d number of eigenvectors cannot be provided.\nOnly %d are possible.',...
    						 Q.DATARED_N,size(U,2)));
  end
  U    = U( :, 1:Q.DATARED_N )';
  H    = h_x_h( U, H );
  Hd   = h_x_h( U, Hd );
  clear U
  f_y  = 1:Q.DATARED_N;
  za_y = ones( Q.DATARED_N, 1 );  

  %=== Delete the temporary directory
  delete_tmp_dir( Q.TMP_AREA, tmpdir )
  
end