%------------------------------------------------------------------------
% NAME:     qcls_batch
%
%           Mscript acting as shell over qpcls to invert the
%           spectra created by qp_rnd_atmxsensor defined in
%           the structure Qr following the inversions definitions
%           given in Qo. Notice that a perfectly caharacterized
%           inversion is achieved when Qr = Qo.
%           As non linear inversions take long processing times,
%           individual inversions are saved so if the mscript
%           is halted for any reason it is possible to resume 
%           the calculation without loosing the previous inversions.
%           The output of qpcls_batch is saved in [Qo.OUT,'.saveoem'].  
%
% FORMAT:   [X,G,y_corrected,y_fitted] = qpcls_batch( Qo, Qr, [do_restart,
%                                                          do_print )]
%
% OUT:      X{i,j}         is X{j} as in qpcls for inversion i
%                          from rnd set.
%           G{i}           is G as in qpcls for inversion i   
%           y_corrected{i} is y_corrected fro inversion i    
%           y_fitted{i}    is y_fitted for inversion i
%
% IN:       Qo              structure defining OEM inversion
%           Qr              structure defining rnd set 
%
% OPTIONAL: do_restart     Flag:
%                          - 0 (default) if previous inversions exist, it
%                            reads them. 
%                          - 1 ignore previous inversions and start to
%                            invert first spectra.
%           do_print       Flag
%                          - 1 (default) plot error and individual
%                            inversions
%                          - 0 do not generate plots 
% TEMPLATE: As in qpcls 
%------------------------------------------------------------------------

% HISTORY: 2001.01.04  Created by Carlos Jimenez


function [X,G,y_corrected,y_fitted,Xiter] = qcls_batch( Qo, Qr, do_restart,do_print )


if ~exist('do_restart')
     do_restart=0;
end

if ~exist('do_print')
     do_print=1;
end


%=== previous inversions?
if ~do_restart

  aux = find(Qo.OUT=='/');
  aux = Qo.OUT(1:aux(length(aux))-1);
  addpath(aux);
  if exist([Qo.OUT,'.saveoem'],'file')==2
     load([Qo.OUT,'.saveoem'],'-mat')
     ex_inv = 1;
  else
     ex_inv = 0;
  end
  rmpath(aux);

else

     ex_inv = 0;

end

%=== load spectra
load([Qr.OUT,'.ysensor'],'-mat')


%=== prepare inversion
if ex_inv

  load([Qo.OUT,'.saveoem'],'-mat')
  i_b = size(X,1)+1;  

else

  i_b = 1;

end

i_e = size(Y,2);


if i_b<i_e

  for i= i_b:i_e

    out(1,1);
    out(1,['RETRIEVAL NUMBER ',num2str(i)]);
    out(1,-1);

    %=== inverting

    [Xi,Gi,yi_corrected,yi_fitted,Xiiter] = qpcls( Qo, Y(:,i) );

    for j=1:size(Xi,1)
       X{i,j}=Xi{j};
    end

    G{i}=Gi;
    y_corrected{i} = yi_corrected;
    y_fitted   {i} = yi_fitted;
    Xiter      {i} = Xiiter;


    %=== displaying last inversion
    if do_print
      plot_inv(Qo,Xi,i);
    end

    %=== displaying statistics
    if i>1 & do_print
      plot_sta(Qo,X,i);
    end

    save([Qo.OUT,'.saveoem'],'X','G','y_corrected','y_fitted','Xiter');
    
  end


else

    %=== displaying statistics
    if i_e>1 & do_print
      plot_sta(Qo,X,i_e);
    end

end



return


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


%-----------------------------------------------------------------------
function plot_inv(Q,X,l)
%
%   plot last inversion
%
%-----------------------------------------------------------------------


% -- getting tags
itag = sstring_length(Q.RETRIEVAL_TAGS);

for i=1:itag
   
  TAG = sstring_get_i(Q.RETRIEVAL_TAGS,i); 
  TAG = param2tag(TAG);
  load([Q.OUT,'.',TAG],'-mat')
  eval(['xt   = SP_',TAG,'(:,',num2str(l),');'])
  xf   = X{i}.xfrac;
  xap  = X{i}.apriori;
  xt   = xt ./ xap;
  z    = X{i}.z/1e3;

  figure(i)
  fig_size(10,10)
  plot(xf,z,'k',xt,z,'b--');
  legend('OEM','TRUE',0)
  lx = min([min(xf) min(xt)]);
  hx = max([max(xf) max(xt)]);
  ly = min(z);
  hy = max(z);
  axis([lx hx ly hy])
  xlabel([TAG,' profile / ',TAG,' a priori  [-]'])
  ylabel('Altitude [km]')
  aux = find(Q.OUT=='/');
  aux = Q.OUT(aux(length(aux))+1:length(Q.OUT));
  duma = find(aux=='_');
  for i=1:length(duma)
    aux(duma(i))='-';
  end
  title(['Inversion ',num2str(l),' from ',aux])


end


%-----------------------------------------------------------------------
function plot_sta(Q,X,l)
%
%   plot statistics until last inversion
%
%-----------------------------------------------------------------------


% -- getting tags
itag = sstring_length(Q.RETRIEVAL_TAGS);

for i=1:itag
   
  TAG = sstring_get_i(Q.RETRIEVAL_TAGS,i); 
  TAG = param2tag(TAG);
  load([Q.OUT,'.',TAG],'-mat')

  z    = X{1,1}.z/1e3;
  sp   = zeros(length(z),l);

  for j=1:l

    eval(['xt   = SP_',TAG,'(:,',num2str(j),');'])
    xap  = X{j,i}.apriori;
    xt   = xt ./ xap;
    sp(:,j) = X{j,i}.xfrac-xt;

  end

  m_z = 100*mean(sp');
  s_z = 100*std(sp');


  figure(gcf+1)
  fig_size(10,10)

  subplot(1,2,1)
  plot(m_z,z,'k');
  lx = min(m_z);
  hx = max(m_z);
  ly = min(z);
  hy = max(z);
  axis([lx hx ly hy])
  xlabel('Bias [%]')
  ylabel('Altitude [km]')

  subplot(1,2,2)
  plot(s_z,z,'k');
  lx = min(s_z);
  hx = max(s_z);
  axis([lx hx ly hy])
  xlabel('Std [%]')

  aux = find(Q.OUT=='/');
  aux = Q.OUT(aux(length(aux))+1:length(Q.OUT));
  duma = find(aux=='_');
  for i=1:length(duma)
    aux(duma(i))='-';
  end
  title([TAG,' retrieval error from ',aux,'                                          '])
  plot1x2

end

%-----------------------------------------------------------------
function TAG = param2tag(TAG);
%----------------------------------------------------------------

    if      strcmp(TAG(1:2),'H2') 
      if      strcmp(TAG(1:3),'H2O') 
       TAG = 'H2O';
      end
    elseif strcmp(TAG(1:2),'O2') 
       TAG = 'O2';
    elseif strcmp(TAG(1:2),'O2') 
       TAG = 'N2';
    elseif strcmp(TAG(1:2),'CO') 
      if strcmp(TAG(1:3),'CO2') 
       TAG = 'CO2';
      end
    end 

return


%--------------------------------------------------------------------
function plot1x2
%--------------------------------------------------------------------
%
a       = 0.02;
w	= 0.80;
y	= 0.10;
%
h	= 0.42;
x1	= 0.10;
x2	= 1*(h+a)+x1;
%
subplot(1,2,1)
set(gca,'Pos',[x1 y h w]);
%
subplot(1,2,2)
set(gca,'Pos',[x2 y h w]);
set(gca,'YTickLabel','');
%
return
%