%------------------------------------------------------------------------
% 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 )]
%
% OUT:      -
% IN:       Q              Setting structure.
% OPTIONAL: message        variable to print out a random set
%                          name in the screen 
%                        
%------------------------------------------------------------------------
%
% HISTORY: 2001.08.20  By Carlos Jimenez/ Patrick Eriksson
%
function  qp_rnd_sensor( Q , message)


global REPORT_LEVEL
r_l = Q.QP_LEVEL;

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

%=== start printing
out(1,1);
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));


%=== creating thermal noise

if Q.MEASNOISE_DO

  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' );


end

%=== creating calibration noise

if Q.CALINOISE_DO

  disp('Sorry, calibration noise not implemented yet.')

end

%
out(1,-1);