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


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



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



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

if do_reduction
  %
  [H,f_y,za_y,Hd] = qp_Hd( Q, H, Hd, f_y, za_y );
  %
end



%=== Save if not return call
if nargout == 0
  save( [Q.OUT,'.h'], 'H' );
  save( [Q.OUT,'.hd'], 'Hd' );
  % Are these variables neeeded and how shall they be saved ?
  %write_datafile( [Q.MODEDIR,'/f_y.ab'], f_y, 'VECTOR' );
  %write_datafile( [Q.MODEDIR,'/za_y.ab'], za_y, 'VECTOR' );
  %write_datafile( [Q.MODEDIR,'/f_sensor.ab'], f_sensor, 'VECTOR' );
  %write_datafile( [Q.MODEDIR,'/za_sensor.ab'], za_sensor, 'VECTOR' );
end