%------------------------------------------------------------------------
% NAME:     qpi_K
%
%           Inline function to include calculation of a K matrix in
%           the control file.
%
% FORMAT:   qpi_K( Q, fid, klevels )
%
% OUT:      -
% IN:       Q        Setting structure
%           fid      File identifier for the control file being created.
%           klevels  Vector with the do-levels to consider.
%                    See further qp_iter1.m.
%------------------------------------------------------------------------

% HISTORY: 2001.03.28  Created by Patrick Eriksson
% HISTORY: 2001.10.13  KLEVELS introduced by Patrick Eriksson


function qpi_K( Q, fid, klevels )


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


%=== Determine the length of x ================================================
%
nx = 0;
ni = 0;


%= Species
%
if in_klevels( 3, klevels )
  %
  nspecies = sstring_length( Q.RETRIEVAL_TAGS );
  %
  for i = 1:nspecies
    name    = sstring_get_i( Q.SPECIES_KGRIDS, i );
    grid    = read_datafile( [ Q.RETRIEVDEF_DIR, '/', name ], 'VECTOR');
    nx      = nx + length(grid);
    ni      = ni + 1;
  end
end


%= Temperature
%
if in_klevels( Q.TEMPERATURE_DO, klevels )
  name    = Q.TEMPERATURE_KGRID;
  grid    = read_datafile( [ Q.RETRIEVDEF_DIR, '/', name ], 'VECTOR');
  nx      = nx + length(grid);
  ni      = ni + 1;
end


%= Pointing
%
if in_klevels( Q.POINTING_DO, klevels )
  nx = nx + 1;
  ni = ni + 1;
end


%= Continuum absorption
%
if in_klevels( Q.CONTABS_DO, klevels )
  name    = Q.CONTABS_KGRID;
  grid    = read_datafile( [ Q.RETRIEVDEF_DIR, '/', name ], 'VECTOR');
  nx      = nx + length(grid)*(Q.CONTABS_ORDER+1);
  ni      = ni + Q.CONTABS_ORDER + 1;
end


%= Ground emission
%
if in_klevels( Q.EGROUND_DO, klevels )
  f_mono  = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.F_MONO], 'VECTOR' );
  [ fp, findex, single_e ] = qp_eground_lims( f_mono, Q.EGROUND_LIMITS );
  if single_e
    nx      = nx + 1;
  else
    nx      = nx + length(f_mono);
  end
  ni      = ni + 1;
end




%=== Write to cfile ===========================================================

if nx

  %= Allocate memory for Kx
  %
  fprintf(fid, 'kxAllocate( y ){\n   ni = %d\n   nx = %d\n}\n', ni, nx ); 
  
  
  %= Calculate the abs. LOS WFs (if needed)
  if in_klevels( 3, klevels ) | in_klevels( Q.CONTABS_DO, klevels ) | ...
	(in_klevels( Q.TEMPERATURE_DO, klevels ) & ~Q.HSE_RETRIEVAL_ON)
    fprintf(fid, 'absloswfsCalc{\n}\n');
  end
  
  %= Species
  if in_klevels( 3, klevels )
    %
    for i = 1:nspecies
      name = [ Q.RETRIEVDEF_DIR, '/', sstring_get_i( Q.SPECIES_KGRIDS, i ) ];
      fprintf(fid, 'VectorReadAscii( k_grid ) {\n  "%s"\n}\n', name );
      name = sstring_get_i( Q.RETRIEVAL_TAGS, i );
      fprintf(fid, 'kSpeciesSingle{\n   tg   = "%s"\n   unit  = "frac"\n}\n', ...
									   name);
      fprintf(fid, 'kxPutInK{\n}\n');
    end
  end
  
  
  %= Temperature
  if in_klevels( Q.TEMPERATURE_DO, klevels )
    name = [ Q.RETRIEVDEF_DIR, '/', Q.TEMPERATURE_KGRID ];
    fprintf(fid, 'VectorReadAscii( k_grid ) {\n  "%s"\n}\n', name );
    fprintf(fid, 'kTemp{\n   hse  = %d\n   fast = 1\n}\n', Q.HSE_RETRIEVAL_ON);
    fprintf(fid, 'kxPutInK{\n}\n');
  end
  
  
  %= Pointing
  %
  if in_klevels( Q.POINTING_DO, klevels )
    fprintf(fid, 'kPointingOffSet{\n   delta = %.2e\n}\n',Q.POINTING_DELTA);
    fprintf(fid, 'kxPutInK{\n}\n');
  end
  
  
  %= Continuum absorption 
  %
  if in_klevels( Q.CONTABS_DO, klevels )
    lims = qp_contabs_freqs(Q);
    name = fullfile( Q.RETRIEVDEF_DIR, Q.CONTABS_KGRID );
    fprintf(fid, 'VectorReadAscii( k_grid ) {\n  "%s"\n}\n', name );
    fprintf(fid, 'kContAbs{\n');
    fprintf(fid, '   order  = %d\n', Q.CONTABS_ORDER );
    fprintf(fid, '   f_low  = %.10e\n', lims(1) );
    fprintf(fid, '   f_high = %.10e\n', lims(2) );
    fprintf(fid, '   l_unit = "km"\n}\n' );
    fprintf(fid, 'kxPutInK{\n}\n');
  end
  
  
  %= Ground emission
  %
  if in_klevels( Q.EGROUND_DO, klevels )
    [ fp, findex, single_e ] = qp_eground_lims( f_mono, Q.EGROUND_LIMITS );
    fprintf(fid, 'kEground{\n');
    fprintf(fid, '   single_e = %d\n}\n', single_e );
    fprintf(fid, 'kxPutInK{\n}\n');
  end
  
  
  %= Convert to RJ temperatures if emission observations
  %
  if Q.EMISSION_ON
    fprintf(fid, 'MatrixTRJ( kx, kx ) { \n}\n');
  end
  
  
  %= Save as kx (with additional variables)
  %
  fprintf(fid, 'MatrixWriteBinary( kx ) { \n""\n}\n');  
  fprintf(fid, 'ArrayOfIndexWriteBinary( kx_lengths ){ \n""\n}\n');  
  fprintf(fid, 'MatrixWriteBinary( kx_aux ){ \n""\n}\n');  
  fprintf(fid, 'ArrayOfStringWriteBinary( kx_names ){ \n""\n}\n');

end

return



function b = in_klevels(level,klevels)

  b = any( level == klevels );

return