%------------------------------------------------------------------------ % NAME: qp_Sxb % % Creates the covariance matrix (corresponding to x or b) for % arbitrary do_level. % % FORMAT: S = qp_Sx( 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). 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 ) % for i = 1:sstring_length(Q.RETRIEVAL_TAGS) % out(2,['Doing ',sstring_get_i(Q.RETRIEVAL_TAGS,i)]); % r_grid = [ Q.RETRIEVDEF_DIR, '/', sstring_get_i(Q.SPECIES_KGRIDS,i) ]; r_grid = read_datafile( r_grid, 'VECTOR' ); % s_file = [ 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( [Q.RETRIEVDEF_DIR,'/',Q.TEMPERATURE_KGRID],'VECTOR'); % Sloc = ... sFromFile( [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 %= Continuum absorption in primary band % if in_klevels( Q.CONTABS_DO, klevels ) % out(2,['Doing continuum absorption fit of order ',int2str(Q.CONTABS_ORDER)]); % Stmp = sparse(0,0); % % %= Calculate absorption for reference species 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 ); QE = []; QE.ABS_SAVE = 1; QE.ABS_EXIT = 1; tdlocal = temporary_directory( Q.TMP_AREA ); template = which( 'basic.tmplt' ); [cfile,basename] = qtool( Qtmp, tdlocal, template, QE ); qp_arts( Qtmp, cfile ); Abs = read_artsvar( basename, 'abs' ); clear Qtmp % %= Determine frequency of the fit points if length(Q.CONTABS_LIMITS) ~= 2 error('The length of Q.CONTABS_LIMITS must be 2.'); end if Q.CONTABS_ORDER == 0 fp = sum(Q.CONTABS_LIMITS) / 2; else fp = linspace( Q.CONTABS_LIMITS(1), Q.CONTABS_LIMITS(2), Q.CONTABS_ORDER+1 ); end % %= Read definition of the covariance matrix D = read_datafile( [Q.RETRIEVDEF_DIR,'/',Q.CONTABS_COVMAT], 'AOMATRIX'); % %= Loop frequency points and create covmat definition file in tmpdir p_abs = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.P_ABS], 'VECTOR'); f_mono = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.F_MONO], 'VECTOR'); np = length( p_abs ); Abs = [Abs(:,1),Abs,Abs(:,np)]; np = np + 2; % sx_file = sprintf('%s/sx_cont.aa',tdlocal); % r_grid = read_datafile( [Q.RETRIEVDEF_DIR,'/',Q.CONTABS_KGRID],'VECTOR'); r_grid = log10( r_grid ); % for i = 0:Q.CONTABS_ORDER D2 = []; D2{1} = D{1}; for j = 2:length(D) D2{j} = zeros( np, 3 ); D2{j}(:,1) = [ max(D{j}(:,1)); log10( p_abs ); min(D{j}(:,1)) ]; D2{j}(:,2) = interp1( D{j}(:,1), D{j}(:,2), ... D2{j}(:,1) ) .* interp1( f_mono, Abs, fp(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_tmp_dir( Q.TMP_AREA, tdlocal ); % 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( [Q.CALCGRIDS_DIR,'/',Q.F_MONO], 'VECTOR' ); [ fp, findex ] = qp_eground_lims( f_mono, Q.EGROUND_LIMITS ); % Sloc = sFromFile( [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 % za_grid = [ Q.SENSOR_DIR, '/', Q.ZA_SENSOR ]; za_grid = read_datafile( za_grid, 'VECTOR' ); % for i = 1:n % s_file = sstring_get_i( Q.POLYFIT_COVMATS, i ); % Sloc = sFromFile( [Q.RETRIEVDEF_DIR,'/',s_file], za_grid ); % nr = length( za_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 ]; 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 % za_grid = [ Q.SENSOR_DIR, '/', Q.ZA_SENSOR ]; za_grid = read_datafile( za_grid, 'VECTOR' ); % for j = 1 : (length(Q.PPOLYFIT_LIMITS)-1) % for i = 1:n % s_file = sstring_get_i( Q.PPOLYFIT_COVMATS, i ); % Sloc = sFromFile( [Q.RETRIEVDEF_DIR,'/',s_file], za_grid ); % nr = length( za_grid ); 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 %------------------------------------------------------------------------------ function b = in_klevels(level,klevels) b = any( level == klevels ); return