%------------------------------------------------------------------------ % NAME: qp_Sxc % % Creates the covariance matrix (corresponding to x or b) for % arbitrary do_level. Different to qp_Sx covaraince matrices % can be given as input for species and temperature. % % FORMAT: S = qp_Sxc( Q, Ss, St) % % OUT: S Created covariance matrix. % IN: Q Setting structure. % Ss Species covariance matrices % - Ss{i} structure where i % indexes the species as given in Q % - they have to be given in % the retrieval pressure grid % Q.SPECIES_KGRIDS % St Temperature covariance matrix % - given in grid % Q.TEMPERATURE_KGRID %------------------------------------------------------------------------ % HISTORY: 2003.05.01 Created by Carlos Jimenez (based on qp_Sx) function S = qp_Sxc( Q, Ss, St ) out(1,1); out(1,'Setting up Sx.'); Sx = qp_Sxaux( Q, Ss, St, 3 ); %=== Save % %= Sx is always stored. %= For CLS, the inverse of the diagonal elements of Sx is stored as Dxinv %= For OEM, the inverse of Sx is stored as Sxinv % if nargout == 0 % save( [Q.OUT,'.sx'], 'Sx' ); % if qpcls_methods( Q.RETRIEVAL_METHOD ) % np = size( Sx, 1 ); % Dxinv = spdiags( 1./spdiags(Sx,0), 0, np, np); save( [Q.OUT,'.dxinv'], 'Dxinv' ); % if strcmp( lower( Q.RETRIEVAL_METHOD ), 'oem' ) % Sxinv = inv( Sx ); save( [Q.OUT,'.sxinv'], 'Sxinv' ); % end end else S = Sx; end out(1,-1); return function S = qp_Sxaux( Q, Ss, St, 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' ); % Sloc = Ss{i}; 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 = St; % 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 %= Pointing % % 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( Q, 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 %------------------------------------------------------------------------------ function b = in_klevels(level,klevels) b = any( level == klevels ); return