%------------------------------------------------------------------------
% NAME:     qp_x2arts
%
%           Maps a state vector to ARTS variables.
%
%           The function can be used for iterative retrievals to apply the
%           present state.
%
%           If the RECALC_ABS=0, the absorption is scaled to match the
%           retrieved species profiles. If RECALC_ABS=1, only the species
%           profiles are stored and the absorption must be recalculed in
%           next iteration.
%
%           The function is developed for qpcls but can maybe be useful
%           for other purposes.
%
% FORMAT:   qp_x2arts( Q, x, basename, Kx, kx_names, kx_index, kx_aux, 
%                          recalc_abs, logspecies, p_abs, vmrs, Abs, Absp, za )
%
% IN:       Q           Q structure.
%           x           The state vector.
%           basename    Basename for the ARTS calculations.
%           Kx          Weighting function matrix.
%           kx_names    As the corresponding ARTS variable.
%           kx_index    Start and stop index for each retrieval identity.
%                       A 2 column matrix.
%           kx_aux      As the corresponding ARTS variable.
%           recalc_abs  Flag to determine if absorption will be recalculated
%           logspecies  Flag for retrieval of the log. of species.
%           p_abs       As ARTS variable p_abs.
%           f_mono      As ARTS variable f_mono.
%           vmrs        As ARTS variable vmrs.
%           Abs         As ARTS variable abs.
%           Absp        As ARTS variable abs_per_tg.
%           za          As ARTS variable za_pencil.
%------------------------------------------------------------------------

% HISTORY: 2001.01.04  Started by Patrick Eriksson.
%          2003.11.04  Updated by Sho Tsujimaru.

function qp_x2arts( Q, x, basename, Kx, kx_names, kx_index, kx_aux, recalc_abs, logspecies, p_abs, f_mono, vmrs, Abs, Absp, za )


yc       = 0;
ispecies = 0;
nf       = length( f_mono );
Cabs     = [];                % Matrix for continuum absorption


i = 0;
%
while i < size(kx_names,1)

  i = i + 1;

  [group,name,par1] = qp_kinfo( kx_names{i} );

  ind = kx_index(i,1):kx_index(i,2);
   
 
  %=== Species
  %
  if strcmp( group, 'Species' )
    %
    ispecies = ispecies + 1;
    %
    xp = interpp( kx_aux(ind,1), x(ind), p_abs );
    %
    if logspecies( ispecies )
      xp = exp( xp );
    end
    %
    if Q.CLS_ABS_SPECIES_ON
      vmrs_old = vmrs;
      vmrs(ispecies,:) = xp';
      if ~recalc_abs
        % Calculate fraction of state vector x -- if any element of
        % vmrs is zero, recalc_abs has to be 1!!!
        xp = xp./vmrs_old';
      end
    else
      vmrs(ispecies,:) = vmrs(ispecies,:) .* xp'; 
    end
    
    %
    if ~recalc_abs
      %
      dAbs = Absp{ispecies} .* ( ones(nf,1) * (xp-1)' );
      Absp{ispecies} = Absp{ispecies} + dAbs; 
      Abs = Abs + dAbs;
      %
    end


  %=== Temperature
  %
  elseif strcmp( group, 'Temperature' )
    %
    t_abs = interpp( kx_aux(ind,1), x(ind), p_abs );
    write_artsvar( basename, 't_abs', t_abs, 'b' );


  %=== Continuum fit
  %
  elseif strcmp( group, 'Continuum' )
    %
    [ flimits, fp, f_mono ] = qp_contabs_freqs( Q );
    %
    indf = find( f_mono>=flimits(1) & f_mono<=flimits(2) );
    indp = find( p_abs<=kx_aux(ind(1),1) & p_abs>=kx_aux(last(ind),1) );
    % 
    np = length(fp);
    nx = length(ind);
    Xv = zeros( np, nx );
    %
    for ii = 1 : length(fp)
      ind      = kx_index(i+ii-1,1):kx_index(i+ii-1,2);
      Xv(ii,:) = 1e-3 * x(ind)';                      % conversion 1/km to 1/m
    end
    %
    fp = fp' /1e9;
    fm = f_mono/1e9;
    %
    Xe = zeros( length(indf), nx );
    %
    for ii = 1:nx
      Xe(:,ii) = polyval( polyfit( fp, Xv(:,ii), Q.CONTABS_ORDER ), fm(indf) );
    end
    %
    Cabs            = zeros( nf, length(p_abs) );
    Cabs(indf,indp) = interpp( kx_aux(ind,1), Xe', p_abs(indp) )';   
    %
    i = i + np - 1;
    %
    clear fp fm Xv Xe indf indp


  %=== Pointing off-set
  %
  elseif strcmp( name, 'Pointing: off-set' )
    %
    za = za + x(ind);
    write_artsvar( basename, 'za_pencil', za, 'b' );


  %=== Frequency off-set
  %
  %  MHz is used here as internal unit.
  % 
  elseif strcmp( name, 'Frequency: off-set' )
    %
    f_mono = f_mono + x(ind) * 1e6;
    write_artsvar( basename, 'f_mono', f_mono, 'b' );


  %=== Ground emission
  %
  elseif strcmp( group, 'Ground emission' )
    %
    x(ind) = check_minmax( x(ind), Q.EGROUND_MINMAX );
    %
    e_ground = read_artsvar( basename, 'e_ground' );
    f_mono   = read_artsvar( basename, 'f_mono' );
    [ fp, findex ] = qp_eground_lims( f_mono, Q.EGROUND_LIMITS );
    %
    for ii = 1 : length(fp)
      is = findex(ii,1) : findex(ii,2); 
      e_ground(is) = x(ind(ii));
    end
    write_artsvar( basename, 'e_ground', e_ground, 'b' );


  %=== Baseline fits
  %
  elseif strcmp( group, 'Baseline' )
    %
    yc = yc + Kx(:,ind) * x(ind);


  %=== Quantities handled by qp_x2y are included here so it can be checked
  %=== that all quantities are handled.
  %
  elseif strcmp( group, 'Baseline' )
    %

  elseif strcmp( name, 'Reference load temperature(s)' )
    %

  elseif strcmp( name, 'Proportional calibration error' )
    %

  %
  % Newly appended part. 
  %
  
  %= Pressure Shift
  %
  elseif strcmp( name, 'Pressure shift' )
     %
     % Scaling factor due to unit change: MHz/hPa --> Hz/Pa
     scfac = 1e4; 
     psf = x(ind) * scfac;
     
     fname  =  fullfile( Q.SPECTRO_DIR, Q.PSF_FILE );
     Lpsf   =  read_linefile(fname);
     n_Lpsf =  length(Lpsf);

     if ~n_Lpsf
       error('No line is found in PSF_FILE.'); 
     end

     fname = fullfile( Q.SPECTRO_DIR, Q.LINEFILE );
     L     = read_linefile( fname );

     %  
     fname0 = fullfile( Q.SPECTRO_DIR, 'backup_linefile' );
     if ~exist( fname0, 'file' ) 
        fname0 = fullfile( Q.SPECTRO_DIR, Q.LINEFILE );
     end
     %    

     L0     = read_linefile( fname0 );
	
     nl=[];
     for jj = 1:length(Lpsf) 
       n = find_line_num( L0, Lpsf{jj}.f );
       nl = [nl, n];    
     end
    
     if n_Lpsf ~= length(nl)
       error('The number of lines is not consistent.'); 
     end

     for jj = 1:n_Lpsf
       mm =  nl(jj);
       L{mm}.psf = psf;
     end

     write_linefile(fname, L); 

  %= Frequency Calibration
  %
  %  MHz is used as internal unit. 
  %
  elseif strcmp( name, 'Freq. calibration' )
     %
     if Q.FREQ_CALIB_ORDER == 0 
       %
       f_mono = f_mono + x(ind) * 1e6;
       %  
     elseif Q.FREQ_CALIB_ORDER > 0
       %
       fpoints = kx_aux(ind, 1);     
       npoints = length(fpoints);    
       df    = x(ind) * 1e6;         

       [m,n]=size(f_mono);
       f_tmp = zeros(m,n); 
       %
       for ipoint=1:npoints
         %
	 f_dist = ones(m,n); 
         %
         for ip=1:npoints
           if ip ~= ipoint         
	     f_dist ... 
	       = f_dist.*(f_mono-fpoints(ip))/(fpoints(ipoint)-fpoints(ip)); 
           end 
         end
         %
         f_tmp = f_tmp + f_dist * df(ipoint);
       end
       f_mono = f_mono + f_tmp; 
       %            
     else 
       error('Polynomial order should not be negative.');
     end
     write_artsvar( basename, 'f_mono', f_mono, 'b' );
 
  else
     error(sprintf( 'Unknown retrieval identity: %s', name ));
  end
end % while

%=== Save profile/absorption data
%
write_artsvar( basename, 'vmrs', vmrs, 'b' );
%
if ~recalc_abs
  if ~isempty( Cabs )
    Abs = Abs + Cabs;
  end
  write_artsvar( basename, 'abs', Abs, 'b' );
  write_artsvar( basename, 'abs_per_tg', Absp, 'b' );

else
  if isempty( Cabs )
    Cabs = zeros( nf, length(p_abs) );
  end
  write_artsvar( basename, 'abs0', Cabs, 'b' );

end
  

function x = check_minmax( x, minmax );

  if isempty( minmax )
    return
  end

  if length(minmax) ~= 2
    error( ...
        'The length of a vector with min and max values must have length 2.');
  end

  i = find ( x < minmax(1) );
  if ~isempty(i)
     x(i) = minmax(1);
  end

  i = find ( x > minmax(2) );
  if ~isempty(i)
     x(i) = minmax(2);
  end

return

%------------------------------------------------------------------------------

function n  = find_line_num( L, freq )

nn = 0;

for i = 1 : length(L)
  if L{i}.f == freq 
       n = i;
      nn =  nn + 1;
  end
end

if nn < 1
  error('No line is found!');
elseif nn > 1
  error(sprintf('A line is identified with more than %d lines.', nn));
end

return