%------------------------------------------------------------------------ % NAME: qp_H % % Sets up H and Hd. The Hd matrix is treated in a seperate % function (qp_Hd). % % See the README file for options. % % If the function is called without output arguments, H and Hd % are stored to files. With output arguments, no files are % created. % % The calculations can be done for arbitrary P_ABS, F_MONO and % ZA_PENCIL. If these variables not are given, they are read from % Q.OUT. % % FORMAT: q_H( Q ) % or % [H,f_y,za_y,f_sensor,za_sensor] = % qp_H( Q [,p_abs,f_mono,za_pencil] ) % % OUT: (see above) % IN: - % OPTIONAL: Q The setting structure. % p_abs % f_mono % za_pencil % % TEMPLATE: - %------------------------------------------------------------------------ % HISTORY: 2001.03.25 Created by Patrick Eriksson. % 2001.09.02 Updated by Carlos Jimenez function [H,f_y,za_y,f_sensor,za_sensor] = ... qp_H(Q,p_abs,f_mono,za_pencil) %=== Shall results be saved % do_save = ~nargout; %=== Read f_mono and za_pencil if not input variables if nargin == 1 f_mono = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.F_MONO ), 'VECTOR' ); za_pencil = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.ZA_PENCIL ), ... 'VECTOR' ); end %=== Init the H matrices [H,f_y,za_y,Hd,f_sensor,za_sensor] = hInit( f_mono, za_pencil ); %=== Antenna pattern if Q.ANTENNA_ON % za_obs = read_datafile( fullfile( Q.SENSOR_DIR, Q.ANTENNA_ZA ), 'VECTOR' ); % fname = fullfile( Q.SENSOR_DIR, Q.ANTENNA_FILE ); % move = any( Q.ANTENNA_MOVE > 0 ); % [H,f_y,za_y,za_sensor] = hAntennaFromFileAdv( H, f_sensor, za_sensor, ... za_obs, fname, Q.ANTENNA_ORDER, Q.ZA_ORDER, 0, 0, move, Q.ANTENNA_MOVE ); end %=== Some stuff that is independent of frequency switch % if Q.BACKEND_ON & Q.BACKEND_IF if Q.DSB_FPRIMARY <= Q.DSB_LO error('With Q.BACKEND_IF=1, DSB_FPRIMARY > DSB_LO is required.') end end if ~Q.FREQSWITCH_ON % Remember to implement any change also for the frequency switch part %=== Mixer and sideband filter if Q.DSB_ON % fname = fullfile( Q.SENSOR_DIR, Q.DSB_FILE ); % [H,f_y,za_y,f_sensor] = hMixerFromFile( H, f_sensor, za_sensor, ... Q.DSB_LO, Q.DSB_FPRIMARY, fname, Q.DSB_ORDER, Q.F_ORDER ); end %=== Backend if Q.BACKEND_ON % f_obs = read_datafile( fullfile(Q.SENSOR_DIR,Q.BACKEND_FREQS), 'VECTOR' ); % if Q.BACKEND_IF f_sensor = f_sensor - Q.DSB_LO; end % fname = fullfile( Q.SENSOR_DIR, Q.BACKEND_FILE ); % [H,f_y,za_y,f_sensor] = hBackendFromFile( H, f_sensor, za_sensor, ... f_obs, fname, Q.BACKEND_ORDER, Q.F_ORDER); end %=== Frequency switch else if ~Q.BACKEND_ON error('Backend must be on when doing frequency switch.'); end if ~Q.BACKEND_IF error('Backend frequencies must be in IF when doing frequency switch.'); end if length( za_sensor ) > 1 error('Frequency switch is only implemented for 1 zenith angle.'); end %= Mixer part % lo1 = Q.DSB_LO + Q.FREQSWITCH_THROW/2; lo2 = Q.DSB_LO - Q.FREQSWITCH_THROW/2; % if Q.DSB_ON % fname = fullfile( Q.SENSOR_DIR, Q.DSB_FILE ); [H1,f_y,za_y,fsensor1] = hMixerFromFile( 1, f_sensor, za_sensor, ... lo1, Q.DSB_FPRIMARY, fname, Q.DSB_ORDER, Q.F_ORDER ); [H2,f_y,za_y,fsensor2] = hMixerFromFile( 1, f_sensor, za_sensor, ... lo2, Q.DSB_FPRIMARY, fname, Q.DSB_ORDER, Q.F_ORDER ); fsensor1 = fsensor1 - lo1; fsensor2 = fsensor2 - lo2; else H1 = eye( length( f_sensor ) ); H2 = H1; fsensor1 = f_sensor - lo1; fsensor2 = f_sensor - lo2; end %= Position of backend channels % f_obs = read_datafile( fullfile(Q.SENSOR_DIR,Q.BACKEND_FREQS), 'VECTOR' ); df = f_obs(2) - f_obs(1); ch_shift = round( Q.FREQSWITCH_THROW / df ); % if any( diff(f_obs) ~= df ) error('The backend channels must be equally spaced.'); end if abs(rem(Q.FREQSWITCH_THROW,df)) > 1e-5 error('The frequency throw is not a multiple of the channel spacing'); end if isodd( ch_shift ) error('The frequency throw must equal an even number of channel spacings'); end %= Backend % fname = fullfile( Q.SENSOR_DIR, Q.BACKEND_FILE ); % [H1,f_y,za_y,f_sensor] = hBackendFromFile( H1, fsensor1, za_sensor, ... f_obs, fname, Q.BACKEND_ORDER, Q.F_ORDER); H2 = hBackendFromFile( H2, fsensor2, za_sensor, ... f_obs, fname, Q.BACKEND_ORDER, Q.F_ORDER); %= Take difference % n = length( f_sensor ); % rows = repmat( (1:n)', 2, 1); cols = 1:(2*n)'; if ~Q.FREQSWITCH_REVERSE wgts = [ repmat(1,n,1); repmat(-1,n,1) ]; else wgts = [ repmat(-1,n,1); repmat(1,n,1) ]; end %= Create total H matrix H = sparse(rows,cols,wgts,n,n*2) * [H1;H2] * H; [f_y,za_y] = h_fix_ys( f_sensor, za_sensor ); clear H1 H2 fsensor1 fsensor2 rows cols wgts n end %=== saving always f_sensor and za_sensor % NOTE: needed by qp_Se, qp_Sxb, qp_assemble_K % save( [Q.OUT,'.f_sensor'], 'f_sensor'); save( [Q.OUT,'.za_sensor'], 'za_sensor'); %=== saving sensor without reduction % if do_save Hs = H; save( [Q.OUT,'.hs'], 'Hs' ); clear Hs end %=== Data reduction part ====================================================== if Q.BINNING_ON | Q.KRED_ON | Q.LRED_ON | Q.BINVIEW_ON % %= Sx is needed if KRED_ON or LRED_ON if Q.KRED_ON | Q.LRED_ON qp_Sx(Q); end [H,f_y,za_y,Hd] = qp_Hd( Q, H, Hd, f_y, za_y ); % end %=== Save if not return call and H is not empty % NOTE : H empty correspond to a H and Hd with reduction % specific for each retrieval point, where Hs and % Hds are saved as structures or individual Hs % amd Hds in the right mscript, e.g. hRedLimb % if do_save & ~isempty(H) save( [Q.OUT,'.h'], 'H' ); save( [Q.OUT,'.hd'], 'Hd' ); end % %= Save frequency and zenith angle grids % if do_save save( [Q.OUT,'.f_y'], 'f_y'); save( [Q.OUT,'.za_y'], 'za_y'); end