%------------------------------------------------------------------------
% NAME:     qpoem_iter1
%
%           Calculates spectrum and weighting function matrix (Kx) for
%           the first iteration of an OEM inversion. The H matrix is 
%           applied and is also an output of the function.
%
%           If DO_NONLIN=1, a number of quantities are stored (see the
%           used template file) for usage during later iterations.
%
% FORMAT:   [fy,H,Kx,kx_names,kx_index,kx_aux,cfile,basename] = ...
%                               qpoem_iter1(Q,tmpdir [,do_nonlin,recalc_abs] )
%
% OUT:      fy         Forward model spectrum (including H).
%           H          The sensor/data-reduction transfer matrix.
%           Kx         The state vector weighting function matrix.
%           kx_names   The name of each retrieval identity.
%           kx_index   Index of first and last element in x for each
%                      retrieval identity.
%           kx_aux     As kx_aux in ARTS.
%           cfile      Name on the control file used.
%           basename   Basename for the performed ARTS run.
% IN:       Q          Setting structure.
%           tmpdir     Name/path of temporary directory (shall exist).
% OPTIONAL: do_nonlin  Flag for linear/non-linear calculations.
%           recalc_abs Flag for doing recalculation of absorption.
% TEMPLATE: oem_iter1.tmplt
%------------------------------------------------------------------------

% HISTORY: 2001.03.28  Started by by Patrick Eriksson.


function [fy,H,Hd,Kx,kx_names,kx_index,kx_aux,cfile,basename] = ...
                                    qpoem_iter1(Q,tmpdir,do_nonlin,recalc_abs)


if ~exist('do_nonlin','var'),  do_nonlin  = 0; end
if ~exist('recalc_abs','var'), recalc_abs = 0; end



%### Load H and Hd ############################################################
%
load( [Q.OUT,'.h'], '-mat' );
load( [Q.OUT,'.hd'], '-mat' );



%### Create cfile #############################################################
%
QE.DO_NONLIN  = do_nonlin;
QE.RECALC_ABS = recalc_abs;
template      = which( 'oem_iter1.tmplt' );
%
[cfile,basename] = qtool( Q, tmpdir, template, QE );



%### Run ARTS and load results ################################################
%
qp_arts( Q, cfile ); 


%### Load y ###################################################################
%
fy = read_artsvar( basename, 'y' );


%### Get Kx ###################################################################
%
[Kx,kx_names,kx_index,kx_aux] = qp_Kx( Q, basename, H, Hd, fy );


%### Apply H on y #############################################################
%
fy = H*fy;