%------------------------------------------------------------------------ % NAME: qp_errorplot % % Calculates, and optionally plots, the error due to the error % sources set to DO-levels 1 and 2. % % This function complements qp_invchar and works in the same way. % % FORMAT: [PTZ,E,Et] = qp_invchar( Q, Dy, Kx, kx_names, kx_index, kx_aux % [,do_plots, mr_limit] ) % % OUT: PTZ Matrix with pressure, temperature and altitude. % E Error for each source. % Et Name on the error sources as cell array. % IN: Q Setting structure. % Dy Contribution function matrix % Kx Weighting function matrix. % kx_names As the corresponding ARTS variable. % kx_index Start and stop index for each retrieval identity. % A 2 column matrix. % kx_aux As the corresponding ARTS variable. % OPTIONAL: do_plots Flag for doing plotting. Default is 1. % mr_limit Do not include altitudes with a measurement response % below this limit (only for plots). Default is 0.3. % To set MR_LIMIT=0, can result in crash. %------------------------------------------------------------------------ % HISTORY: 2002.01.21 Created by Patrick Eriksson function [PTZ,E,Et] = qp_errorplot(Q,Dy,Kx,kx_names,kx_index,kx_aux,do_plots,mr_limit) if ~exist('do_plots'), do_plots = 1; end if ~exist('mr_limit'), mr_limit = 0.3; end %=== This cell structure defines existing error sources % T = { 'MEASNOISE_DO', 'CALINOISE_DO', 'TEMPERATURE_DO', 'POINTING_DO', ... 'CONTABS_DO', 'EGROUND_DO', 'POLYFIT_DO', 'PPOLYFIT_DO' }; %=== Some sizes % np = size( Kx, 2 ); nrq = size( kx_names, 1 ); %=== Check that Kx and Dy have consistent sizes % if size(Dy,2)~=size(Kx,1) | size(Dy,1)~=np error('Sizes of Kx and Dy do not match. Recalculation needed?') end %=== Load H matrices and include Dy % % Dy is here included in H and Hd to save memory. % In this way, the observation error is obtained directly and the results % is stored in what will be the error covariance matrix. % out(1,'Loading H matrices ...'); % load( [Q.OUT,'.h'], '-mat' ); load( [Q.OUT,'.hd'], '-mat' ); %=== Calculate measurement response % out(1,'Calculating measurement response ...'); % A = Dy * Kx; % measres = zeros(np,1); % for i = 1:nrq ind = kx_index(i,1) : kx_index(i,2); measres(ind) = sum( A(ind,ind)' )'; end %=== Read ptz % PTZ = read_datafile( Q.APRIORI_PTZ, 'MATRIX', 1 ); %=== Init error variables % n = 0; E = []; Et = []; %=== Make copy of Q, and set DO=0 for all error sources % Q0 = Q; % for i = 1 : length(T) Q = setfield( Q, T{i}, 0 ); end %=== Treat thermal noise seperately % if Q0.MEASNOISE_DO==1 | Q0.MEASNOISE_DO==2 out(1,'Doing measurement noise ...'); % Q.MEASNOISE_DO = 2; % S = qp_Se( Q, H, Hd, 2 ); % n = n + 1; E(:,n) = sqrt( diag( Dy*S*Dy' ) ); Et{n} = 'Measurement noise';; % Q.MEASNOISE_DO = 0; end % if Q0.CALINOISE_DO==1 | Q0.CALINOISE_DO==2 out(1,'Doing calibration noise ...'); % Q.CALINOISE_DO = 2; % S = qp_Se( Q, H, Hd, 2 ); % n = n + 1; E(:,n) = sqrt( diag( Dy*S*Dy' ) ); Et{n} = 'Calibration noise'; % Q.CALINOISE_DO = 0; end %=== Remaining error sources (that are based on Kb) % Q = Q0; % Q.MEASNOISE_DO = 0; Q.CALINOISE_DO = 0; % out(1,'Calculating Sb ...'); Sb = qp_Sxb( Q, [1 2] ); % if size(Sb,1) > 0 % out(1,'Calculating Kb ...'); % tmpdir = temporary_directory( Q.TMP_AREA ); [fy, Kb, kb_names, kb_index] = qp_iter1( Q, H, Hd, tmpdir, [1 2], 0, 0 ); delete_tmp_dir( Q.TMP_AREA, tmpdir ); % for i = 1 : size(kb_index) % [group,name] = qp_kinfo( kb_names{i} ); out(1,sprintf('Doing %s ...',lower(name))); % ind = kb_index(i,1):kb_index(i,2); DK = Dy * Kb(:,ind); % n = n + 1; E(:,n) = sqrt( diag( DK * Sb(ind,ind) * DK' ) ); Et{n} = name; % end end %=== Return if not plots shall be made === %=== % %=== if ~do_plots %=== Return ? return %=== end %=== out(1,'Preparing plots ...'); %=== First figure - 1 if isempty(get(0,'Children')) % To handle the case when no window fig = 0; % is opened else fig = gcf - 1; end %=== Loop the retrieval identities ============================================ % ne = size(E,2); % for i = 1:nrq %= Names of present retrieval identity [group,name] = qp_kinfo( kx_names{i} ); %= Index in x of present retrieval identity ind = kx_index(i,1):kx_index(i,2); %= Get index for points with measurement response above limit ip = find( measres(ind) >= mr_limit ); %= Do something only if there are points > mr_limit if isempty(ip) out(1,sprintf('Retrieval identity: %s', name )); out(1,' all data below limit on measurement response'); else %= Create figure fig = fig + 1; figure(fig); fig_size(15,15,'cm'); %=== Scalar identity ====================================================== if length( ind ) == 1 row = 0.025; fsize = 10; y0 = 0.97; xoff = 0.02; textw = 0.4; set(gcf,'Unit','Cent'); set(gca,'Position',[-0.01 -0.01 1 1]); set(gca,'Box','Off'); set(gca,'TickLength',[0 0]) irow = 1; text2( xoff, y0-irow*row, sprintf('%s',name),fsize*1.4,1); irow = irow+2; for j = 1 : ne text2(xoff,y0-irow*row,Et{j},fsize,0); text2(xoff+textw,y0-irow*row,sprintf(': %.2e',E(ind,j)),fsize,0); irow = irow + 1; end %=== Vector identity ====================================================== else ip = ind(min(ip)):ind(max(ip)); %= Some plotting quantities % if strcmp( group, 'Species' ) % z(ind) = interpp( PTZ(:,1), PTZ(:,3), kx_aux(ind,1) ); % yfac = 1000; ytext = 'Altitude [km]'; yunit = 'km'; xfac = 100; xunit = '%'; % elseif strcmp( group, 'Temperature' ) % z(ind) = interpp( PTZ(:,1), PTZ(:,3), kx_aux(ind,1) ); % yfac = 1000; ytext = 'Altitude [km]'; yunit = 'km'; xfac = 1; xunit = 'K'; % elseif strcmp( group, 'Continuum' ) % z(ind) = interpp( PTZ(:,1), PTZ(:,3), kx_aux(ind,1) ); % yfac = 1000; ytext = 'Altitude [km]'; yunit = 'km'; xfac = 1e-3; xunit = '1/m'; % Internal unit for x is here 1/km % elseif strcmp( group, 'Ground emission' ) % ne = length( Q.EGROUND_LIMITS ); f_mono = read_datafile( [Q.CALCGRIDS_DIR,'/',Q.F_MONO], 'VECTOR' ); z(ind) = qp_eground_lims( f_mono, Q.EGROUND_LIMITS ); % yfac = 1e9; ytext = 'Frequency [GHz]'; yunit = 'GHz'; xfac = 1; xunit = '-'; % elseif strncmp( group, 'Baseline', 8 ) % global EARTH_RADIUS % z(ind) = za2geomtan( EARTH_RADIUS, Q.PLATFORM_ALTITUDE, kx_aux(ind,1)); % %= Revert IND and IP as tangent altitudes are coming in wrong order ind = flipdim(ind,2); ip = flipdim(ip,2); % yfac = 1e3; ytext = 'Geometrical ztan [km]'; yunit = 'km'; xfac = 1; xunit = 'K'; % else error(sprintf( 'Unknown retrieval identity: %s', name )); end sleg = []; h = []; for j = 1 : ne hv = pline( E(ip,j)*xfac, z(ip)/yfac, j ); hold on sleg = strvcat( sleg, Et{j} ); h = [h hv]; end axis([0 max(max(E(ip,:)))*xfac*1.2 0.99*min(z(ip))/yfac ... 1.01*max(z(ip))/yfac]) title2( name, 14 ); xlabel( ['Error [',xunit,']'] ); ylabel( ytext ); legend( h, sleg, 0 ); hold off end end end