%------------------------------------------------------------------------
% 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.
%           
%           Use this function in a sequence
%                     qp_rnd_atm          
%                     qp_rnd_sensor         
%                     qp_rnd_atmxsensor
% 
%           The rnd realizations are stored in files looking like:  
%            
%            out.ClO.ab     example of rnd realizations of species
%                           saved in readable format for ARTS batch
%                           calculations.  
%            out.ClO        as above but saved in standard matlab format.
%            
%            out.t_abs.ab   if Q.DO_TEMPERATURE > 0 rnd realizations of 
%                           temperature saved in readable format for ARTS 
%                           batch calculations.  
%            out.t_field    as above but saved in standard matlab format.
%
%            out.z_abs.ab   if Q.DO_TEMPERATURE > 0 and Q.HSE_IN_ON = 1
%                           corresponding realizations of altitude
%                           fullfilling hse saved in readable format
%                           for ARTS batch calculations.   
%            out.z_field    as above but saved in standard matlab format.
%                           
%            out.za_pencil.ab  if Q.POINTING_DO > 0 realizations of zenith
%                           angle in readable format for ARTS batch
%                           calculations
%            out.za_field   as above but saved in standard matlab format.
%             
%            
%            out.p_abs                           
%            out.t_abs                           
%            out.z_abs
%            out.za_pencil                           
%            out.vmrs       a priori pressure, temperature, altitude, pencil
%                           beam zenith angle and species vmr values in 
%                           standard matlab format.
%                           Notice the following:
%                             - if Q.DO_TEMPERATURE > 0 the temperature stored 
%                               is the mean value of all the temperature
%                               realizations, i.e. the mean value of 
%                               out.t_field. This might be different from
%                               the apriori temperature defined in Q if e.g.
%                               different SETUP atmosphers are used.
%                             - if Q.DO_TEMPERATURE > 0 and Q.HSE_IN_ON = 1
%                               the same applies to the altitudes, i.e it
%                               is the mean value of out.z_field.
%                             - if Q.POINTING_DO > 0 the same applies to
%                               the zenith angles, i.e. it is the 
%                               mean value of out.za_field.
%                             - the same also applies to vmrs, if different
%                               SETUP atmospheres are used the mean value
%                               of the species realizations is stored, which
%                               might be different from the a priori in Q.
%
%            out.f_mono     monochromatic frequencies in standard matlab
%                           format.
%
%            out.yo         reference spectra corresponding to apriori
%                           state in standard matlab format. Notice          
%                           that it might be different from the mean
%                           spectrum of out.ybatch if the SETUP states
%                           are not states around the APRIORI, but is the
%                           right *linearization state* as the Kx and
%                           consequently the Hd will be calculated around
%                           this APRIORI. 
%            out.ybatch     radiance realizations in standard matlab format.
%
%
% FORMAT:   qp_rnd( Q , [message , Qr)]
%
% OUT:      -
% IN:       Q              Setting structure.
% OPTIONAL: message        variable to print out a random set
%                          name in the screen 
%           Qr             If Qr is given the rnd realizations are not done
%                          but copied from Qr. This can be used e.g. to share
%                          rnd realizations between Q and Qr when Q and
%                          Qr have different sensor.             
%
% TEMPLATE: rnd.tmplt
%------------------------------------------------------------------------
%
% HISTORY: 2001.08.15  By Carlos Jimenez/ Patrick Eriksson
%
function  qp_rnd_atm( Q , message, Qr)


global REPORT_LEVEL
global EARTH_RADIUS
r_l = Q.QP_LEVEL;

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

if ~exist('Qr')
   do_qr = 0;;
else
   do_qr = 1;
end

%=== start printing
out(1,1);

if ~do_qr

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,'"'];

elseif isempty(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));

% -------   running around apriori state before doing realizations
%           to obtain reference spectra. 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 the Kx and
%           consequently the Hd were calculated around
%           this APRIORI. So care has to be taken by the 
%           user to assure this
yo     = qp_ympb( Q );
yo     = yo'; 
save( [Q.OUT,'.yo'], 'yo' );
clear yo

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(fullfile( 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);
    if Q.HSE_IN_ON
      z_grid  = PTZ(:,3);
      z_grid  = interpp(p_grid,z_grid,p_abs);    
    end

  end
  % -- always as needed by qp_rnd_noise 
  za = read_datafile(fullfile( 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

    TAG = sstring_get_i(Q.RETRIEVAL_TAGS,i); 
    TAGn = TAG;   
 
    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,');'])

    % -- if parameterizations as retrieval tags, reading the 
    %    corresponding profiles
    %      NOTE: 2 tag names needed as ybatch in arts needs to see
    %            the main molecule O3 while for the rnd operations
    %            we need to add the isotope to save the realizations
    %            for different isotopes. The same scheme is not applied
    %            to freq range as the far from the band won't be part
    %            of the realizations.
  
    TAGarts = param2arts(TAG);
    TAG     = param2tag(TAG);

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

    % -- getting covariance matrix

    REPORT_LEVEL = 0;
    S = sFromFile(fullfile( Q.RETRIEVDEF_DIR, COVMAT ), log10(p_abs) );
    REPORT_LEVEL = r_l;

    % -- generating realizations

    if ( length(Q.CLS_SPECIES_POS_ON)==1 & Q.CLS_SPECIES_POS_ON )  |  ...
       ( length(Q.CLS_SPECIES_POS_ON)>1 & Q.CLS_SPECIES_POS_ON(i) )

      %-- @ lognormal pdf --> x = log(@) is normal
      %   then the S read here is assumed to be Sx defined
      %   in the normal space, created e.g. by sForLogNorm
      %   from a given S@.
 
      %-- obtaining mu for the given sigma in Sx,
      %   these are the mu and sigma of the N(mu,sigma).
      %   Then exp(N(mu,sigma)) gives the
      %   lognorm(mu,sigma) of statistics S@
      %
      mu    = log(1) - 0.5 * (diag(S));    

      %-- N(0,sigma) by Chloleski decomp method 
      %    
      R =  chol(S)' * randn(n_pg, Q.NUMBER_DO_B); 

      %--adding mu to have N(mu,sigma)
      %
      mu = repmat(mu,1,Q.NUMBER_DO_B);
      R = mu + R;

      %-- obtaining lognorm(mu,sigma)
      %   as exp(N(mu,sigma)
      R  = exp(R);

      %-- refining so mean is really one (1)
      %
      if size(R,2) == 1
        Rm = R;
      else
        Rm = (mean(R'))';
      end
      Rm = repmat(Rm,1,Q.NUMBER_DO_B);
      R = 1 + R - Rm;
      clear Rm

      % -- un-normalizing
      sp_ap = repmat(sp_ap,1,Q.NUMBER_DO_B);     
      R     = sp_ap .* R;

      % -- putting negative values to zero
      %    NOTE: even using lognorm, as we guarantied
      %          that the mean is the true mean (1)
      %          for a small number of realizations
      %          negative values can be obtained.
      Rm = reshape(R,1,size(R,1)*size(R,2));
      Rm(find(Rm<0)) = 1e-24;
      R = reshape(Rm,size(R,1),size(R,2)); 
      clear Rm

    else
 
      %-- normal pdf
      %   NOTE: using Chloleski decomp method 
      %         with relative std and mean 0
      R =  chol(S)' * randn(n_pg, Q.NUMBER_DO_B); 

      % -- removing mean value to make real mean 0
      %    adding one to have mean 1 and realtive std
      if size(R,2) == 1
        Rm = R;
      else
        Rm = (mean(R'))';
      end
      Rm = repmat(Rm,1,Q.NUMBER_DO_B);
      R = 1 + R - Rm;
      clear Rm

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

    end

    % -- saving species
    if j==1
      eval(['SP_',TAG,'=R;'])
    else
      eval(['SP_',TAG,'= [ SP_',TAG,' R];'])
    end


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

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

    % -- creating filename for keeping rnd data
    % fname = [Q.OUT,'.', TAG,'.ab'];

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

  % -- for do_z depending in temperature and hse
  Q.ALTITUDE_DO = 0;

  if Q.TEMPERATURE_DO

    % -- getting covariance matrix
    REPORT_LEVEL = 0;
    S = sFromFile(fullfile( 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');

    %=== Writing altitude realizations to fulfill hse


    if Q.HSE_IN_ON
 
      Q.ALTITUDE_DO = 1;

      % -- getting h2o profile

      PSP     = read_datafile([Q.APRIORI_VMR_B,'.H2O.aa'],'MATRIX',1); 
      h2o_abs = PSP(:,2);
      h2o_p   = PSP(:,1);
      h2o_abs = interpp(h2o_p,h2o_abs,p_abs); 

      A = zeros(length(p_abs),Q.NUMBER_DO_B);

      for a=1:Q.NUMBER_DO_B
         A(:,a)= hseCalc(p_abs,R(:,a),z_grid,h2o_abs,Q.HSE_PREF,Q.HSE_ZREF,EARTH_RADIUS,1);
      end

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

      % -- saving altitude
      if j==1
        Z = A;
      else
        Z = [ Z A ];
      end

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

  end

  %=== Writing pointing realizations

  if Q.POINTING_DO

    n_za    = length(za);

    % -- generating realizations

    if Q.POINTING_COR

      if strcmp(Q.POINTING_PDF,'gaussian')
        R = ones(n_za,1)*(Q.POINTING_STDV * randn(Q.NUMBER_DO_B, 1))'; 
      elseif strcmp(Q.POINTING_PDF,'uniform')
        R = ones(n_za,1)*(Q.POINTING_STDV * (rand(Q.NUMBER_DO_B, 1)-0.5))'; 
      else
        error('Not recognized Q.POINTING_PDF');      
      end

    else
      if strcmp(Q.POINTING_PDF,'gaussian')
        R = (Q.POINTING_STDV * randn(Q.NUMBER_DO_B, n_za))'; 
      elseif strcmp(Q.POINTING_PDF,'uniform')
        R = (Q.POINTING_STDV * (rand(Q.NUMBER_DO_B, n_za)-0.5))'; 
      else
        error('Not recognized Q.POINTING_PDF');      
      end
    end


    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' );  % y before
 

  % -- saving intermediate spectra to be put together 
  %    after spectrum in qp_rnd_atmxsensor
  save( [Q.OUT,num2str(j),'.ybatch'], 'y' );
  clear Y

end


%=== Saving data

out(2,'Saving spectra');

% -- saving spectra
Y = [];
for j=1:n_ap
  load( [Q.OUT,num2str(j),'.ybatch'], '-mat' );
  Y =[Y y];
  clear y
  delete([Q.OUT,num2str(j),'.ybatch'])  
end


save( [Q.OUT,'.ybatch'], 'Y' );
clear Y

% -- saving species
for i=1:n_sp

  TAG = sstring_get_i(Q.RETRIEVAL_TAGS,i); 
  TAG = param2tag(TAG);
  eval(['save ',Q.OUT,'.',TAG,' SP_',TAG,' -mat'] );

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

end

% -- saving temperature and z if hse
if Q.TEMPERATURE_DO
  % NOTE: called t_field to distinguish from t_abs
  save( [Q.OUT,'.t_field'], 'T' );    
  if Q.HSE_IN_ON
    save( [Q.OUT,'.z_field'], 'Z' );         
  end
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, the same aplies to z_abs.


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
  if Q.HSE_IN_ON
    z_abs = mean(Z');
  else
    z_abs = read_artsvar( basename,'z_abs');
  end
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 );
%


else % do_qr



out(1,['Copying ',message,' sets of spectra-observation variables']);


% -- saving original ybatch
%    
load([Qr.OUT,'.ybatch'], '-mat' );
save( [Q.OUT,'.ybatch'], 'Y' );
clear Y
load([Qr.OUT,'.yo'], '-mat' );
save( [Q.OUT,'.yo'], 'yo' );
clear yo

% -- saving species

p_sp            = find(Qr.RETRIEVAL_TAGS  =='"');
n_sp            = size(p_sp,2)/2;
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 = Qr.RETRIEVAL_TAGS(',pa,':',pb,');'])
  TAG = param2arts(TAG);
  load([Qr.OUT,'.',TAG], '-mat' );
  eval(['save ',Q.OUT,'.',TAG,' SP_',TAG]);

end


% -- saving temperature
if Q.TEMPERATURE_DO
  load([Qr.OUT,'.t_field'], '-mat' );
  save([Q.OUT,'.t_field'], 'T' );   
  clear T 
end

% -- saving pointing
if Q.POINTING_DO
  load([Qr.OUT,'.za_field'], '-mat' );
  save( [Q.OUT,'.za_field'], 'P' );
  clear 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

load([Qr.OUT,'.t_abs'], '-mat' );
save( [Q.OUT,'.t_abs'], 't_abs' );

load([Qr.OUT,'.z_abs'], '-mat' );
save( [Q.OUT,'.z_abs'], 'z_abs' );

load([Qr.OUT,'.vmrs'], '-mat' );
save( [Q.OUT,'.vmrs'], 'vmrs' );


load([Qr.OUT,'.p_abs'], '-mat' );
save( [Q.OUT,'.p_abs'], 'p_abs' );

load([Qr.OUT,'.za_pencil'], '-mat' );
save( [Q.OUT,'.za_pencil'], 'za_pencil' );

load([Qr.OUT,'.f_mono'], '-mat' );
save( [Q.OUT,'.f_mono'], 'f_mono' );

end %do_qr

out(1,-1);

return

%------------------------------------------------------------------
%
%                             SUBFUNCTIONS
%  
%-----------------------------------------------------------------

function TAG = param2tag(TAG);

isbar =find(TAG == ',');

%== leaving 1st species

if ~isempty(isbar)

 %== checking for isotopes
 TAG = TAG( 1:( isbar(1)-1 ));

end

%== leaving isotope info but removing freq range
%   or leaving continuum, both cases replacing 
%   '-' by '_'

isbar =find(TAG == '-');


if ~isempty(isbar)

  if length(isbar) > 1  % isotope with frequencies
     
    %=  remove freq range
    TAG = TAG(1:( isbar(2)-1) );       

  end
  TAG(isbar(1)) = '_'; 

end

return

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

function TAG = param2arts(TAG);

isbar =find(TAG == ',');

if ~isempty(isbar)

 %== checking for isotopes
 TAG = TAG( 1:( isbar(1)-1 ));

else

  isbar =find(TAG == '-');

  if ~isempty(isbar)

    %== checking for isotopes
    TAG = TAG(1:( isbar(1)-1) );

  end

end

return