%------------------------------------------------------------------------ % NAME: qp_assemble_K % % Loads a weighting function matrix calculated by ARTS and adds % weighting functions for quantities not handled by ARTS. % % FORMAT: [ K, k_names, k_index, k_aux ] = ... % qp_assemble_K( Q, basename, H, Hd, y, bl, % klevels [, x, iteration, QE] ) % % OUT: K Weighting function matrix. % k_names Name of each retrieval identity. % k_index Index of first and last element in x for each % retrieval identity. % k_aux As kx_aux in ARTS. % IN: Q Setting structure. % basename Basename for the performed ARTS run. % H The sensor/data-reduction transfer matrix. % Hd The data-reduction transfer matrix. % y Measurement spectra. % bl Fit of baseline ripple. % klevels Vector with the do-levels to consider. % See further qp_iter1.m. % OPTIONAL: x The state vector. If not given, it is assumed that % the a priori state is valid. Accordingly, this % argument shall be given for iterations > 1. % iteration If not given, it is assumed to be 1. % QE Extra Q structure. %------------------------------------------------------------------------ % HISTORY: 2001.08.31 Created by Patrick Eriksson. % 2003.11.08 Updated by Sho Tsujimaru. function [K,k_names,k_index,k_aux,Kb,kb_names] = qp_assemble_K(Q,basename,H,Hd,y,bl,klevels,x,iteration, QE) %=== Set default values for optional arguments % if ~exist( 'iteration' ) iteration = 1; end % if ~exist( 'QE' ) QE = []; QE.ABS_SAVE = 0; QE.ABS_EXIT = 0; end % % Note that if the order is changed here, qp_Sxb must be changed % accordingly. % %=== Load K variables (K is always stored as Kx) % if exist( [basename,'.kx.ab'], 'file' ) K = read_artsvar( basename, 'kx' ); k_aux = read_artsvar( basename, 'kx_aux' ); k_names = read_artsvar( basename, 'kx_names' ); k_lengths = read_artsvar( basename, 'kx_lengths' ); %=== Post-processing for some variables % %= Ground emission if in_klevels( Q.EGROUND_DO, klevels ) % f_mono = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.F_MONO), 'VECTOR' ); [ fp, findex, single_e ] = qp_eground_lims( f_mono, Q.EGROUND_LIMITS ); ne = length(fp); nf = length( f_mono ); % if ~single_e & length(fp)~=nf % iground = 1; i0 = 0; while ~strcmp( qp_kinfo(k_names{iground}), 'Ground emission' ) i0 = i0 + k_lengths(iground); iground = iground + 1; end % Kground = zeros( size(K,1), ne ); % for i = 1 : ne Kground(:,i) = sum( K(:,i0+(findex(i,1):findex(i,2)))' )'; end % K = [ K(:,1:i0), Kground, K(:,(i0+1+nf):size(K,2)) ]; % eap = interp1( f_mono, k_aux((i0+1):(i0+nf),2), fp ); k_aux = [ k_aux(1:i0,:); [fp,eap]; k_aux((i0+1+nf):size(K,2),:) ]; % k_lengths(iground) = ne; % k_names{iground} = 'Ground emission'; end end %---------------------------------------------------------------------------- %=== Apply H % K = H * K; else K = []; k_aux = []; k_names = []; k_lengths = []; end % === part corresponding to spectroscopic weighting function %=== Load Kb variables (Kb is always stored as Kb) if exist( [basename,'.k.ab'], 'file' ) Kb = read_artsvar( basename, 'k' ); kb_names = read_artsvar( basename, 'k_names' ); S_S = read_artsvar( basename, 'S_S' ); %=== Apply H % Kb = H * Kb; % save the variables %------------------- if isempty(Q.SAVE); Q.SAVE = Q.OUT;% additional field for figures to be saved; end strsave=[' save ', Q.SAVE, '.Kb Kb kb_names' ]; eval(strsave); strsave=[' save ', Q.SAVE, '.S_S S_S' ]; eval(strsave); else Kb = []; kb_names = []; S_S = []; end %============================================================================= %============================================================================= %=== Parts of K not covered by ARTS %============================================================================= %============================================================================= if in_klevels( Q.POLYFIT_DO, klevels ) load( [Q.OUT,'.f_sensor'], '-mat'); load( [Q.OUT,'.za_sensor'], '-mat'); nf = length( f_sensor ); nza = length( za_sensor ); np = Q.POLYFIT_ORDER + 1; %= Scale frequencies fsc = -1 + 2*(f_sensor-f_sensor(1)) / (last(f_sensor)-f_sensor(1)); Kpol = zeros( nf*nza, nza*np ); col = 0; findex = 1:nf; kpost_names = cell( np, 1 ); for i = 0:Q.POLYFIT_ORDER if i == 0 k = ones(nf,1); else k = fsc.^i; end for j = 1:nza col = col + 1; rindex = findex + (j-1)*nf; Kpol(rindex,col) = k; end kpost_names{i+1} = sprintf('Baseline polynomial: %d',i); end Kpost = Hd * Kpol; clear Kpol kpost_aux = [ repmat(vec2col(za_sensor),np,1), zeros(np*nza,1) ]; kpost_lengths = repmat( nza, np, 1 ); else Kpost = []; kpost_names = []; kpost_aux = []; kpost_lengths = []; end if in_klevels( Q.PPOLYFIT_DO, klevels ) load( [Q.OUT,'.f_sensor'], '-mat'); load( [Q.OUT,'.za_sensor'], '-mat'); nf = length( f_sensor ); nza = length( za_sensor ); np = Q.PPOLYFIT_ORDER + 1; npp = length(Q.PPOLYFIT_LIMITS) - 1; Kpol = zeros( nf*nza, nza*np*npp ); col = 0; ktmp_names = cell( np*npp, 1 ); % for l = 1:npp if l < npp pind = find( (f_sensor>=Q.PPOLYFIT_LIMITS(l)) & ... (f_sensor=Q.PPOLYFIT_LIMITS(l)) & ... (f_sensor<=Q.PPOLYFIT_LIMITS(l+1)) ); end if isempty(pind) s = sprintf(... 'The PPOLYFIT range %.4f-%.4f GHz does not cover any frequency',... Q.PPOLYFIT_LIMITS(l)/1e9,Q.PPOLYFIT_LIMITS(l+1)/1e9); error(s); end %= Scale frequencies fsc = -1 + 2*(f_sensor(pind)-f_sensor(pind(1))) / ... (last(f_sensor(pind))-f_sensor(pind(1))); for i = 0:Q.PPOLYFIT_ORDER if i == 0 k = ones(length(pind),1); else k = fsc.^i; end for j = 1:nza col = col + 1; rindex = pind + (j-1)*nf; Kpol(rindex,col) = k; end ktmp_names{ (l-1)*np + i + 1 } = ... sprintf('Baseline piecewise polynomial: %d/%d',i,l); end end Kpost = [ Kpost, Hd * Kpol ]; clear Kpol kpost_names = [ kpost_names; ktmp_names ]; kpost_aux = [ kpost_aux; repmat(vec2col(za_sensor),np*npp,1), ... zeros(np*npp*nza,1) ]; kpost_lengths = [ kpost_lengths; repmat( nza, np*npp, 1 ) ]; end if in_klevels( Q.TB_REFLOADS_DO, klevels ) if length( Q.TB_REFLOADS_NOMINAL ) ~= 2 error('TB_REFLOADS_NOMINAL must be a vector of length 2.'); end if any( Q.TB_REFLOADS_NOMINAL < 0 ) error('TB_REFLOADS_NOMINAL values should >= 0.'); end % Nominal brightness temperatures y1n = Hd * repmat(Q.TB_REFLOADS_NOMINAL(1),size(Hd,2),1); y2n = Hd * repmat(Q.TB_REFLOADS_NOMINAL(2),size(Hd,2),1); % Present estimated brightness temperatures if exist( 'x', 'var' ) i1 = size(K,2) + size(Kpost,2) + 1; y1 = Hd * repmat(x(i1),size(Hd,2),1); y2 = Hd * repmat(x(i1+1),size(Hd,2),1); else y1 = y1n; y2 = y2n; end % Use k1 as partial?esult k1 = (y2n-y1n) ./ ( (y2-y1).^2 ); k2 = (y1-(y-bl)) .* k1; k1 = (y-bl-y2) .* k1; Kpost = [ Kpost, k1, k2 ]; clear k1 k2 kpost_names = [ kpost_names; {'Reference load temperature(s)'} ]; kpost_aux = [ kpost_aux; 0, y1n; 0, y2n ]; kpost_lengths = [ kpost_lengths; 2 ]; end if in_klevels( Q.PROPCAL_DO, klevels ) if ~isscalar( Q.PROPCAL_TB_REF ) error('PROPCAL_TB_REF must be a scalar (length 1).'); end if Q.PROPCAL_TB_REF < 0 error('PROPCAL_TB_REF should be >= 0.'); end Kpost = [ Kpost, y - bl - Hd*repmat(Q.PROPCAL_TB_REF,size(Hd,2),1) ]; kpost_names = [ kpost_names; {'Proportional calibration error'} ]; kpost_aux = [ kpost_aux; 0, 0 ]; kpost_lengths = [ kpost_lengths; 1 ]; end % Newly appended part. %= Pressure Shift % % MHz/hPa is used as a unit. % if in_klevels( Q.PSF_DO, klevels ) % out(2, 'Calculating a part of K not covered by ARTS: Kpsf'); % dpsf = Q.PSF_DELTA; fname = fullfile( Q.SPECTRO_DIR, Q.PSF_FILE ); Lpsf = read_linefile(fname); n_Lpsf = length(Lpsf); if ~n_Lpsf error('No line found in PSF_FILE.'); else mname = Lpsf{1}.name; end % Spectroscopic databasse. if iteration < 2 fname0 = fullfile( Q.SPECTRO_DIR, Q.LINEFILE ); else fname0 = fullfile( Q.SPECTRO_DIR, 'backup_linefile' ); if ~exist( fname0, 'file' ) error('Backup line file cannot be found.'); end end L0 = read_linefile(fname0); % Database to be perturbed. fname = fullfile( Q.SPECTRO_DIR, 'dist_linefile' ); copyfile(fullfile( Q.SPECTRO_DIR, Q.LINEFILE ), fname); L = read_linefile(fname); for i = 1:n_Lpsf n = find_line_num( L0, Lpsf{i}.f ); psf = L{n}.psf; L{n}.psf = L{n}.psf + dpsf; end write_linefile(fname, L); % tmpdir = fileparts(basename); % tmpdir for perturbation calculation. tmpdir_dist = temporary_directory( Q.TMP_AREA ); eval(['!cp -r ', tmpdir, '/* ', tmpdir_dist]); Qdist = Q; Qdist.LINEFILE = 'dist_linefile'; % Creating template file. if iteration < 2 template = which( 'basic.tmplt' ); else template = which( 'cls_iter2_y.tmplt' ); end % Pertuebation calculation. [ cfile_dist, basename_dist ] = qtool( Qdist, tmpdir_dist, template, QE ); qp_arts( Qdist, cfile_dist ); % ydist = read_artsvar( basename_dist, 'y' ); y0 = read_artsvar( basename, 'y' ); %= Calculating jacobian for PSF. % Scaling factor due to unit change: Hz/Pa --> MHz/hPa. scfac = 1e-4; Kpsf = H * (ydist -y0) / (dpsf*scfac); clear ydist; clear Qdist; delete(fname); delete_tmp_dir( Q.TMP_AREA, tmpdir_dist ); Kpost = [ Kpost, Kpsf ]; kpost_names = [ kpost_names; {sprintf('Pressure shift: %s', mname)}]; kpost_aux = [ kpost_aux; 0, psf * scfac]; kpost_lengths = [ kpost_lengths; 1 ]; end %= Frequency Calibration % % MHz is used as internal unit. % if in_klevels( Q.FREQ_CALIB_DO, klevels ) % out(2, 'Calculating a part of K not covered by ARTS: Kfcal'); % df = Q.FREQ_CALIB_DELTA; % Scaling factor due to unit change: Hz --> MHz scfac = 1e-6; y0 = read_artsvar( basename, 'y' ); y0 = H * y0; %= Determine frequency of the fit points. f_mono = read_artsvar( basename, 'f_mono' ); flimits(1) = f_mono(1); flimits(2) = last( f_mono ); % if Q.FREQ_CALIB_ORDER == 0 fpoints = sum(flimits) / 2; else fpoints = linspace( flimits(1), flimits(2), Q.FREQ_CALIB_ORDER + 1 ); end % npoints = length(fpoints); Kfcal=[]; %= WF calculation % if Q.FREQ_CALIB_ORDER == 0 % f_dist = f_mono + df; Kfcal = calc_Kfcal( Q, basename, H, y0, f_dist, df, ... scfac, iteration, QE ); % elseif Q.FREQ_CALIB_ORDER > 0 % [m,n]=size(f_mono); % for ipoint=1:npoints f_dist = ones(m,n); % for ip=1:npoints if ip ~= ipoint f_dist ... = f_dist.*(f_mono-fpoints(ip))/(fpoints(ipoint)-fpoints(ip)); end end % f_dist = f_dist * df; f_dist = f_dist + f_mono; Ktmp = calc_Kfcal( Q, basename, H, y0, f_dist, df, ... scfac, iteration, QE ); Kfcal = [Kfcal, Ktmp]; % end % else error('Polynomial order should not be negative.'); end % Kfcal_aux = [ fpoints', fpoints' ]; Kfcal_aux(:,2) = 0; Kpost = [ Kpost, Kfcal ]; kpost_names = [ kpost_names; {sprintf('Freq. calibration')} ]; kpost_aux = [ kpost_aux; Kfcal_aux ]; kpost_lengths = [ kpost_lengths; Q.FREQ_CALIB_ORDER + 1 ]; end %= Sinusoidal baseline % if in_klevels( Q.SINEFIT_DO, klevels ) load( [Q.OUT,'.f_sensor'], '-mat'); load( [Q.OUT,'.za_sensor'], '-mat'); nf = length( f_sensor ); nza = length( za_sensor ); ns = length(Q.SINEFIT_PERIODS); %% This could be changed for reading of a file instead. p = Q.SINEFIT_PERIODS; Ksin = zeros( nf*nza, nza*ns*2 ); kpost_sin = cell( ns*2, 1); for it=1:ns vk = vec2col(2*pi*f_sensor/p(it)); k = [sin(vk) cos(vk)]; Ksin(:,2*it-1:2*it) = repmat(k, nza, 1); kpost_sin{2*it-1} = ... sprintf('Baseline sine fit: %g', Q.SINEFIT_PERIODS(it)); kpost_sin{2*it} = ... sprintf('Baseline cosine fit: %g', Q.SINEFIT_PERIODS(it)); end %% Kpost = [Kpost, Hd*Ksin]; Kpost = [Kpost, Ksin]; kpost_names = [kpost_names; kpost_sin]; za_double = reshape(repmat(vec2row(za_sensor),2,1),nza*2,1); kpost_aux = [kpost_aux; repmat(za_double,ns,1), zeros(ns*nza*2,1) ]; kpost_lengths = [kpost_lengths; repmat( nza, ns*2, 1)]; end %============================================================================= %=== Put parts together %============================================================================= %=== K % if ~isempty( Kpost ) % K = [ K, Kpost ]; % end %=== Help variables % if ~isempty( Kpost ) k_aux = [ k_aux; kpost_aux ]; k_lengths = [ k_lengths; kpost_lengths ]; k_names = [ k_names; kpost_names ]; end % nrq = size( k_names, 1 ); % k_index = [ [1;1+cumsum(k_lengths(1:(nrq-1)))], cumsum(k_lengths) ]; return %------------------------------------------------------------------------------ function b = in_klevels(level,klevels) b = any( level == klevels ); return %------------------------------------------------------------------------------ function n = find_line_num( L, freq ) nn = 0; for i = 1 : length(L) if L{i}.f == freq n = i; nn = nn + 1; end end if nn < 1 error('No line found!'); elseif nn > 1 error(sprintf('A line is identified with more than %d lines.', nn)); end return %------------------------------------------------------------------------------ function K = calc_Kfcal( Q, basename, H, y0, f_dist, df, scfac, iteration, QE ) % Creating temporary directory where perturbation calculation will be done. tmpdir = fileparts(basename); tmpdir_dist = temporary_directory( Q.TMP_AREA ); eval(['!cp -r ', tmpdir, '/* ', tmpdir_dist]); Qdist = Q; % if iteration < 2 template = which( 'basic.tmplt' ); fname = fullfile(Q.CALCGRIDS_DIR, 'f_mono_dist.aa'); write_datafile(fname, f_dist, 'vector'); Qdist.F_MONO = 'f_mono_dist.aa'; else template = which( 'cls_iter2_y.tmplt' ); end % [ cfile_dist, basename_dist ] = qtool( Qdist, tmpdir_dist, template, QE ); % In nonlinear iteration, f_mono is read not from CalcGrids but from tmpdir. if iteration >= 2 write_artsvar( basename_dist, 'f_mono', f_dist, 'b'); end qp_arts( Qdist, cfile_dist ) ; ydist = read_artsvar( basename_dist, 'y' ); ydist = H * ydist; % K = (ydist-y0) / (df * scfac); clear ydist; clear Qdist; delete_tmp_dir( Q.TMP_AREA, tmpdir_dist ); % return