%------------------------------------------------------------------------
% 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.
%

if ~find(strcmpi(Q.CLS_WFS_UNIT, {'vmr', 'frac'}))
  error('Allowed values for Q.VMR_WFS_UNITS are vmr and frac');
end
  
%=== The unit of the k matrix
if Q.CLS_ABS_SPECIES_ON | strcmpi(Q.CLS_WFS_UNIT, 'vmr') 
  k_unit = 'vmr';
else
  k_unit = 'frac';
end

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


%= Species
%
if in_klevels( 3, klevels )  &  ~isempty( Q.RETRIEVAL_TAGS )
  %
  nspecies = sstring_length( Q.RETRIEVAL_TAGS );
  %
  for i = 1:nspecies
    name    = sstring_get_i( Q.SPECIES_KGRIDS, i );
    grid    = read_datafile( fullfile( 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( fullfile( 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


%= Frequency
%
if in_klevels( Q.FREQUENCY_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( fullfile( 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( fullfile( 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 )  &  ~isempty( Q.RETRIEVAL_TAGS )
    %
    for i = 1:nspecies
      name = fullfile( 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  = "%s"\n}\n',...
						       		         name, k_unit);
      fprintf(fid, 'kxPutInK{\n}\n');
    end
  end
  
  
  %= Temperature
  if in_klevels( Q.TEMPERATURE_DO, klevels )
    name = fullfile( Q.RETRIEVDEF_DIR, Q.TEMPERATURE_KGRID );
    fprintf(fid, 'VectorReadAscii( k_grid ) {\n  "%s"\n}\n', name );
    fprintf(fid, 'kTemp{\n   hse  = %d\n   fast = %d\n}\n', ...
                                      Q.HSE_RETRIEVAL_ON, Q.TEMPERATURE_FAST );
    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


  %= Frequency
  %
  %  MHz is used here as internal unit.
  %
  if in_klevels( Q.FREQUENCY_DO, klevels )
    fprintf(fid, 'kFrequencyOffSet{\n   delta = %.2e\n',...
                                                        Q.FREQUENCY_DELTA/1e6);
    fprintf(fid, '   f_unit = "MHz"\n}\n' );
    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)
  %
  QE = [];
  QE = qp_hdf( Q, QE );
  %
  fprintf(fid, 'MatrixWrite%s( kx ) { \n""\n}\n', QE.FORMAT );  
  fprintf(fid, 'ArrayOfIndexWrite%s( kx_lengths ){ \n""\n}\n', QE.FORMAT );  
  fprintf(fid, 'MatrixWrite%s( kx_aux ){ \n""\n}\n', QE.FORMAT );  
  fprintf(fid, 'ArrayOfStringWrite%s( kx_names ){ \n""\n}\n', QE.FORMAT );

end


 %= Spectro Parameters
 %= Check if the calculation of the spectroscopic parameters is requested
 if ( Q.SPECTRO_DO ==1 ) 
    fprintf(fid, 'kbInit{\n}\n');
   fprintf(fid, 'kSpectro{\n do_intens=%d\n',  Q.INTENS_ON);
   fprintf(fid, 'do_position=%d\n',  Q.POSITION_ON);
   fprintf(fid, 'do_agam=%d\n',  Q.AGAM_ON);
   fprintf(fid, 'do_sgam=%d\n',  Q.SGAM_ON);
   fprintf(fid, 'do_nair=%d\n',  Q.NAIR_ON);
   fprintf(fid, 'do_nself=%d\n',  Q.NSELF_ON);
   fprintf(fid, 'do_pSift=%d\n}\n',  Q.PSHIFT_ON);
   
  %= Convert to RJ temperatures if emission observations
  %
  if Q.EMISSION_ON
    fprintf(fid, 'MatrixTRJ( k, k ) { \n}\n');
  end 
   
   fprintf(fid, 'MatrixWrite%s( k) { \n""\n}\n', QE.FORMAT );
   fprintf(fid, 'ArrayOfStringWrite%s( k_names ){ \n""\n}\n', QE.FORMAT );
   fprintf(fid, 'MatrixWrite%s(S_S) { \n""\n}\n', QE.FORMAT );
 end 
  % 
return



function b = in_klevels(level,klevels)

  b = any( level == klevels );

return