%------------------------------------------------------------------------
% NAME:     qp_rnd_sensor
%
%           A function to generate the rand fields associated with a Q 
%           structure defining the sensor. The fields needed in Q to 
%           generate the rand fields are stated in /Qpack/README.
%           So far only measurement noise.
%
%
% FORMAT:   qp_rnd_sensor( 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.   
%                        
%------------------------------------------------------------------------
%
% HISTORY: 2001.08.20  By Carlos Jimenez/ Patrick Eriksson
%
function  qp_rnd_sensor( Q , message, Qr)


global REPORT_LEVEL
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


  %=== creating thermal noise

  out(1,['Generating ',message,' rnd noise for sensor']);


  %=== loading sensor variables

  out(2,'Getting sensor data');

  p_abs       = read_datafile([Q.CALCGRIDS_DIR,'/',Q.P_ABS],    'VECTOR'); 
  f_mono      = read_datafile([Q.CALCGRIDS_DIR,'/',Q.F_MONO],   'VECTOR'); 
  za_p        = read_datafile([Q.CALCGRIDS_DIR,'/',Q.ZA_PENCIL],'VECTOR'); 

  REPORT_LEVEL = 0;
  [H,f_y,za_y,f_sensor,za_sensor] = ...
                     qp_H(Q,p_abs,f_mono,za_p,0);
  REPORT_LEVEL = r_l; 


  %=== Number of realizations
  n_real  =  sum(Q.NUMBER_DO);

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

  out(2,'Doing measurement noise');

  % -- getting covariance matrix
  % -- NOTE: Doing each za at one time to avoid a covariance
  %          matrix of dimension f_y x f_y, even play with sparse
  %          matrices lis very slow

  REPORT_LEVEL = 0;
  S = sFromFile([Q.SENSOR_DIR,'/',Q.MEASNOISE_COVMAT],f_sensor);
  REPORT_LEVEL = r_l;

  % -- generating realizations

  n_f   = length(f_sensor);
  n_za  = length(f_y)/n_f;


  for i= 1:n_za

    out(2,['  Zenith angle ',num2str(i)]); 
    r = ((chol(S))' * randn(n_f,n_real)); 
    
    rm = (mean(r'))';
    rm = repmat(rm,1,size(r,2));
    r = rm-r;

    if i==1
      TN = r;
    else
      TN = [TN;r]; 
    end

  end

  % -- saving thermal noise
  save( [Q.OUT,'.measnoise'], 'TN' );
  clear TN


  %=== creating constant baseline offset

  if Q.POLYFIT_DO > 0

    if Q.POLYFIT_ORDER == 0

      out(1,['Generating ',message,' baseline offset order 1 for sensor']);

      % -- getting covariance matrix

      REPORT_LEVEL = 0;
      S = sFromFile([Q.SENSOR_DIR,'/',sstring_get_i(Q.POLYFIT_COVMATS,1)],za_sensor);
      REPORT_LEVEL = r_l;

      n_za = length(za_sensor);

      % -- random values per za

      ra   = (chol(S))' * randn( n_za, n_real); 

      % -- extending to whole freq vector
      
      PO    = zeros( n_za * n_f, n_real ); 
      for j = 1:n_za
   
        out(2,['  Zenith angle ',num2str(j)]); 
        ind      = (1 + ( j -1 ) * n_f ) : ( j * n_f);
        PO(ind,:) = repmat( ra(j,:), n_f, 1);        

      end

      % -- saving polyfit
      save( [Q.OUT,'.polyfit'], 'PO' );
      
    else
  
      disp('Sorry, baseline offsets of order > 1 not implemented yet.')

    end

  end


else %do_qr

  out(1,['Copying ',message,' rnd noise for sensor']);

  load([Qr.OUT,'.measnoise'], '-mat' );
  save( [Q.OUT,'.measnoise'], 'TN' );
  clear TN

  if Q.POLYFIT_DO > 0 & Q.POLYFIT_ORDER == 1

    load([Qr.OUT,'.polyfit'], '-mat' );
    save( [Q.OUT,'.polyfit'], 'PO' );
    clear PO

  end

end %do_qr


out(1,-1);