%------------------------------------------------------------------------
% NAME:     qp_rnd
%
%           A function to generate spectra corresponidng to a Q 
%           structure. The fields needed in Q to generate the
%           random spectra are stated in /Qpack/README
%           
%           NOTE: PROVISIONAL function !!, it has to be replaced
%                 by a new function generating the random inputs
%                 with mscripts, so far using arts methods 
%                 implemented by Patrick Eriksson to generate
%                 the files containing the rnd inputs.
%           
%           Extra fields needed in Q structure are
%
%    Q.RND_DIR
%      Directory where to keep the random files
%    Q.RND_NAME
%      Root name of the random files keeping the data
%    Q.RND_NUMBER
%      Number of random spectra. If there is more than one Q.RND_VMR_DIRS
%      see below) it has to be a column vector with the numbers corresponding
%      to the number of random spectra for each atmosphere e.g [100; 80] for
%      2 atmospheres.
%    Q.RND_VMR_DIRS
%      Atmospheres used as different a priori for the random calculations.
%      It can be used to train around different a priori species profiles.
%    Q.RND_PTZ
%      The corresponding ptz grids. Notice that is not a problem if they 
%      are different as ARTS makes an interpolation to a common grid,
%      given here as Q.P_ABS. Consequently the random species profiles
%      are calculated in Q.P_ABS.
%    Q.RND_SPEC_DO
%      Boolean to do species.
%    Q.RND_SPEC_TAGS
%      Species tags to do random realizations, only four at present.
%    Q.RND_OTHER_TAGS
%      Rest of species tags, these are not disturbed. 
%    Q.RND_SPEC_COVMATS
%      Covariance matrices with the variability of the species.
%    Q.RND_TEMP_DO      
%      Boolean to do temperature
%    Q.RND_TEMP_COVMATS  
%      Covariance matrices with temperature variability.
%    Q.RND_POIN_DO 
%      Boolean to do pointing, disturbing zenith angle
%    Q.RND_POIN_LIM
%      Standard deviation of the pointing variability, the same for
%      all the zenitth angles of the observation.
%    Q.RND_THN_DO
%      Boolean to do thermal noise, no correlated and identical for all channels    
%    Q.RND_THN
%      Standar deviation of the noise
%
%
% FORMAT:   qp_rnd( Q , [message)
%
% OUT:      -
% IN:       Q              Setting structure.
% OPTIONAL: message        variable to print out what random set
%                          is being done
%
% TEMPLATE: ycalc.tmplt
%------------------------------------------------------------------------
%
% HISTORY: 2001.07.118  Started by Carlos Jimenez
%
function  qp_rnd( Q , message)
%
%=== start printing
out(1,1);
out(1,['Doing ',message,' sets of spectra-observation variables']);
%
%=== checking not implemented options
if Q.RND_FREQ_DO
  error('Frequency random realizations not implemented yet, sorry');
end
if Q.RND_ALT_DO
  error('Altitude random realizations not implemented yet, sorry');
end
%
%=== Create temporary dir
%
tmpdir    = temporary_directory( Q.TMP_AREA );
%
%=== Get full path of template file and create cfile
%
template = which( 'rnd.tmplt' );
%
%=== Run ARTS for the different rnd calculations  
%
% Running as many times as different a priori
p_ap            = find(Q.RND_VMR_DIRS=='"');
p_tz            = find(Q.RND_PTZ=='"');
n_ap            = size(p_ap,2)/2;
if n_ap==0, n_ap=1; end
%
% Number and species to do realizations
p_sp            = find(Q.RND_SPEC_TAGS  =='"');
p_cov           = find(Q.RND_SPEC_COVMATS=='"');
Q.RND_SPEC_N    = size(p_sp,2)/2;
Q.RND_SPEC_FILE = [];
%
for i=1:Q.RND_SPEC_N
  is = num2str(i);
  pa = num2str(p_sp(2*(i-1)+1));
  pb = num2str(p_sp(2*(i-1)+2));
  eval(['Q.RND_SPEC_TAGS',is,'=Q.RND_SPEC_TAGS(',pa,':',pb,');'])
  pa = num2str(p_cov(2*(i-1)+1)+1);
  pb = num2str(p_cov(2*(i-1)+2)-1);
  eval(['Q.RND_SPEC_COVMAT',is,'=Q.RND_SPEC_COVMATS(',pa,':',pb,');'])
  if i==Q.RND_SPEC_N
    Q.RND_SPEC_FILE = [ Q.RND_SPEC_FILE '""'];
  else
    Q.RND_SPEC_FILE = [ Q.RND_SPEC_FILE '"",'];
  end
end
%
%=== Running for each a priori
%
% To put together the realizations from each a priori
Y =   []; T =   []; N =   []; P   = [];
SP1 = []; SP2 = []; SP3 = []; SP4 = [];
%
for i=1:n_ap
  %
  % number of realizations for this a priori
  Q.RND_N  =  Q.RND_NUMBER(1,i);
  %
  % a priori for this running
  if n_ap==1
    Q.RND_VMR_DIR=Q.RND_VMR_DIRS;
    Q.RND_PTZI=Q.RND_PTZ;
  else
    pa = num2str(p_ap(2*(i-1)+1));
    pb = num2str(p_ap(2*(i-1)+2));
    eval(['Q.RND_VMR_DIR=Q.RND_VMR_DIRS(',pa,':',pb,');'])
    pa = num2str(p_tz(2*(i-1)+1));
    pb = num2str(p_tz(2*(i-1)+2));
    eval(['Q.RND_PTZI=Q.RND_PTZ(',pa,':',pb,');'])
  end
  %
  % -- Running arts
  [cfile,basename] = qtool( Q, tmpdir, template );  
  qp_arts( Q, cfile );
  %
  %  -- Load spectra and transform with sensor matrix
  load( [Q.OUT,'.h'], '-mat' );
  y = read_artsvar( basename, 'ybatch' );
  y = h_x_h(H,y);
  Y=[Y y];
  clear H y
  %
  % -- Load species input 
  %
  if Q.RND_SPEC_DO 
    for l=1:Q.RND_SPEC_N
       ls = num2str(l);      
       eval(['TAG=Q.RND_SPEC_TAGS',ls,';'])
       TAG=TAG(2:length(TAG)-1);
       SP = read_datafile([Q.RND_DIR,Q.RND_NAME,'.',TAG,'.ab'],'MATRIX');
       eval(['SP',ls,'=[SP',ls,' SP];'])
    end
  clear SP
  end
  %
  %  Load temperature input
  if Q.RND_TEMP_DO
    t = read_datafile([Q.RND_DIR,Q.RND_NAME,'.t_abs.ab'],'MATRIX');
    T = [T t];
    clear t
  end
  %
  %  Load pointing input
  if Q.RND_POIN_DO
    p = read_datafile([Q.RND_DIR,Q.RND_NAME,'.za_pencil.ab'],'MATRIX');
    P = [P p];
    clear p
  end
  %
  % Load noise
   if Q.RND_POIN_DO
     n = read_datafile([Q.RND_DIR,Q.RND_NAME,'.noise.ab'],'MATRIX');
     N = [N n];
     clear n
   end   
  %
end
%
%=== Saving data
%
% -- saving spectra
save( [Q.RND_DIR,Q.RND_NAME,'.ybatch'], 'Y' );
%
% -- saving noise
if Q.RND_POIN_DO
  save( [Q.RND_DIR,Q.RND_NAME,'.noise'],  'N' );
end
%
% -- saving species
if Q.RND_SPEC_DO
  for l=1:Q.RND_SPEC_N
       ls = num2str(l);
       eval(['TAG=Q.RND_SPEC_TAGS',ls,';'])
       TAG=TAG(2:length(TAG)-1);
       eval(['SP=SP',ls,';'])
       save( [Q.RND_DIR,Q.RND_NAME,'.',TAG], 'SP' );
  end
end
%
% -- saving temperature 
if Q.RND_TEMP_DO
    save( [Q.RND_DIR,Q.RND_NAME,'.t_abs'], 'T' );
end
%
% -- saving pointing
if Q.RND_POIN_DO
    save( [Q.RND_DIR,Q.RND_NAME,'.za_pencil'], 'P' );
end
%
% -- saving a prori p_abs
%    NOTE: notice that p_abs is common for all the atmospheres used
%    and the retrieval vertical coordinate. The retrieval can be given
%    in altitude if temperature and H2O are provided for the measurement.
p_abs = read_artsvar( basename,'p_abs');
save( [Q.RND_DIR,Q.RND_NAME,'.p_abs'], 'p_abs' );
%
z_abs = read_artsvar( basename,'z_abs');
save( [Q.RND_DIR,Q.RND_NAME,'.z_abs'], 'z_abs' );
%
%=== Delete the temporary directory
delete_tmp_dir( Q.TMP_AREA, tmpdir );
%
out(1,-1);
%