%------------------------------------------------------------------------ % NAME: qp_Sx % % Creates the Sx matrix, and derived quantities, corresponding % to the settings in Q. % % If the function is called without output argument (the first % version below), Sx, Sxinv and Dxinv are stored in files. % If there exist an output argument (that is, Dx), no saving % is done. % % Sxinv is the inverse of Sx. % Dx is the diagonal elements of Sx. % Dxinv is the inverse of Dx. % % FORMAT: qp_Sx( Q ) % or % Dx = qp_Sx( Q ) % % OUT: Dx Diagonal elements of Sx as a sparse matrix. % IN: Q Setting structure. % OPTIONAL: - % % TEMPLATE: - %------------------------------------------------------------------------ % HISTORY: 2001.03.25 Created by Patrick Eriksson. % 2001.08.26 Modified by Carlos Jimenez function Dx = qp_Sx( Q ) global REPORT_LEVEL r_l = Q.QP_LEVEL; %=== Create a temporary directory % tmpdir = temporary_directory( Q.TMP_AREA ); out(1,1); out(1,'Setting up Sx.'); %=== Create and append covariance matrices %= Species % p_sp = find(Q.RETRIEVAL_TAGS =='"'); p_gr = find(Q.SPECIES_KGRIDS =='"'); p_cov = find(Q.SPECIES_COVMATS =='"'); n_sp = size(p_sp,2)/2; for i=1:n_sp is = num2str(i); pa = num2str( 1 + p_sp(2*(i-1)+1)); pb = num2str(-1 + p_sp(2*(i-1)+2)); eval(['TAG = Q.RETRIEVAL_TAGS(',pa,':',pb,');']) pa = num2str( 1 + p_cov(2*(i-1)+1)); pb = num2str(-1 + p_cov(2*(i-1)+2)); eval(['COVMAT = Q.SPECIES_COVMATS(',pa,':',pb,');']) pa = num2str(p_gr(2*(i-1)+1)+1); pb = num2str(p_gr(2*(i-1)+2)-1); eval(['p_grid=Q.SPECIES_KGRIDS(',pa,':',pb,');']) f_grid=[Q.RETRIEVDEF_DIR,'/',p_grid]; eval(['r_grid = read_datafile(''',f_grid,''',','''VECTOR''',');']); % -- getting covariance matrix REPORT_LEVEL = 0; S = sFromFile([Q.RETRIEVDEF_DIR,'/',COVMAT],log10(r_grid)); REPORT_LEVEL = r_l; eval(['Sx_',is,'= S;']) % -- length and position of Sx if i==1 l_s = [1 length(r_grid)]; else l_s = [l_s l_s(length(l_s))+1 l_s(length(l_s))+length(r_grid)]; end end %= Pointing % if Q.POINTING_DO == 3 Sx_p = power(Q.POINTING_STDV,2); l_s= [l_s l_s(length(l_s))+1 l_s(length(l_s))+1]; end %= Temperature % if Q.TEMPERATURE_DO == 3 f_grid=[Q.RETRIEVDEF_DIR,'/',Q.TEMPERATURE_KGRID]; eval(['r_grid = read_datafile(''',f_grid,''',','''VECTOR''',');']); % -- getting covariance matrix REPORT_LEVEL = 0; Sx_t = sFromFile([Q.RETRIEVDEF_DIR,'/',Q.TEMPERATURE_COVMAT],log10(r_grid)); REPORT_LEVEL = r_l; l_s = [l_s l_s(length(l_s))+1 l_s(length(l_s))+length(r_grid)]; end %= Continuum absorption in primary band % if Q.CONTABS_DO == 3 disp('Continuuum absorption cannot be part of the state vector yet, sorry.') end %= Polynomial fit of baseline ripple % if Q.POLYFIT_DO == 3 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 ]; eval(['za_grid = read_datafile(''',za_grid,''',','''VECTOR''',');']); % for i = 1:n % file = sstring_get_i( Q.POLYFIT_COVMATS, i ); % % -- getting covariance matrix REPORT_LEVEL = 0; S = sFromFile([Q.RETRIEVDEF_DIR,'/',file],za_grid); eval(['Sx_pf',num2str(i),'=S;']) REPORT_LEVEL = r_l; % l_s = [l_s l_s(length(l_s))+1 l_s(length(l_s))+length(za_grid)]; % end end end %=== Append covariance matrices l_sl = l_s(length(l_s)); Sx = sparse(zeros(l_sl,l_sl)); c_l = 1; % -- species for i=1:n_sp ind = l_s(c_l):(l_s(c_l+1)); eval(['Sx(ind,ind) = Sx_',num2str(i),';']); c_l = c_l + 2; end % -- pointing if Q.POINTING_DO == 3 ind = l_s(c_l):(l_s(c_l+1)); Sx(ind,ind) = Sx_p; c_l = c_l + 2; end % -- temperature if Q.TEMPERATURE_DO == 3 ind = l_s(c_l):(l_s(c_l+1)); Sx(ind,ind) = Sx_t; c_l = c_l + 2; end % -- polynomial fit if Q.POLYFIT_DO == 3 n = Q.POLYFIT_ORDER + 1; % for i = 1:n ind = l_s(c_l):(l_s(c_l+1)); eval(['Sx(ind,ind) = Sx_pf',num2str(i),';']);; c_l = c_l + 2; end end %=== Calculate the other related matrices % np = size( Sx, 1 ); Dx = spdiags( spdiags(Sx,0), 0, np, np); % Dxinv = spdiags( 1./spdiags(Sx,0), 0, np, np); Sxinv = inv( Sx ); %=== Save % if nargout == 0 save( [Q.OUT,'.dxinv'], 'Dxinv' ); save( [Q.OUT,'.sxinv'], 'Sxinv' ); save( [Q.OUT,'.sx'], 'Sx' ); end %=== Delete the temporary directory % delete_tmp_dir( Q.TMP_AREA, tmpdir ); out(1,-1); return % -------------------------------------------------- %= Old Continuum absorption in primary band % left to be fixed if 0 % %= Calculate absorption for reference species QE = []; QE.ABS_TAGS = Q.CONTABS_REF_SPECIES; QE.PTZ = sstring_get_i( Q.SETUP_PTZ, 1 ); QE.VMRS = sstring_get_i( Q.SETUP_VMR_DIRS, 1 ); tmpdir2 = temporary_directory( Q.TMP_AREA ); template = which( 'calc_abs.tmplt' ); [cfile,basename] = qtool( Q, tmpdir2, template, QE ); qp_arts( Q, cfile ); Abs = read_artsvar( basename, 'abs' ); delete_tmp_dir( Q.TMP_AREA, tmpdir2 ); % %= Determine frequency of the fit points fp = qp_contabs_freqs( Q ); % %= 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; for i = 0:Q.CONTABS_ORDER s = sprintf('%s/sx%d.aa',tmpdir,i); eval(['sx_file',int2str(i),' = ''',s,''';']) 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( s, D2, 'AOMATRIX' ); end % gr_file = [ Q.RETRIEVDEF_DIR, '/', Q.CONTABS_KGRID ]; % fprintf(fid, 'VectorReadAscii( k_grid ) {\n "%s"\n}\n', gr_file ); fprintf(fid, 'VectorCalcLog10( k_grid, k_grid ) {\n}\n'); for i = 0:Q.CONTABS_ORDER fprintf(fid, 'sFromFile( k_grid ) {\n "%s"\n}\n', ... eval(['sx_file',int2str(i)]) ); fprintf(fid, 'sxAppend{\n}\n'); end end end % ------------------------------------------------------------------