%------------------------------------------------------------------------------
% NAME:     qpopt_linefile
%
%           Selects transitions to include in a linefile by pure numerical
%           force.
%
%           This function performs rather detailed calculations, but can
%           be relatively slow for molecules with large number of transitions,
%           such as HNO3.
%
%           The calculations are done as:
%           1. All transitions of the molecule inside DF_SEARCH from the 
%              primary and image bands are loaded.
%           2. If DF_COMBINE has a value, transitions are combined.
%           3. For each transition the following is done:
%              a. If the transition is inside the primary band, the line
%                 frequency is put into F_MONO. Otherwise, thye closest
%                 end of the primary band is included.
%              b. If the transition is closer to the image band than the
%                 primary band, F_MONO is extended following the procedure in
%                 a but for the image band.
%              c. Spectra are calculted for ZA_PENCIL and F_MONO and it is
%                 checked if any value is >= TB_LIMIT. Values for the image
%                 band are multiplicated with SIDEBAND_RATIO before comparision
%                 to TB_LIMIT.
%              d. If no value is above TB_LIMIT, c is repeated for next set-up
%                 atmosphere. Otherwise, go to e. 
%              e. If TB_OUT is an output argument, c is performed for all set-
%                 up atmospheres and the maximum value is styored in TB_MAX.
%                 Otherwise, take next transition.
%           4. Points 1-3 are performed for all molecules.
%           
%           If the function is called without output argument, the linefile
%           is stored to the corresponding file specified by Q. 
%
%           Note that for weak transitions, the printed brightness temperature
%           (requires Qpack report level 3) can be negative due to a non-
%           perfect compensation of cosmic background radiation.
%
%           Handled line formats (LFORMAT shall equal the name given below):
%
%           VERDANDI 
%              Molecules are given in the Verdandi tag format, for example
%              31 for O3-666.
%
% FORMAT:   [ L, tb_max ] = qpopt_linefile( Q, f_range, lo, sideband_ratio,
%              za_pencil, tb_limit, df_search, molecules, ldir, lformat, 
%              [ df_combine, nsetup ] )
%
% OUT:      L                Structure array with line data (see 
%                            WRITE_LINEFILE).
%           tb_max           Max intensity found for the transitions in L for
%                            all considered set-up atmospheres.
%                            The intensities are given as contribution to
%                            the primary band. 
% IN:       Q                Qpack struct.
%           f_range          Frequency range of primary band as [flow,fhigh].
%           lo               LO frequency. Setting LO=[] means that no
%                            image band shall be considered.
%           sideband_ratio   Relative contribution from the sideband.
%           za_pencil        Zenith angle grid to be used.
%           tb_limit         Intensity threshold (in K) for transitions .
%                            The threshold is specified for the primary band.
%                            Contributions to the image band is multiplicated
%                            with SIDEBAND_RATIO.
%           df_search        How far away from primary and image band that
%                            transitions search shall be extended.
%           molecules        The molecules (isotopes) for which search shall
%                            be performed. The way the molecules are specified
%                            follows LFORMAT.
%           ldir             Directory with spectroscopic files. One file per
%                            item in MOLECULES.
%           lformat          String describing used line format. See above.
% OPTIONAL: df_combine       If this argument is defined and not empty, 
%                            "identical" transitions are com,bined using 
%                            COMBINE_LINES where DF_COMBINE is passed on as
%                            DF.
%           nsetup           If given, only the first NSETUP set-up atmospheres
%                            are considered, thus giving faster but less
%                            accurate result.
%------------------------------------------------------------------------------

% HISTORY: 2002.01.15  Created by Patrick Eriksson


function [L,tb_max] = qpopt_linefile(Q,f_range,lo,sideband_ratio,za_pencil,tb_limit,df_search,molecules,ldir,lformat,df_combine,nsetup)


%=== Handle optional arguments
%
if ~exist('df_combine','var')
  df_combine = [];
end
%
if ~exist('nsetup','var')
  nsetup = Inf;
end


%=== Check that line format is OK
%
if ~any( strcmp( upper(lformat), {'VERDANDI'} ) )
  error(sprintf('Unknown line format (%s)',lformat));
end


%=== Check remaining input
%
if length( f_range ) ~= 2
  error('The argument F_RANGE must have length 2.');
end
%
if ~isscalar( sideband_ratio )
  error('The argument SIDEBAND_RATIO must be a scalar.');
end
%
if ~isscalar( tb_limit )
  error('The argument TB_LIMIT must be a scalar.');
end
%
if ~isscalar( df_search )
  error('The argument DF_SEARCH must be a scalar.');
end


%=== Some frequencies
%
if ~isempty( lo )
  f_image  = sort( 2 * lo - f_range );  % Image band
  f_search = [ vec2row(f_range), vec2row(f_image) ];
  f_search = [ min(f_search)-df_search, max(f_search)+df_search ];
else
  f_search = f_range + [-df_search df_search];
end


%=== Write messages
%
out(1,1);
out(1,'Searching for transitions above given threshold ...');
out(1,sprintf('Primary band : %.3f - %.3f GHz',f_range/1e9));
if isempty( lo )
  out(1,'Image band not considered.');
else
  out(1,sprintf('Image band   : %.3f - %.3f GHz',f_image/1e9));
end
out(1,sprintf('Threshold    : %.4f K',tb_limit));
out(1,0);


%=== Convert cosmic background radiation to Tb
%
global COSMIC_BG_TEMP
%
f_cbgr  = f_range(2);  % This is the highest frerquency that can be in F_MONO
tb_cbgr = invrayjean( f_cbgr, planck( f_cbgr, COSMIC_BG_TEMP ) );


%=== Some variables for cfiles
%
tmpdir    = temporary_directory( Q.TMP_AREA );
%
nsetup = min( [ sstring_length( Q.SETUP_VMR ), nsetup ] );
%
template = which( 'basic.tmplt' );
%
QE.ABS_SAVE = 0;
QE.ABS_EXIT = 0;
%
Q.OTHER_TAGS = '"H2O","N2"';


%=== Copy continua file to temporary dir and use this as SPECTRO_DIR
%
copyfile( fullfile(Q.SPECTRO_DIR,Q.CONTINUA), fullfile(tmpdir,Q.CONTINUA) ); 
%
Q.SPECTRO_DIR = tmpdir;
%
linefile = fullfile(tmpdir,Q.LINEFILE);


%=== Use temporary dir as CALCGRIDS_DIR
%
write_datafile( fullfile(tmpdir,Q.ZA_PENCIL), za_pencil, 'VECTOR' );
%
copyfile( fullfile(Q.CALCGRIDS_DIR,Q.P_ABS), fullfile(tmpdir,Q.P_ABS) ); 
%
Q.CALCGRIDS_DIR = tmpdir;


%=== Loop isotopomers
%
L      = cell( 1000, 1 );
nout   = 0;
tb_max = zeros( 1000, 1 );
%
for imo = 1 : length( molecules )

  if strcmp( upper(lformat), 'VERDANDI' )
    tagname = vtag2arts( molecules(imo) );
    out(1,['Doing ',tagname,' (Verdandi tag ',int2str(molecules(imo)), ...
                                                                     ') ...']);
    mo  = floor( molecules(imo) / 10 );
    iso = round( rem( molecules(imo), 10 ) );
    if mo < 10
      lfile = fullfile(ldir,sprintf('v0%d.%d',mo,iso));
    else
      lfile = fullfile(ldir,sprintf('v%d.%d',mo,iso));
    end
    file_ok( lfile, Q, tmpdir );
    Ltest = read_verdandi( lfile, f_search );
  end

  out(2,sprintf(...
               'found %d transitions inside range of interest',length(Ltest)));

  %= Combine lines 
  %
  if ~isempty( df_combine )
    out(2,'combining lines');
    nin = length( Ltest );
    Ltest = combine_lines( Ltest, df_combine );
    if length(Ltest) < nin
      out(2,sprintf('new number of transitions is %d ',length(Ltest)));
    else
      out(2,'no lines combined');
    end
  end

  %= Create cfiles
  %
  out(2,'creating control files for ARTS');
  %
  Q.RETRIEVAL_TAGS = ['"',tagname,'"'];
  %
  for ii = 1 : nsetup
    Q.APRIORI_VMR    = sstring_get_i( Q.SETUP_VMR, ii );
    Q.APRIORI_PTZ    = sstring_get_i( Q.SETUP_PTZ, ii );
    [cfile,basename] = qtool( Q, tmpdir, template, QE );
    basenames{ii}    = sprintf('%s%d',basename,ii);
    cfiles{ii}       = sprintf('%s.arts',basenames{ii});
    copyfile( cfile, cfiles{ii} );
  end


  %= Wait bar
  %
  if out(2)
    hbar = waitbar(0,['Doing transitions for ',tagname]); 
  end


  %= Loop transitions
  %
  nlines = length(Ltest);
  %
  for i = 1 : nlines

    %= Write to linefile
    %
    write_linefile( linefile, Ltest(i) ,0 ,0 );

    %= Put in closest primary frequency in f_mono
    %
    if Ltest{i}.f < f_range(1)
      f_mono = f_range(1);
    elseif Ltest{i}.f > f_range(2)
      f_mono = f_range(2);
    else
      f_mono = Ltest{i}.f;
    end

    %= Add closest image frequency if frequency is closer to image band
    %= than primary band.
    %
    if ~isempty( lo ) & ...
                  min(abs(Ltest{i}.f-f_image)) < min(abs(Ltest{i}.f-f_range))
      if Ltest{i}.f < f_image(1)
        f_mono = [ f_mono, f_image(1) ];
      elseif Ltest{i}.f > f_image(2)
        f_mono = [ f_mono, f_image(2) ];
      else
        f_mono = [ f_mono, Ltest{i}.f ];
      end
      tb_fac = [1 sideband_ratio];
    else
      tb_fac = 1;
    end

    %= Create fole for F_MONO
    %
    write_datafile( fullfile(Q.CALCGRIDS_DIR,Q.F_MONO), f_mono, 'VECTOR' );

    %= Loop set-up atmospheres
    %
    saved = 0;
    %
    for ii = 1 : nsetup

      %= Run ARTS and load spectra
      %
      qp_arts( Q, cfiles{ii} );
      %
      y = read_artsvar( basenames{ii}, 'y' );

      Y = reshape( y-tb_cbgr, length( f_mono ), length( za_pencil ) );

      tb = max(max(Y').*tb_fac);

      if tb > tb_limit
        if ~saved
          nout      = nout + 1;
          L{nout}   = Ltest{i};      
          saved     = 1;
          if nargout == 1
            break;
          end
        end
        if tb > tb_max(nout)
          tb_max(nout) = tb; 
        end
      end
    end

    if saved
      out(3,sprintf('+ %.3f saved (%.3f K)',Ltest{i}.f/1e9,tb_max(nout)));
    else
      out(3,sprintf('- %.3f not saved (%.1e K)',Ltest{i}.f/1e9,tb));
    end

    if out(2)
      waitbar( i/nlines, hbar ); 
    end

  end

  delete(hbar);drawnow
  out(3,' ');

end


%=== Remove excess part of TB_MAX
%
if nout < length(L)
  L      = L(1:nout);
  tb_max = tb_max(1:nout);
end


%=== Delete the temporary directory
%
delete_tmp_dir( Q.TMP_AREA, tmpdir );


%=== Save if no output argument
%
if ~nargout
  %
  write_linefile( fullfile(Q.SPECTRO_DIR,'/',Q.LINEFILE), L , 0, 0 );
  %
end


out(1,-1);

return


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

function file_ok( lfile, Q, tmpdir )
  %
  if ~exist( lfile, 'file' )
    delete_tmp_dir( Q.TMP_AREA, tmpdir );
    error(sprintf('Could not open file %s',lfile ));
  end
return