%------------------------------------------------------------------------ % NAME: qp_Hd % % Sets up Hd and includes this into H. % % See the README file for options. % % FORMAT: [H,f_y,za_y,Hd] = qp_Hd(Q,H,Hd,f_y,za_y) % % OUT: H % f_y % za_y % Hd % IN: Q The setting structure. % H % Hd % f_y % za_y %------------------------------------------------------------------------ % HISTORY: 2001.03.29 Created by Patrick Eriksson. % 2001.07.30 Adding red with Kx by Carlos Jimenez % 2001.09.02 Adding red for limb by Carlos Jimenez function [H,f_y,za_y,Hd] = qp_Hd(Q,H,Hd,f_y,za_y) %=== Average of subsequent spectra if Q.BINVIEW_ON % [H,f_y,za_y,Hd] = hBinView( H, f_y, za_y, Hd, Q.BINVIEW_DATA ); % end %=== Averaging of neighbouring channels (binning) if Q.BINNING_ON % [H,f_y,za_y,Hd] = hBinFile( Hd, f_y, za_y, Q.BINNING_FILE ); % end if Q.KRED_ON %=== Data reduction by Kx % [H,f_y,za_y,Hd] = hRedKx(Q,H,f_y,za_y,Hd); % elseif Q.LRED_ON %=== Data reduction by Limb % if Q.KRED_ON error('LRED and KRED both on?'); end [H,f_y,za_y,Hd] = hRedLimb(Q,H,f_y,za_y,Hd); % end return %------------------------------------------------------------------------ % % SUBFUNCTIONS % %------------------------------------------------------------------------ %------------------------------------------------------------------------ % NAME: hRedKx % % Includes reduction based on Kx into Hd. The basic idea is to % do a Hotelling transformation of the spectral space. But as % the spectral variability can be decomposed as: % % Sy = KxSxKx' + KbSbKb' + Se % % and we are interested in capture the variability of the % parameters to be retrieved (x), the transformation is % calculated only from that variability: % % KxSxKx' = E V E' % % where E are the eigenvectors and V the eigenvalues. % % The Hotteling transformation is derived as: % % 1. Notice that % KxSxKx' = Kx sqrt(Sx) sqrt(Sx)' Kx' = Kx sqrt(Sx) (Kx sqrt(Sx))' % where sqrt(Sx) is one of the possible square roots of Sx. % % 2. Then a SVD % Kx sqrt(Sx) = U M V' % where U and V are left and right singular vectors. % % 3. The eigenvectors are then U as from (1): % E V E' = U M V' (U M V')' = U M M U' % % Different assumptions are implemented for sqrt(Sx): % % 1.Q.KRED_DEPTH = 0 % sqrt(Sx) = I % i.e. Sx is assumed to be the identity matrix % % 2.Q.KRED_DEPTH = 1 % sqrt(Sx) = sqrt(diag(Sx)) % i.e. Sx is assumed to be diagonal % % 1.Q.KRED_DEPTH = 2 % sqrt(Sx) = chol(Sx) % i.e. full Sx is used and the sqrt is done by a choleski % decomposition. % % More info in: % Eriksson, P., Jimenez, C., Búhler, S. and Murtagh, D. "A % Hotteling Transformation approach for rapid inversion of % atmospheric spectra", J. Quant Rad. Spectroscopy, pp % 2001. % % % % FORMAT: [H,f_y,za_y,Hd] = hRedKx(Q,Hd,f_y,za_y,Hd) % % RETURN: Q Q structure % H total H matrix after binning % f_y new frequency vector % za_y new zenith angle vector % Hd data reduction H matrix after binning % IN: H total H matrix before reduction % f_y input frequency vector % za_y input zenith angle vector % Hd data reduction H matrix before reduction % % % HISTORY: 01.07.18 Created by Carlos Jimenez/ Patrick Eriksson %------------------------------------------------------------------------ function [H,f_y,za_y,Hd] = hRedKx(Q,H,f_y,za_y,Hd) out(1,1); out(1,'Setting up reduction matrix based on Kx (H)'); %=== Get Kx tmpdir = temporary_directory( Q.TMP_AREA ); [fy,Kx] = qp_iter1(Q,H,Hd,tmpdir,3); %=== Scale KX by elements of SX if Q.KRED_DEPTH == 2 load( [Q.OUT,'.sx'], '-mat' ); np = size( Sx, 1 ); Kx = h_x_h(Kx,chol(Sx)); elseif Q.KRED_DEPTH == 1 load( [Q.OUT,'.dxinv'], '-mat' ); np = size( Dxinv, 1 ); Kx = h_x_h(Kx,sqrt(inv(Dxinv))); end %=== Calculate the reduction transfer matrix [U,S,V] = svd( Kx, 0 ); clear S V Kx if Q.KRED_N > size(U,2) error(sprintf( ... '%d number of eigenvectors cannot be provided.\nOnly %d are possible.',... Q.KRED_N,size(U,2))); end U = U( :, 1:Q.KRED_N )'; H = h_x_h( U, H ); Hd = h_x_h( U, Hd ); clear U f_y = 1:Q.KRED_N; za_y = ones( Q.KRED_N, 1 ); %=== Delete the temporary directory delete_tmp_dir( Q.TMP_AREA, tmpdir ); out(1,-1); return %------------------------------------------------------------------------ % NAME: hRedLimb % % Includes reduction based on Limb into Hd. By Limb it means % a special reduction based on combining the reduction based % on Kx with a prior optimization of the spectral input % from a limb observation regarding each specific point of a % retrieval grid. This reduction is useful for a retrieval % method where each element of the retrieved species state % vector is obtained independently, so the spectral input % in the limb observation can be optimized for each element. % So far we use it with the net retrievals from the package % Npack. % % The basic idea is to decide which elements for the spectral % input are left for each retrieval point, and then to proceed % with a normal Hotelling transformation of the new spectral % space. In practice is done by doing the Hot transf in a Kx % padded with zeros to leave only the rows (part of the scan) % and columns (the retrieval tag) corresponding to each % retrieval tag and point. % % Notice that now for each point of the retrieval grid there % will be a reduction matrix, also saved in Hd as described % below. The same applies to the subsequent H. % % The parameters for this reduction are given as usual as % part of Q, see README. % % For details of the Hotelling transformation see hRedKx. % % % FORMAT: [H,f_y,za_y,Hd] = hRedLimb(Q,Hd,f_y,za_y,Hd) % % RETURN: Q Q structure % H [] empty to signal in qp_H that H should not be % saved again as it is saved here. [Q.OUT,'.h'] % will contain Hij matrices corresponding to % the tag i and ret point j. % f_y new frequency vector % za_y new zenith angle vector % Hd [] empty as in H, same saving strategy % % IN: H total H matrix before reduction % f_y input frequency vector % za_y input zenith angle vector % Hd data reduction H matrix before reduction % % TEMPLATE: iter1.tmplt % % HISTORY: 01.08.14 Created by Carlos Jimenez %------------------------------------------------------------------------ function [H,f_y,za_y,Hd] = hRedLimb(Q,H,f_y,za_y,Hd) out(1,1); out(1,'Setting up reduction matrix optimizing the limb scan (H)'); out(2,'Calculating Kx'); tmpdir = temporary_directory( Q.TMP_AREA ); % -- Kx without sensor [fy,kx,kx_names,kx_index,kx_aux] = qp_iter1(Q,1,1,tmpdir,3); % f_mono = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.F_MONO ), 'VECTOR' ); p_abs = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.P_ABS ), 'VECTOR'); za_pencil = read_datafile( fullfile( Q.CALCGRIDS_DIR, Q.ZA_PENCIL ), ... 'VECTOR'); % -- z_tan global EARTH_RADIUS z_tan = za2geomtan(EARTH_RADIUS,Q.PLATFORM_ALTITUDE,za_pencil); % -- z_abs ptz = read_datafile(Q.APRIORI_PTZ,'MATRIX'); z_abs = interpp(ptz(:,1),ptz(:,3),p_abs); t_abs = interpp(ptz(:,1),ptz(:,2),p_abs); if Q.HSE_IN_ON % -- obtaining a priori and interpolating to p_abs PSP = read_datafile([Q.APRIORI_VMR,'.H2O.aa'],'MATRIX',1); h2o_abs = interpp(PSP(:,1),PSP(:,2),p_abs); z_abs = hseCalc(p_abs,t_abs,z_abs,h2o_abs,Q.HSE_PREF,Q.HSE_ZREF,EARTH_RADIUS,1); end % nrq = sstring_length(Q.RETRIEVAL_TAGS); nf = length(f_mono); nza = length(za_pencil); clear za_pencil %=== Setting new H and Hd % % To accomodate the new H and Hd - one for each % ret point and ret grid, matrices will be saved % on [Q.OUT,'.h'] and [Q.OUT,'.hd'] as usual % but with one matrix for tag and ret point % following % Hij % Hdij i=tag number j=ret point %=== Do data reduction from each point of each ret grid out(2,'Doing data reduction'); h = H; hd = Hd; clear H Hd % -- loading retrieval grid p_gr = find(Q.LRED_KGRIDS =='"'); n_sp = size(p_gr,2)/2; for i=1:n_sp is = num2str(i); pa = num2str(p_gr(2*(i-1)+1)+1); pb = num2str(p_gr(2*(i-1)+2)-1); eval(['p_grid=Q.LRED_KGRIDS(',pa,':',pb,');']) f_grid=[Q.RETRIEVDEF_DIR,'/',p_grid]; eval(['p_grid',is,'= read_datafile(''',f_grid,''',','''VECTOR''',');']); end % -- are the grid equal and whole species? % then reduction common for all species % done only once if ~Q.LRED_INDTAGS do_one = 1; for i=1:n_sp-1 if ~strcmp(sstring_get_i(Q.LRED_KGRIDS,i),sstring_get_i(Q.LRED_KGRIDS,i+1)) do_one = 0; end end if do_one nrq = 1; end end for i = 1:nrq is = num2str(i); out(2,['Tag ',is]); eval(['p_grid = p_grid',is,';']); z_grid = interpp(p_abs,z_abs,p_grid); l_grid = length(z_grid); % -- leaving only eigenvectord for i species if Q.LRED_INDTAGS ind = kx_index(i,1):kx_index(i,2); kxi = kx(:,ind); else kxi = kx; end for j = 1:l_grid js = num2str(j) out(2,[' Retrival altitude: ',num2str(z_grid(j)/1e3),' km']); [h_n,hd_n] = limbKx(Q,h,hd,kxi,z_grid(j),z_tan,nf,nza); eval(['H',is,js,'=h_n;']) if i==1 & j==1 eval(['save ',Q.OUT,'.h H',is,js]) else eval(['save ',Q.OUT,'.h H',is,js,' -append']) end eval(['clear H',is,js,' h_n;']) eval(['Hd',is,js,'=hd_n;']) if i==1 & j==1 eval(['save ',Q.OUT,'.hd Hd',is,js]) else eval(['save ',Q.OUT,'.hd Hd',is,js,' -append']) end eval(['clear Hd',is,js,' hd_n;']) end end %=== Saving aux var nrq and l_grid ind = []; for i=1:nrq eval(['l = length(p_grid',is,');']); ind = [ind l]; end eval(['save ',Q.OUT,'.hd ind -append']) eval(['save ',Q.OUT,'.h ind -append']) %=== To signal to qp_H that H and Hd should not been % saved as they are saved here. H = []; Hd = []; %=== New f_y and za_y f_y = 1:Q.LRED_N; za_y = ones( Q.LRED_N, 1 ); %=== Delete the temporary directory delete_tmp_dir( Q.TMP_AREA, tmpdir ); out(1,-1); % %-------------------------------------------- ---------------------------- function [H,Hd] = limbKx(Q,H,Hd,Kx,z,ztan,nf,nza); %------------------------------------------------------------------------ %=== Looking for z + HISTEP > z_tan > z - LOSTEP hi = Q.LRED_HISTEP; lo = Q.LRED_LOSTEP; z_ind = find( ( (z+hi)>ztan ) & ( ztan > (z-lo) )); % -- special cases if isempty(z_ind) % %-- when z is lower or higher than any z_tan then % kept as many scans as REDL_HISTEP or LOSTEP % from the lowest or higher z_tan if z < ztan(length(ztan)) % z_ind = find( (ztan(length(ztan))+hi)>ztan ); % elseif z > ztan(1) % z_ind = find( ztan > (ztan(1)-lo)); % end % end %=== Building the Kx corresponding to the optimized grid % -- 1st leaving only the Kx corresponding to the nrq species % (already done here) % -- 2nd leaving only the Kx right channels % (padding with zeros the rows which channels are excluded) nind = length(z_ind); ind = (1+ nf * (z_ind(1)-1)):(nf * z_ind(length(z_ind))); Ka = Kx(ind,:); Kx = zeros(size(Kx)); Kx(ind,:)= Ka; %=== Apply sensor to Kx Kx = h_x_h(H,Kx); %=== Scale Kx by elements of Sx if Q.LRED_DEPTH == 2 load( [Q.OUT,'.sx'], '-mat' ); np = size( Sx, 1 ); Kx = h_x_h(Kx,chol(Sx)); elseif Q.LRED_DEPTH == 1 load( [Q.OUT,'.dxinv'], '-mat' ); np = size( Dxinv, 1 ); Kx = h_x_h(Kx,sqrt(inv(Dxinv))); end %=== Calculate the reduction transfer matrix [U,S,V] = svd( Kx, 0 ); clear S V Kx if Q.LRED_N > size(U,2) error(sprintf( ... '%d number of eigenvectors cannot be provided.\nOnly %d are possible.',... Q.LRED_N,size(U,2))); end %=== saving results U = U( :, 1:Q.LRED_N )'; H = sparse(H); H = h_x_h( U, H ); Hd = sparse(Hd); Hd = h_x_h( U, Hd );