%------------------------------------------------------------------------
% NAME:     qp_assemble_K
%
%           Loads a weighting function matrix calculated by ARTS and adds
%           weighting functions for quantities not handled by ARTS.
%           
%
% FORMAT:   [ K, k_names, k_index, k_aux ] = ... 
%                                  qp_assemble_K( Q, basename, H, Hd, klevels )
%
% OUT:      K          Weighting function matrix.
%           k_names    Name of each retrieval identity.
%           k_index    Index of first and last element in x for each
%                      retrieval identity.
%           k_aux      As kx_aux in ARTS.
% IN:       Q          Setting structure.
%           basename   Basename for the performed ARTS run.
%           H          The sensor/data-reduction transfer matrix.
%           Hd         The data-reduction transfer matrix.
%           klevels    Vector with the do-levels to consider.
%                      See further qp_iter1.m.
%------------------------------------------------------------------------

% HISTORY: 2001.08.31  Created by by Patrick Eriksson.


function [K,k_names,k_index,k_aux] = qp_assemble_K(Q,basename,H,Hd,klevels)


%
% Note that if the order is changed here, qp_Sxb must be changed 
% accordingly.
%



%=== Load K varaibles(K is always stored as Kx)
%
K         = read_artsvar( basename, 'kx' );
k_aux     = read_artsvar( basename, 'kx_aux' );
k_names   = read_artsvar( basename, 'kx_names' );
k_lengths = read_artsvar( basename, 'kx_lengths' );



%=== Post-processing for some variables
%
%= Ground emission
if in_klevels( Q.EGROUND_DO, klevels )
  %
  f_mono = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.F_MONO], 'VECTOR' );
  [ fp, findex ] = qp_eground_lims( f_mono, Q.EGROUND_LIMITS );
  %
  if length(fp) ~= length(f_mono)
    %
    nf = length( f_mono );
    %
    i0 = 1;
    while ~strcmp(  qp_kinfo(k_names{i0}), 'Ground emission' )
      i0 = i0 + 1;
    end
    %
    error('Implementation not finished');
    %
  end
end
%-----------------------------------------------------------------------------



%=== Apply H
%
K = H * K;





%=============================================================================
%=============================================================================
%=== Parts of K not covered by ARTS
%=============================================================================
%=============================================================================


if in_klevels( Q.POLYFIT_DO, klevels )

  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 in_klevels( Q.PPOLYFIT_DO, klevels )

  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

  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






%=============================================================================
%=== Put parts together
%=============================================================================


%=== K
%
if ~isempty( Kpost )
  %
  K = [ K, Hd*Kpost ];
  %
end


%=== Help variables
%
if ~isempty( Kpost )
  k_aux     = [ k_aux; kpost_aux ];
  k_lengths = [ k_lengths; kpost_lengths ];
  k_names   = [ k_names; kpost_names ];
end
%
nrq    = size( k_names, 1 );
%
k_index = [ [1;1+cumsum(k_lengths(1:(nrq-1)))], cumsum(k_lengths) ];


return
%------------------------------------------------------------------------------


function b = in_klevels(level,klevels)

  b = any( level == klevels );

return