%------------------------------------------------------------------------
% NAME:     qp_rnd_atm
%
%           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.  Notice
%           that in the calculation the sensor matrix is not
%           applied, i.e the function generates monochromatic
%           pencil beam spectra.
%           
%
%
% FORMAT:   qp_rnd( Q , [message )]
%
% OUT:      -
% IN:       Q              Setting structure.
% OPTIONAL: message        variable to print out a random set
%                          name in the screen 
%                        
%
% TEMPLATE: rnd.tmplt
%------------------------------------------------------------------------
%
% HISTORY: 2001.08.15  By Carlos Jimenez/ Patrick Eriksson
%
function  qp_rnd_atm( Q , message)


global REPORT_LEVEL
r_l = Q.QP_LEVEL;

%=== checking inputs
if ~exist('message')
   message='';
end

%=== start printing
out(1,1);
out(1,['Doing ',message,' sets of spectra-observation variables']);


%=== checking not implemented options
if Q.CONTABS_DO
  disp('Warning: continuum random realizations not implemented, the will not be done.')
end

%=== if a priori as only atmosphere for rnd realizations 
if ~isfield(Q,'SETUP_VMR')

     Q.SETUP_VMR = ['"',Q.APRIORI_VMR,'"'];
     Q.SETUP_PTZ = ['"',Q.APRIORI_PTZ,'"'];

end

%=== creating rnd realizations of _DO variables


% -- number of different a priori
p_ap            = find(Q.SETUP_VMR=='"');
p_tz            = find(Q.SETUP_PTZ=='"');
n_ap            = size(p_ap,2)/2;

% -- checking right number
if length(Q.NUMBER_DO)~=n_ap
  error('Q.NUMBER_DO do not have as many elements as number of SETUP atmospheres');
end

% -- number and species to do realizations
p_sp            = find(Q.RETRIEVAL_TAGS  =='"');
p_cov           = find(Q.SPECIES_COVMATS =='"');
n_sp            = size(p_sp,2)/2;


% -- setting seed of random generator
rand('state',sum(100*clock));


for j=1:n_ap

  % -- a priori 
  pa = num2str(1  + p_ap(2*(j-1)+1));
  pb = num2str(-1 + p_ap(2*(j-1)+2));
  eval(['Q.APRIORI_VMR_B=Q.SETUP_VMR(',pa,':',pb,');'])
  pa = num2str(1  + p_tz(2*(j-1)+1));
  pb = num2str(-1 + p_tz(2*(j-1)+2));
  eval(['Q.APRIORI_PTZ_B=Q.SETUP_PTZ(',pa,':',pb,');'])


  out(2,['Starting rnd realizations for a priori number ',num2str(j)]);

  % -- p, t and za grid
  PTZ     = read_datafile(Q.APRIORI_PTZ_B,'MATRIX',1); 
  p_grid  = PTZ(:,1);
  p_abs   = read_datafile([Q.CALCGRIDS_DIR,'/',Q.P_ABS],'VECTOR'); 
  n_pg    = length(p_abs);
  if Q.TEMPERATURE_DO
    t_grid  = PTZ(:,2);
    t_grid  = interpp(p_grid,t_grid,p_abs);
  end
  % -- always as needed by qp_rnd_noise 
  za = read_datafile([Q.CALCGRIDS_DIR,'/',Q.ZA_PENCIL],'VECTOR');
  
  % -- number of realizations for this a priori
  Q.NUMBER_DO_B  =  Q.NUMBER_DO(1,j);


  %=== Writing species realizations

  for i=1:n_sp

    is = num2str(i);
    pa = num2str( 1 + p_sp(2*(i-1)+1));
    pb = num2str(-1 + p_sp(2*(i-1)+2));
    eval(['TAG = Q.RETRIEVAL_TAGS(',pa,':',pb,');'])

    pa = num2str( 1 + p_cov(2*(i-1)+1));
    pb = num2str(-1 + p_cov(2*(i-1)+2));
    eval(['COVMAT = Q.SPECIES_COVMATS(',pa,':',pb,');'])

    % -- obtaining a priori and interpolating to p_abs
    PSP     = read_datafile([Q.APRIORI_VMR_B,'.',TAG,'.am'],'MATRIX',1); 
    sp_ap   = PSP(:,2);
    sp_ap   = interpp(p_grid,sp_ap,p_abs);


    % -- getting covariance matrix
    REPORT_LEVEL = 0;
    S = sFromFile([Q.RETRIEVDEF_DIR,'/',COVMAT],log10(p_abs));
    REPORT_LEVEL = r_l;

    % -- generating realizations

    R = ( (chol(S)') * randn(n_pg, Q.NUMBER_DO_B)); 

    % -- removing mean value so mean of set is zero

    Rm = (mean(R'))';
    Rm = repmat(Rm,1,Q.NUMBER_DO_B);
    R = 1 + R -Rm;
    clear Rm

    for l=1:Q.NUMBER_DO_B
      R(:,l) = sp_ap .* R(:,l);
    end

    % -- saving species
    if j==1
      eval(['SP_',TAG,'=R;'])
    else
      eval(['SP_',TAG,'= [ SP_',TAG,' R];'])
    end
 
    % -- creating filename
    fname = [Q.OUT,'.', TAG,'.ab'];
    out(2,['  Creating ',num2str(Q.NUMBER_DO_B),' profiles for ',TAG]);

    % -- saving realizations
    write_datafile(fname,R,'MATRIX');

    % -- saving Q:FILE
    if i==1
      Q.FILE_TAGS = '""'; 
    else
      Q.FILE_TAGS = [ Q.FILE_TAGS ',""' ];
    end
     
  end

  %=== Writing temperature realizations

  if Q.TEMPERATURE_DO

    % -- getting covariance matrix
    REPORT_LEVEL = 0;
    S = sFromFile([Q.RETRIEVDEF_DIR,'/',Q.TEMPERATURE_COVMAT],log10(p_abs));
    REPORT_LEVEL = r_l;


    % -- generating realizations 
    R = ((chol(S)') * randn(n_pg,Q.NUMBER_DO_B)); 
    for i=1:Q.NUMBER_DO_B
      R(:,i) = t_grid +  R(:,i);
    end

    % -- creating filename
    fname = [Q.OUT,'.t_abs.ab'];
    out(2,['  Creating ',num2str(Q.NUMBER_DO_B),' profiles for temperature']);

    % -- saving temperature
    if j==1
      T = R;
    else
      T = [ T R ];
    end

    % -- saving realizations
    write_datafile(fname,R,'MATRIX');
 
  end

  %=== Writing pointing realizations

  if Q.POINTING_DO

    n_za    = length(za);

    % -- generating realizations
 
    R = (Q.POINTING_STDV * randn(Q.NUMBER_DO_B, n_za))'; 
    for i=1:Q.NUMBER_DO_B
      R(:,i) = za +  R(:,i);
    end

    % -- creating filename
    fname = [Q.OUT,'.za_pencil.ab'];
    out(2,['  Creating ',num2str(Q.NUMBER_DO_B),' zenith angles scans.']);

    % -- saving pointing
    if j==1
      P =R;
    else
      P = [ P R ];
    end

    % -- saving realizations only offset
    R  =  R(1,:) - za(1);

    write_datafile(fname,R,'MATRIX');
 
  end


  %=== Running arts to get spectra

  out(2,'Generating spectra');
  tmpdir    = temporary_directory( Q.TMP_AREA );
  template = which( 'rnd.tmplt' );
  [cfile,basename] = qtool( Q, tmpdir, template );  
  qp_arts( Q, cfile );

  %  -- Load spectra 
  y = read_artsvar( basename, 'ybatch' );

  % -- saving spectra
  if j==1
    Y =y;
  else
    Y = [ Y y ];
  end
  clear y

end

%=== Saving data

out(2,'Saving spectra');

% -- saving spectra
save( [Q.OUT,'.ybatch'], 'Y' );
clear Y

% -- NOTE:  The difference to the reference spectra is
%           also saved. The reference spectra
%           corresponds to the APRIORI state, notice
%           that might not be the close to the mean
%           spectrum if the SETUP states are not
%           states around the APRIORI, but is the
%           right *linearization state* as th Kx are 
%           calculates around this APRIORI.

yo = (read_artsvar( basename, 'y' ))';
save( [Q.OUT,'.yo'], 'yo' );

% -- saving species
for i=1:n_sp

  is = num2str(i);
  pa = num2str( 1 + p_sp(2*(i-1)+1));
  pb = num2str(-1 + p_sp(2*(i-1)+2));
  eval(['TAG = Q.RETRIEVAL_TAGS(',pa,':',pb,');'])
  eval(['save ',Q.OUT,'.',TAG,' SP_',TAG,' -mat'] );

  % -- saving mean starte if more than one a priori
  if n_ap > 1
    eval(['vmrs{',is,'} = mean(SP_',TAG,''');']);
  end

end

% -- saving temperature
if Q.TEMPERATURE_DO
  % NOTE: called t_field to distinguish from t_abs
  save( [Q.OUT,'.t_field'], 'T' );    
end

% -- saving pointing
if Q.POINTING_DO
  save( [Q.OUT,'.za_field'], 'P' );
end

% -- saving a prori p_abs t_abs z_abs and vmrs
%    NOTE: notice that p_abs is common for all the atmospheres used.
%    t_abs is only common if temp was not disturbed, if it was disturbed
%    the mean temperature of T is saved. z_abs is obtained fullfilling
%    hydrostatique equlibrium. vmrs is also the mean state if more than
%    one apriori is used


if ~Q.TEMPERATURE_DO
  t_abs = read_artsvar( basename,'t_abs');
else
  t_abs = mean(T');
end
save( [Q.OUT,'.t_abs'], 't_abs' );

if ~Q.TEMPERATURE_DO
  z_abs = read_artsvar( basename,'z_abs');
else
  %!!! to think 
  z_abs = read_artsvar( basename,'z_abs');
end
save( [Q.OUT,'.z_abs'], 'z_abs' );


% -- if more that one apriori, already mean state calculated
if n_ap==1
  vmrs = read_artsvar( basename,'vmrs');
end  
save( [Q.OUT,'.vmrs'], 'vmrs' );


% - save also really needed
save( [Q.OUT,'.p_abs'], 'p_abs' );
za_pencil = za;
save( [Q.OUT,'.za_pencil'], 'za_pencil' );
f_mono = read_artsvar( basename, 'f_mono' );
save( [Q.OUT,'.f_mono'], 'f_mono' );


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