%------------------------------------------------------------------------ % NAME: qp_Sxb % % Creates the covariance matrix (corresponding to x or b) for % arbitrary do_level. % % FORMAT: S = qp_Sxb( Q, klevels ) % % OUT: S Created covariance matrix. % IN: Q Setting structure. % klevels Vector with the do-levels to consider. % See further qp_iter1.m. %------------------------------------------------------------------------ % HISTORY: 2001.10.13 Created by Patrick Eriksson (based on old qp_Sx). % 2003.11.10 Updated by Sho Tsujimaru. function S = qp_Sxb( Q, klevels ) % % Note that the order here must be exactly the same as for qpi_K followed % by qp_assemble_K. % S = sparse(0,0); %= Species % if in_klevels( 3, klevels ) & ~isempty( Q.RETRIEVAL_TAGS ) % for i = 1:sstring_length(Q.RETRIEVAL_TAGS) % out(2,['Doing ',sstring_get_i(Q.RETRIEVAL_TAGS,i)]); % r_grid = fullfile( Q.RETRIEVDEF_DIR, sstring_get_i(Q.SPECIES_KGRIDS,i) ); r_grid = read_datafile( r_grid, 'VECTOR' ); % s_file = fullfile( Q.RETRIEVDEF_DIR, sstring_get_i(Q.SPECIES_COVMATS,i) ); % Sloc = sFromFile( s_file, log10(r_grid) ); % nr = length( r_grid ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Sloc ]; % end end %= Temperature % if in_klevels( Q.TEMPERATURE_DO, klevels ) % out(2,'Doing temperature profile'); % r_grid = read_datafile( fullfile( Q.RETRIEVDEF_DIR, Q.TEMPERATURE_KGRID ),... 'VECTOR'); % Sloc = ... sFromFile( fullfile( Q.RETRIEVDEF_DIR, Q.TEMPERATURE_COVMAT ), ... log10(r_grid) ); % nr = length( r_grid ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Sloc ]; % end %= Pointing % if in_klevels( Q.POINTING_DO, klevels ) % out(2,'Doing pointing off-set'); % nx = size( S, 1 ); % S(nx+1,nx+1) = power(Q.POINTING_STDV,2); % end %= Frequency % % MHz is used here as internal unit. % if in_klevels( Q.FREQUENCY_DO, klevels ) % out(2,'Doing frequency off-set'); % nx = size( S, 1 ); % S(nx+1,nx+1) = power(Q.FREQUENCY_STDV/1e6,2); % end %= Continuum absorption in primary band % if in_klevels( Q.CONTABS_DO, klevels ) % out(2,['Doing continuum absorption fit of order ',int2str(Q.CONTABS_ORDER)]); % %= Determine frequency of the fit points [ flimits, fp ] = qp_contabs_freqs(Q); % %= Calculate absorption at FP for reference species. NOTE, absorption in 1/km Qtmp = Q; Qtmp.RETRIEVAL_TAGS = Q.CONTABS_REF_SPECIES; Qtmp.OTHER_TAGS = []; Qtmp.APRIORI_VMR = sstring_get_i( Q.SETUP_VMR, 1 ); Qtmp.APRIORI_PTZ = sstring_get_i( Q.SETUP_PTZ, 1 ); % Abs = 1000 * qp_abs(Qtmp,fp); % clear Qtmp % %= Read definition of the covariance matrix D = read_datafile( fullfile( Q.RETRIEVDEF_DIR, Q.CONTABS_COVMAT ), ... 'AOMATRIX' ); % %= Loop frequency points and create covmat definition file in tmpdir p_abs = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.P_ABS ), 'VECTOR'); np = length( p_abs ); Abs = [Abs(:,1),Abs,Abs(:,np)]; np = np + 2; % sx_file = [ tempname, '_sx.aa' ]; % r_grid = read_datafile( fullfile( Q.RETRIEVDEF_DIR, Q.CONTABS_KGRID ), ... 'VECTOR'); r_grid = log10( r_grid ); % Stmp = sparse(0,0); % for i = 0:Q.CONTABS_ORDER D2 = []; D2{1} = D{1}; for j = 2:length(D) D2{j} = zeros( np, 3 ); D2{j}(:,1) = join_grids( {max(D{j}(:,1)),log10(p_abs),min(D{j}(:,1))},... 1e-6, 0, 1 ); D2{j}(:,2) = interp1( D{j}(:,1), D{j}(:,2), D2{j}(:,1) ) .* Abs(i+1,:)'; D2{j}(:,3) = interp1( D{j}(:,1), D{j}(:,3), D2{j}(:,1) ); end % write_datafile( sx_file, D2, 'AOMATRIX' ); % Sloc = sFromFile( sx_file, r_grid ); % nr = length( r_grid ); nx = size( Stmp, 1 ); % Stmp = [ Stmp, sparse(nx,nr); sparse(nr,nx), Sloc ]; % end % nr = size( Stmp, 1 ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Stmp ]; % delete( sx_file ); % end %= Ground emission % if in_klevels( Q.EGROUND_DO, klevels ) % out(2,sprintf('Doing ground emission (seperated in %d range(s))',... length(Q.EGROUND_LIMITS)-1 )); % f_mono = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.F_MONO), 'VECTOR' ); fp = qp_eground_lims( f_mono, Q.EGROUND_LIMITS ); % Sloc = sFromFile( fullfile( Q.RETRIEVDEF_DIR, Q.EGROUND_COVMAT ), fp ); % nr = length( fp ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Sloc ]; % end %= Polynomial fit of baseline ripple % if in_klevels( Q.POLYFIT_DO, klevels ) % out(2,['Doing polynomial baseline fit of order ',int2str(Q.POLYFIT_ORDER)]); % Stmp = sparse(0,0); % n = Q.POLYFIT_ORDER + 1; % if n ~= sstring_length( Q.POLYFIT_COVMATS ) error(... 'The length of POLYFIT_COVMATS and the polynomial order do not agree'); end % load( [Q.OUT,'.za_sensor'], '-mat'); % for i = 1:n % s_file = sstring_get_i( Q.POLYFIT_COVMATS, i ); % Sloc = sFromFile( fullfile( Q.RETRIEVDEF_DIR, s_file), za_sensor ); % nr = length( za_sensor ); nx = size( Stmp, 1 ); % Stmp = [ Stmp, sparse(nx,nr); sparse(nr,nx), Sloc ]; end % nr = size( Stmp, 1 ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Stmp ]; end %= Piecewise polynomial fit of baseline ripple % if in_klevels( Q.PPOLYFIT_DO, klevels ) % out(2,['Doing piecewise polynomial baseline fit of order ', ... int2str(Q.PPOLYFIT_ORDER)]); % Stmp = sparse(0,0); % n = Q.PPOLYFIT_ORDER + 1; % if n ~= sstring_length( Q.PPOLYFIT_COVMATS ) error(... 'The length of PPOLYFIT_COVMATS and the polynomial order do not agree'); end % load( [Q.OUT,'.za_sensor'], '-mat'); % for j = 1 : (length(Q.PPOLYFIT_LIMITS)-1) % for i = 1:n % s_file = sstring_get_i( Q.PPOLYFIT_COVMATS, i ); % Sloc = sFromFile( fullfile( Q.RETRIEVDEF_DIR, s_file), za_sensor ); % nr = length( za_sensor ); nx = size( Stmp, 1 ); % Stmp = [ Stmp, sparse(nx,nr); sparse(nr,nx), Sloc ]; % end end % nr = size( Stmp, 1 ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Stmp ]; end %= Temperature of reference loads % if in_klevels( Q.TB_REFLOADS_DO, klevels ) % out(2,'Doing reference load temperatures'); % if length( Q.TB_REFLOADS_STDV ) ~= 2 error('TB_REFLOADS_STDV must be a vector of length 2.'); end % nx = size( S, 1 ); % Stmp = diag( power( Q.TB_REFLOADS_STDV,2 ) ); % S = [ S, sparse(nx,2); sparse(2,nx), Stmp ]; end %= Proportional calibration error % if in_klevels( Q.PROPCAL_DO, klevels ) % out(2,'Doing proportional calibration error'); % if Q.TB_REFLOADS_DO error('You should not use PROPCAL and TB_REFLOADS simultaneously.'); end % if ~isscalar( Q.PROPCAL_STDV ) error('PROPCAL_STDV must be a scalar (length 1).'); end % nx = size( S, 1 ); % S(nx+1,nx+1) = power(Q.PROPCAL_STDV,2); end % Newly appended part. %= Pressure Shift % % MHz/HPa is used as internal unit. % if in_klevels( Q.PSF_DO, klevels ) % out(2,'Doing pressure-shift'); % nx = size( S, 1 ); % S(nx+1,nx+1) = power(Q.PSF_STDV/1e4,2); % end %= Frequency Calibration % % MHz is used as internal unit. % if in_klevels( Q.FREQ_CALIB_DO, klevels ) % if Q.FREQUENCY_DO error('Both FREQUENCY_DO and FREQ_CALIB_DO are turned on.'); end % if Q.FREQ_CALIB_ORDER < 0 error('Polynomial order should not be negative.'); end % out(2,['Doing freq. calibration fit of order ', ... int2str(Q.FREQ_CALIB_ORDER)]); %= Determine frequency of the fit points. f_mono = read_datafile( fullfile(Q.CALCGRIDS_DIR,Q.F_MONO), 'vector' ); flimits(1) = f_mono(1); flimits(2) = last( f_mono ); if Q.FREQ_CALIB_ORDER == 0 fp = sum(flimits) / 2; else fp = linspace( flimits(1), flimits(2), Q.FREQ_CALIB_ORDER + 1 ); end % s_file = fullfile( Q.RETRIEVDEF_DIR, Q.FREQ_CALIB_COVMAT ); % Sloc = sFromFile( s_file, fp ); % % Rescaling due to the change of unit: Hz ---> MHz. % This is essential for avoiding numerical instability. Sloc = Sloc / 1e12; % nr = length( fp ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Sloc]; % end %= Sinusoidal fit of baseline % if in_klevels( Q.SINEFIT_DO, klevels ) n = length(Q.SINEFIT_PERIODS); out(2,['Doing sinusoidal baseline fit with ' int2str(n) ' periods.']); if n ~= sstring_length( Q.SINEFIT_COVMATS ) error(... 'The length of SINEFIT_COVMATS and the number of periods do not agree'); end Stmp = sparse(0,0); load( [Q.OUT,'.za_sensor'], '-mat'); for it = 1:n s_file = sstring_get_i( Q.SINEFIT_COVMATS, it ); Sloc = sFromFile( fullfile( Q.RETRIEVDEF_DIR, s_file), za_sensor ); nStmp = size(Stmp,1); nSloc = size(Sloc,1); %% nr = length( za_sensor ); %% nx = size( Stmp, 1 ); Stmp = [ Stmp, sparse(nStmp,2*nSloc); ... sparse(nSloc,nStmp), Sloc, sparse(nSloc,nSloc); ... sparse(nSloc,nStmp+nSloc), Sloc ]; end % nr = size( Stmp, 1 ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Stmp ]; end if in_klevels( Q.PSINEFIT_DO, klevels ) % out(2,['Doing piecewise sinusodial baseline fit with ', ... int2str(length(Q.PSINEFIT_PERIODS)) ' periods']); % Stmp = sparse(0,0); n = length(Q.PSINEFIT_PERIODS); % if n > sstring_length( Q.PSINEFIT_COVMATS ) error(... 'The length of PSINEFIT_COVMATS and the number of periods does not agree'); end % load( [Q.OUT,'.za_sensor'], '-mat'); % for j = 1 : (length(Q.PSINEFIT_LIMITS)-1) % for i = 1:n % s_file = sstring_get_i( Q.PSINEFIT_COVMATS, i ); % Sloc = sFromFile( fullfile( Q.RETRIEVDEF_DIR, s_file), za_sensor ); % nStmp = size(Stmp,1); nSloc = size(Sloc,1); %% nr = length( za_sensor ); %% nx = size( Stmp, 1 ); Stmp = [ Stmp, sparse(nStmp,2*nSloc); ... sparse(nSloc,nStmp), Sloc, sparse(nSloc,nSloc); ... sparse(nSloc,nStmp+nSloc), Sloc ]; end end % nr = size( Stmp, 1 ); nx = size( S, 1 ); % S = [ S, sparse(nx,nr); sparse(nr,nx), Stmp ]; end %------------------------------------------------------------------------------ function b = in_klevels(level,klevels) b = any( level == klevels ); return