%------------------------------------------------------------------------
% 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.08.31  Created by by Patrick Eriksson.


function [Kx,kx_names,kx_index,kx_aux] = qp_Kx( Q, basename, H, Hd, fy )


%=== Load Kx and check if size matches fy
%
Kx = read_artsvar( basename, 'kx' );
%
if size(Kx,1) ~= length(fy)
  error('Sizes of Kx and F(y) do not match (recalculation needed?).')
end 


%=== Apply H
%
Kx = H*Kx;



%=============================================================================
%=============================================================================
%=== Post-sensor parts og Kx
%=============================================================================
%=============================================================================

if Q.POLYFIT_DO == 3

  f  = read_datafile( [Q.SENSOR_DIR,'/',Q.F_SENSOR], 'VECTOR' );
  za = read_datafile( [Q.SENSOR_DIR,'/',Q.ZA_SENSOR], 'VECTOR' );

  nf  = length( f );
  nza = length( za );
  np  = Q.POLYFIT_ORDER + 1;

  %= Scale frequencies
  fsc = -1 + 2*(f-f(1)) / (last(f)-f(1));

  Kpol = zeros( nf*nza, nza*np );

  col    = 0;
  findex = 1:nf;

  kpost_names = cell( np, 1 );

  for i = 0:Q.POLYFIT_ORDER

    if i == 0
      k = ones(nf,1);
    else
      k = fsc.^i;
    end
  
    for j = 1:nza
      col = col + 1;
      rindex = findex + (j-1)*nf;
      Kpol(rindex,col) = k;
    end

    kpost_names{i+1} = sprintf('Baseline polynomial: %d',i); 
  end

  Kpost = Kpol;
  clear Kpol

  kpost_aux = [ repmat(vec2col(za),np,1), zeros(np*nza,1) ];
  kpost_lengths = repmat( nza, np, 1 );

else
  Kpost = [];
  kpost_names = [];
  kpost_aux = [];
  kpost_lengths = [];
end



if Q.PPOLYFIT_DO == 3

  f  = read_datafile( [Q.SENSOR_DIR,'/',Q.F_SENSOR], 'VECTOR' );
  za = read_datafile( [Q.SENSOR_DIR,'/',Q.ZA_SENSOR], 'VECTOR' );

  nf  = length( f );
  nza = length( za );
  np  = Q.PPOLYFIT_ORDER + 1;
  npp = length(Q.PPOLYFIT_LIMITS) - 1;

  Kpol = zeros( nf*nza, nza*np*npp );

  col    = 0;

  ktmp_names = cell( np*npp, 1 );

  %
  for l = 1:npp 

    pind = find( (f>=Q.PPOLYFIT_LIMITS(l)) & (f<Q.PPOLYFIT_LIMITS(l+1)) );

    %= Scale frequencies
    fsc = -1 + 2*(f(pind)-f(pind(1))) / (last(f(pind))-f(pind(1)));

    for i = 0:Q.PPOLYFIT_ORDER

      if i == 0
        k = ones(length(pind),1);
      else
        k = fsc.^i;
      end
  
      for j = 1:nza
        col = col + 1;
        rindex = pind + (j-1)*nf;
        Kpol(rindex,col) = k;
      end

      ktmp_names{ (l-1)*np + i + 1 } = ...
                     sprintf('Baseline piecewise polynomial: %d/%d',i,l); 

    end
  end
keyboard
  Kpost = [ Kpost, Kpol ];
  clear Kpol

  kpost_names = [ kpost_names; ktmp_names ];
  kpost_aux = [ kpost_aux; repmat(vec2col(za),np*npp,1), zeros(np*npp*nza,1) ];
  kpost_lengths = [ kpost_lengths; repmat( nza, np*npp, 1 ) ];

end



if ~isempty( Kpost )

  Kx = [ Kx, Hd*Kpost ];

end



%=============================================================================
%=============================================================================

%=== Set up help variables
%
kx_aux     = read_artsvar( basename, 'kx_aux' );
kx_names   = read_artsvar( basename, 'kx_names' );
kx_lengths = read_artsvar( basename, 'kx_lengths' );
%
if ~isempty( Kpost )

  kx_aux     = [ kx_aux; kpost_aux ];
  kx_lengths = [ kx_lengths; kpost_lengths ];
  kx_names   = [ kx_names; kpost_names ];
end
%
nrq    = size( kx_names, 1 );
%
kx_index = [ [1;1+cumsum(kx_lengths(1:(nrq-1)))], cumsum(kx_lengths) ];