%------------------------------------------------------------------------
% NAME:     qp_H
%
%           Sets up H and Hd. The Hd matrix is treated in a seperate 
%           function (qp_Hd).
%
%           See the README file for options.
%
%           If the function is called without output arguments, H and Hd
%           are stored to files. With output arguments, no files are
%           created.
%
%           The calculations can be done for arbitrary P_ABS, F_MONO and 
%           ZA_PENCIL. If these variables not are given, they are read from 
%           Q.OUT.
%
% FORMAT:   q_H( Q )
%             or
%           [H,f_y,za_y,f_sensor,za_sensor] = 
%                        qp_H( Q [,p_abs,f_mono,za_pencil,do_reduction] )
%
% OUT:      (see above)
% IN:       -
% OPTIONAL: Q              The setting structure.
%           p_abs
%           f_mono
%           za_pencil
%           do_reduction   Flag to exclude or include all data reduction.
%
% TEMPLATE: -
%------------------------------------------------------------------------

% HISTORY: 2001.03.25  Created by Patrick Eriksson.
%          2001.09.02  Updated by Carlos Jimenez

function [H,f_y,za_y,f_sensor,za_sensor] = ...
                               qp_H(Q,p_abs,f_mono,za_pencil,do_reduction)


%=== Set default values for optional arguments
%
if  ~exist( 'do_reduction', 'var' )
  do_reduction = 1;
end



%=== Shall results be saved
%
do_save = ~nargout;



%=== Read f_mono and za_pencil if not input variables
if nargin == 1
  f_mono    = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.F_MONO], 'VECTOR' );
  za_pencil = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.ZA_PENCIL], 'VECTOR' );
end



%=== Init the H matrices
[H,f_y,za_y,Hd,f_sensor,za_sensor] = hInit( f_mono, za_pencil );



%=== Antenna pattern
if Q.ANTENNA_ON
  %
  za_obs = read_datafile( [Q.SENSOR_DIR,'/',Q.ZA_SENSOR], 'VECTOR' );
  %
  fname = [Q.SENSOR_DIR,'/',Q.ANTENNA_FILE];
  %
  [H,f_y,za_y,za_sensor] = hAntennaFromFileAdv( H, f_sensor, za_sensor, ...
         za_obs, fname, Q.ANTENNA_ORDER, Q.ZA_ORDER, 0, 0, 1, Q.ANTENNA_MOVE );
end



%=== Mixer and sideband filter
if Q.DSB_ON
  %
  fname = [Q.SENSOR_DIR,'/',Q.DSB_FILE];
  %
  [H,f_y,za_y,f_sensor] = hMixerFromFile( H, f_sensor, za_sensor, Q.DSB_LO,...
                               Q.DSB_FPRIMARY, fname, Q.DSB_ORDER, Q.F_ORDER );
end



%=== Backend
if Q.BACKEND_ON
  %
  f_obs = read_datafile( [Q.SENSOR_DIR,'/',Q.F_SENSOR], 'VECTOR' ); 
  %
  fname = [Q.SENSOR_DIR,'/',Q.BACKEND_FILE];
  %
  [H,f_y,za_y,f_sensor] = hBackendFromFile( H, f_sensor, za_sensor, ...
                                     f_obs, fname, Q.BACKEND_ORDER, Q.F_ORDER);
end



%=== saving sensor without reduction
%
if do_save
  Hs = H;
  save( [Q.OUT,'.hs'], 'Hs' );
  clear Hs
end



%=== Data reduction part ======================================================

if do_reduction & (Q.BINNING_ON | Q.KRED_ON |  Q.LRED_ON  | Q.BINVIEW_ON )
  %
  [H,f_y,za_y,Hd] = qp_Hd( Q, H, Hd, f_y, za_y );
  %
end



%=== Save if not return call and H is not empty
%    NOTE : H empty correspond to a H and Hd with reduction
%           specific for each retrieval point, where Hs and
%           Hds are saved as structures or individual Hs 
%           amd Hds in the right mscript, e.g. hRedLimb     
%
if do_save & ~isempty(H)
  save( [Q.OUT,'.h'], 'H' );
  save( [Q.OUT,'.hd'], 'Hd' );
end
%
%= Save frequency and zenith angle grids
%
if do_save
  save( [Q.OUT,'.f_sensor'], 'f_sensor');
  save( [Q.OUT,'.za_sensor'], 'za_sensor');
  save( [Q.OUT,'.f_y'], 'f_y');
  save( [Q.OUT,'.za_y'], 'za_y');
end