%------------------------------------------------------------------------ % 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,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 from Q.OUT ######################################################## % load( [Q.OUT,'.h'], '-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 ); % fy = read_artsvar( basename, 'y' ); Kx = read_artsvar( basename, 'kx' ); % if size(Kx,1) ~= length(fy) error('Sizes of Kx and F(y) do not match (recalculation needed?).') end % fy = H*fy; Kx = H*Kx; % kx_aux = read_artsvar( basename, 'kx_aux' ); kx_names = read_artsvar( basename, 'kx_names' ); kx_lengths = read_artsvar( basename, 'kx_lengths' ); nrq = size( kx_names, 1 ); kx_index = [ [1;1+cumsum(kx_lengths(1:(nrq-1)))], cumsum(kx_lengths) ];