%------------------------------------------------------------------------
% NAME:     qp_iter1
%
%           Calculates spectrum and a weighting function matrix (K).
%
%           If only the spectra are of interest, it is better to use
%           QP_Y.
%
%           Which weighting function matrix that is calculated is controled
%           by KLEVELS (see options below).
%
%           The function creates a control file and a number of data files
%           that should be deleted by the calling function.
%
%           The spectrum vector (fy) and the weighitng function variables
%           are set to be empty if there is an ARTS crash.
%
% FORMAT:   [fy,K,k_names,k_index,k_aux,Kb,kb_names,cfile,basename] = ...
%                       qp_iter1(Q,H,Hd,tmpdir,klevels [,do_nonlin,recalc_abs])
%           
%
% OUT:      fy         Forward model spectrum (including H).
%           K          The state vector weighting function matrix.
%           k_names    The name of each retrieval identity.
%           k_index    Index of first and last element in x for each
%                      retrieval identity.
%           k_aux      As k_aux in ARTS.
%           k_aux_abs  The same as k_aus but this time is in
%           absolute units.
%           Kb         Spectroscopic parameters weighting function matrix.
%           kb_names   The name of each spectroscopic parameter.
%           cfile      Name on the control file used.
%           basename   Basename for the performed ARTS run.
% IN:       Q          Setting structure.
%           H          The sensor/data-reduction transfer matrix.
%           Hd         The data-reduction transfer matrix.
%           tmpdir     Name/path of temporary directory (shall exist).
%           klevels    Vector with the do-levels to consider.
%                      Kx -> klevels = 3
%                      Kb -> klevels = 2
%                      error characterisation -> klevels = [1 2]
% OPTIONAL: do_nonlin  Flag for linear/non-linear calculations.
%           recalc_abs Flag for doing recalculation of absorption.
% TEMPLATE: iter1.tmplt
%------------------------------------------------------------------------

% HISTORY: 2001.03.28  Started by Patrick Eriksson.
%          2001.10.13  KLEVELS introduced by Patrick Eriksson


function [fy,K,k_names,k_index,k_aux,k_aux_abs,Kb,kb_names, cfile, basename] = ...
    qp_iter1(Q,H,Hd,tmpdir,klevels,do_nonlin,recalc_abs);


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


%### Check that KLEVELS has allowed value(s) ##################################
%
if ~( all(klevels==3) | all(klevels==2) | all(klevels==[1 2]) ) 
    error('Allowed selections for KLEVELS are 3, 2 and [1 2]');
end


%### Create cfile #############################################################
%
QE.DO_NONLIN  = do_nonlin;
QE.RECALC_ABS = recalc_abs;
QE            = qp_hdf( Q, QE );
%
QE.KLEVELS    = 0;
for i = 0 : (length(klevels)-1)
    QE.KLEVELS = QE.KLEVELS + 10^i * klevels(length(klevels)-i);
end
%
template      = which( 'iter1.tmplt' );
%
[cfile,basename] = qtool( Q, tmpdir, template, QE );


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

%### Load y ###################################################################
%
yfilename = [ basename, '.y', qp_hdf(Q) ];



if exist( yfilename, 'file' )
    fy = read_artsvar( basename, 'y' );
else
    fy      = [];
    K       = [];
    k_names = [];
    k_index = [];
    k_aux   = [];
    return
end



%### Check size ###############################################################
%
if  max(size(H)) > 1 & length(fy) ~= size(H,2)
    error('Length of y and size of H do not match (recalculate H?).')
end


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


%### Get K ###################################################################
%
[K,k_names,k_index,k_aux,Kb,kb_names] = qp_assemble_K(Q, basename, ...
						  H, Hd, fy,0,klevels);

% get the kx_aux in absolute value. This is necessary for plot.
VMRs = read_artsvar( basename, 'vmrs' ); % vmrs profiles
k_aux_abs= k_aux; % initializing
nrq      = size( k_names, 1 );
ispecies =0;
p_abs= read_datafile(fullfile(Q.CALCGRIDS_DIR, Q.P_ABS),  'Matrix');

for i = 1:nrq
    
    [group,name] = qp_kinfo(  k_names{i} ); 
    
    ind = k_index(i,1):k_index(i,2);
    
    %=== Species
    % only the species are in relative units.
    if strcmp( group, 'Species' ) 
        %
        ispecies     = ispecies + 1;
        k_aux_abs(ind,1)   = k_aux(ind,1);
        k_aux_abs(ind,2)  = interpp( p_abs, VMRs(ispecies,:)', k_aux_abs(ind,1));
    else
        k_aux_abs(ind,1) = k_aux(ind,1);
        k_aux_abs(ind,2) = k_aux(ind,2);
    end
end