%------------------------------------------------------------------------
% NAME:     qp_rnd_atmxsensor.m
%
%           A function to apply a sensor matrix to a rnd spectra
%
%           Use this function in a sequence
%                     qp_rnd_atm          
%                     qp_rnd_sensor         
%                     qp_rnd_atmxsensor
%
%
% FORMAT:   qp_rnd_atmxsensor( Q , [ Qi ) ]
%
% OUT:      -
% IN:       Q              rnd structure.
% OPTIONAL: Qi             inversion structure or external H flag:
%                          - if Qi is a structure, H and Hd are taken
%                          from here and not from Q, allowing different
%                          rnd sets to always be applied the same 
%                          H and Hd. H and Hd from Qi must be
%                          precalculated.
%                          - if Qi is 1 it is used as a flag to not
%                          calculate H and Hd from Q here, so they
%                          must be precalculated by qp_H or coming
%                          from another source. 
%
%------------------------------------------------------------------------
%
% HISTORY: 2001.08.04  Started by Carlos Jimenez

function  qp_rnd_atmxsensor( Q, Qi)

global REPORT_LEVEL
r_l = Q.QP_LEVEL;

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

if exist('Qi')
     % -- checking if structure or flag
     if ~isfield(Qi,'OUT')
       % -- so it's flag
       do_extH = Qi;
       Qi      = Q;
     else
       % -- it's structure, do_extH = 1 as
       %    H-Hd precalculated from Qi
       do_extH = 1;
     end
else
    % -- if no Qi obvioulsy Qi is Q
    Qi = Q;
    % -- and no external H
    do_extH = 0;
end

%=== calculate or read H matrix

if ~do_extH

  out(2,'Calculating sensor matrix');
  REPORT_LEVEL = 1;
  qp_Sx(Qi);
  qp_H(Qi);
  REPORT_LEVEL = r_l;

else

  out(2,'Reading sensor matrix');

end


%=== calculate H matrix

out(2,'Applying sensor matrix to spectra');

% -- NOTE: if limb reduction as there is a H and Hd 
%          matrix for each retrieval grid point Y 
%          becomes an structure with Y{i} the spectra
%          for each retrieval grid point

do_struc = 0;
if isfield(Q,'LRED_ON') 
   if Q.LRED_ON
     do_struc = 1;
   end
end

% -- batch spectra
load( [Q.OUT,'.ybatch'], '-mat' );
% -- reference spectra
load( [Q.OUT,'.yo'], '-mat' );

if ~do_struc

  % -- H
  load( [Qi.OUT,'.h'], '-mat' );

  % -- transforming Y
  Y  = H * Y;
  yo = H * yo';
  clear H

  % -- transforming noise
  if Q.MEASNOISE_DO
    load( [Q.OUT,'.measnoise'], '-mat' );
    load( [Qi.OUT,'.hd'], '-mat' );
    TN = Hd * TN;
    clear Hd
    Y = Y + TN;
  end
  if Q.CALINOISE_DO
    disp('Sorry, calibration option not implemented')
  end
  if Q.POLYFIT_DO > 0 & Q.POLYFIT_ORDER == 0
    load([Q.OUT,'.polyfit'], '-mat' );
    load( [Qi.OUT,'.hd'], '-mat' );
    PO = Hd * PO;
    clear Hd
    Y = Y + PO;
  end


  % -- saving
  save([Q.OUT,'.ysensor'],'Y')
  save([Q.OUT,'.yosensor'],'yo')
  clear Y yo

else

  % -- number of species and r_grid
  eval(['load ',Qi.OUT,'.h ind -mat']);

  for i=1:length(ind)

    is = num2str(i);
    for j=1:ind(i)

        js = num2str(j);
        eval(['load ',Qi.OUT,'.h H',is,js,' -mat']);
        eval(['H = H',is,js,';']) 
        Yaux{i,j}  = H * Y;
        yoaux{i,j} = H * yo';
        eval(['clear H H',is,js])

    end
  end
  
  Y = Yaux;
  yo = yoaux;
  clear Yaux yoaux
  save([Q.OUT,'.yosensor'],'yo')
  clear yo
 

  if Q.MEASNOISE_DO
    
    load( [Q.OUT,'.measnoise'], '-mat' );

    for i=1:length(ind)

      is = num2str(i);
      for j=1:ind(i)

        js = num2str(j);
        eval(['load ',Qi.OUT,'.hd Hd',is,js,' -mat']);
        eval(['Hd = Hd',is,js,';']) 
        TNaux = Hd * TN;
        Y{i,j}  = Y{i,j}  + TNaux;
        eval(['clear Hd Hd',is,js])

      end
    end
 
    clear TNaux TN

  end

  if Q.CALINOISE_DO
    disp('Sorry, calibration option not implemented')
  end

  if Q.POLYFIT_DO > 0 & Q.POLYFIT_ORDER == 0
    
    load( [Q.OUT,'.polyfit'], '-mat' );

    for i=1:length(ind)

      is = num2str(i);
      for j=1:ind(i)

        js = num2str(j);
        eval(['load ',Qi.OUT,'.hd Hd',is,js,' -mat']);
        eval(['Hd = Hd',is,js,';']) 
        POaux = Hd * PO;
        Y{i,j}  = Y{i,j}  + POaux;
        eval(['clear Hd Hd',is,js])

      end
    end
 
    clear POaux PO

  end


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


end



out(1,-1);