%-------------------------------------------------------------- % NAME: qp_plotspectro % % Calculates and plots the spectroscopic error. % Input: Q structure. Here are defined also the parameters which % have to be ploted% % kx_names As the corresponding ARTS variable. % kx_index Start and stop index for each retrieval identity. % A 2 column m% % kx_aux As the corresponding ARTS variable. % kx_aux_abs The same as kx_aux but in absolute values % measresp The precalculated measurement response % s_tot The precalculated total statistical error % Optional: mr_limit % band % % Output: a plot function qp_plotspectro(Q,PTZ, Kx, kx_names, kx_index, kx_aux, measres, Dy, s_tot, mr_limit); % this is just for the title of the plot if ~exist('band'), band = ' ' ; else band=['Band', band]; end if ~exist('mr_limit'), mr_limit = 0.3; end % loading the matrix holding the spectroscopic uncertainties % and weighting function if exist ( [Q.SAVE,'.S_S'], 'file') out(3, '==== Loading the matrix holding the spectroscopic error ==='); load( [Q.SAVE,'.S_S'], '-mat' ); else error(['Matrix ' , Q.SAVE,'.S_S', ' does not exist' ]); end % if exist ( [Q.SAVE,'.Kb'], 'file') out(3, '==== Loading Kb matrix ==='); load( [Q.SAVE,'.Kb'], '-mat' ); else error(['Matrix ' , Q.SAVE,'.Kb', ' does not exist' ]); end % some sizes %----------- % number of spectro parameters nps = size( Kb, 2 ); nss = size( S_S, 1 ); %=== Some sizes % np = size( Kx, 2 ); nrq = size( kx_names, 1 ); % % size of measurement nm1 = size(Dy,2); nm2 = size(Kb,1) ; % check some sizes if size(Dy,2)~=size(Kb,1) error('Sizes of Kb and Dy do not match. Recalculation needed?') end if nps~= nss error('Size se of Kb and S_S do not match. Recalculation needed?') end % separate Kb, kb_names, ... for each parameter out(3, '==== Separate the matrix for each parameter ==='); nr_param=0; % initialize the parameter if Q.INTENS_ON nr_param=nr_param +1; name_param{nr_param}='intensity'; end if Q.POSITION_ON nr_param=nr_param +1; name_param{nr_param}='position'; end if Q.AGAM_ON nr_param=nr_param +1; name_param{nr_param}='agam'; end if Q.SGAM_ON nr_param=nr_param +1; name_param{nr_param}='sgam'; end if Q.NAIR_ON nr_param=nr_param +1; name_param{nr_param}='nair'; end if Q.NSELF_ON nr_param=nr_param +1; name_param{nr_param}='nself'; end if Q.PSHIFT_ON nr_param=nr_param +1; name_param{nr_param}='Psift'; end % number of lines for one specific parameter npss =nps/nr_param; for j=1:nr_param Kb_param = Kb(:, npss*(j-1)+1:npss*j); S_S_param = S_S(npss*(j-1)+1:npss*j, 1); kb_names_param = kb_names(npss*(j-1)+1:npss*j); % get the tag for the species for which the weighting function % is specified to be calculated spec = strread(Q.SPECTRO_TAGS, '%s', 'delimiter', ','); nr_ss= length(spec); for i =1 : nr_ss spec{i} = spec{i}(2:length(spec{i})-1); %just to extract the %quotation marks s = strread(spec{i}, '%s', 'delimiter', '-'); if length(s) > 1 name_spec{i} = [char(s(1)),'-', char(s(2))] ; %take only species and isotope else name_spec{i} =[char(s), '-*']; end end % get the name of the species found in kb_names for i=1:npss s=strread(kb_names_param{i}, '%s', 'delimiter', '/'); dpar = s(2); s2 = strread(s{1}, '%s', 'delimiter', '@'); line_freq = s2(2); s3 = strread(s2{1}, '%s', 'delimiter', '-'); name_spec2{i}=[char(s3(1)),'-', char(s3(2))]; % rename kb_names; this is just to make it nicer kb_names_param{i} =deblank( [char(name_spec2{i}), '@', char(line_freq), '|', char(dpar)]); end % initialize the correlation matrix; this creates a matrix for the uncorrelated case, % e.g., identity matrix. Mcorr =diag(ones(1, npss) , 0); % set the default value for the correlation (0) % this is to handle the case when no correlation is specified % in control file if ~exist('Q.CORRELATION'); Q.CORRELATION =zeros (nr_ss); end % rearange with respect to the species; Kb_new = []; S_S_new = []; kb_names_new = {}; nr_line = 0; % ni =0; % keeps the number of species for which weithing % functions are calculated. If a species is specified in % wtgss but no lines corresponding to the species are % found in the spectral data, then this species is not counted for i =1: nr_ss I =[]; % get the index for the lines beloging to each species % specified in wtgss. for ij =1:npss if strcmp(char(name_spec{i}), char(name_spec2{ij})) I =[I, ij]; end end % do something only lines correspomding to a specific species % are found % rearange with respect to the species; if length(I) ni =ni +1; Kb_spec = Kb_param(:, I); S_S_spec = S_S_param(I,1); kb_names_spec = kb_names_param(I); % Kb_new = [Kb_new,Kb_spec]; S_S_new = [S_S_new;S_S_spec]; kb_names_new(nr_line+1:nr_line+length(I)) = kb_names_spec; % === Calculate the correlation matrix. The idea is to assume one % correlation coefficient for one species; % ---------------------------------------------------------------- out(2,['*** Doing for ', name_spec{i}, ' ***']); out(2,['===============================']); out(3,['=== constructing the corelation matrix for ', name_spec{i}, ' ===']); Mcorr_spec= ones(length(I), length(I))* Q.CORRELATION(i); correlation(ni) = Q.CORRELATION(i); for ij = 1: length(I) Mcorr_spec(ij,ij) =1; end out(3,'=== filling into full correlation matrix ==='); Mcorr(nr_line+1:nr_line+length(I), nr_line+1:nr_line+length(I)) = Mcorr_spec; nr_line =nr_line + length(I); out(3,['=== calculate the total error for ', name_spec{i}, ' ===']); % calculate the covariance matrix just for one molecular % species for m=1:size(Mcorr_spec,1) for ij=1:size(Mcorr_spec,2) S_S_spec_new(m,ij)=Mcorr_spec(m,ij)*S_S_spec(m)*S_S_spec(ij); end end clear S_S_spec; S_S_spec = S_S_spec_new; clear S_S_spec_new; % calculate the total error for the species in question out(3,['=== calculate the total error for ', name_spec{i}, ' ===']); S_spec_tot{ni} =Dy* Kb_spec; S_spec_tot{ni} = S_spec_tot{ni}*S_S_spec*S_spec_tot{ni}'; S_spec_tot{ni} = sqrt(diag(S_spec_tot{ni})); name_spec_real{ni} = name_spec{i}; end clear kb_names_spec S_S_spec Kb_spec; end nr_ss = length(S_spec_tot); % number of species for which weighting function is calculated. % This could differ from wtgss is there species present in wtgss but o lines for the corresponding % species were found in the catalog % calculate the total error for all species; !!!just for a check. S_tot = S_spec_tot{1}.^2; % initialize for i =2 : nr_ss S_tot =(S_tot+S_spec_tot{i}.^2); end S_tot = sqrt(S_tot); % clear some variables clear S_S_param Kb_param kb_names_param ; clear name_spec2 spec name_spec S_S_param = S_S_new; Kb_param = Kb_new; kb_names_param = kb_names_new; clear S_S_new Kb_new kb_names_new ; % calculate the error for individual parameters and lines for i=1:npss S_Si_param{i} = Dy*Kb_param(:,i); S_Si_param{i} = S_Si_param{i} * S_S_param(i).^2 * S_Si_param {i}'; S_Si_param{i} = sqrt(diag(S_Si_param{i})); end % calculating the full covariance matrix for one parameter (intensity, % position, agam, ...) for m=1:size(Mcorr,1) for ij=1:size(Mcorr,2) S_S_param_new(m,ij)=Mcorr(m,ij)*S_S_param(m)*S_S_param(ij); end end clear S_S_param; S_S_param = S_S_param_new; clear S_S_param_new; % calculate the total error for one parameter S_S_tot_param = Dy*Kb_param; S_S_tot_param = S_S_tot_param *S_S_param* S_S_tot_param'; S_S_tot_param = sqrt(diag(S_S_tot_param)); % Making plots out(1, 'Plotting ...'); [line_style, line_color, line_width]= plinecolor(10); %=== Set axes position A_POS = [0.08 0.2 0.4 0.55 0.48 .2 0.5 0.55]; %=== First figure - 1 if isempty(get(0,'Children')) % To handle the case when no window fig = 0; % is opened else fig = gcf; end %=== Loop the retrieval identities ============================================ % for ii = 1:nrq %= Names of present retrieval identity [group,name] = qp_kinfo( kx_names{ii} ); %= Index in x of present retrieval identity ind = kx_index(ii,1):kx_index(ii,2); %= Get index for points with measurement response above limit ip = find( measres(ind) >= mr_limit ); ip =ind; %= Sort the spectroscopic error in decreasing order % % ====== sort for the uncorrelated case ==== for ipar=1:npss dummy(ipar)=max(S_Si_param{ipar}(ind)); end [dummy, I]=sort(dummy); I=fliplr(I); clear dummy; for ipar=1:npss kb_names_param1{ipar}= kb_names_param{I(ipar)}; S_Si_param_new{ipar}(ind) = S_Si_param{I(ipar)}(ind); end %= 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 for iname=1:length(name) if name(iname)==' '|name(iname)==','|name(iname)=='/' name(iname)='_'; end end s = strread(name, '%s', 'delimiter', '-'); if length(s)>1 name = [char(s(1)),'-', char(s(2))] ; %take only species and isotope end %=== Scalar identity %=========================== % creates a table tex if length( ip ) == 1 sotname = [Q.SAVE,'-', deblank(name),'_ALL_', ... num2str(j), '.tex']; SOT=fopen(sotname,'wt'); fprintf(SOT,'\\documentclass[12pt,leqno,twoside,a4paper]{article}\n'); fprintf(SOT,'\\usepackage[]{epsf} \n'); fprintf(SOT,'\\begin{document}\n') fprintf(SOT,'\\begin{table}\n'); fprintf(SOT,'\\begin{tabular}{|ll|} \\hline\n'); round(za2ztan(Q,s_tot(ind)*xfac)) fprintf(SOT,'Retrieval precision& %5.0f \\\\\n', round(za2ztan(Q,s_tot(ind)*xfac))); fprintf(SOT,'Total: %s & %5.0f\\\\\n', ... name_param{j},round(za2ztan(Q,S_S_tot_param(ind)*xfac))); for ipar=1:npss if ipar <=20 fprintf(SOT,' %s & %5.0f\\\\\n', kb_names_param1{ipar}, ... round(za2ztan(Q, S_Si_param_new{ipar}(ind)*xfac))); end end fprintf(SOT,'\\end{tabular}\n'); fprintf(SOT,'\\end{table}\n'); fprintf(SOT,'\\end{document}\n'); fclose(SOT) %=== 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( fullfile( 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'; % elseif strcmp( name, 'Reference load temperatures' ) % z(ind) = Q.TB_REFLOADS_NOMINAL; % yfac = 1; ytext = 'Nominal value [K]'; yunit = 'K'; xfac = 1; xunit = 'K'; % else error(sprintf( 'Unknown retrieval identity: %s', name )); end for iname=1:length(name) if name(iname)==' '|name(iname)==','|name(iname)=='/' name(iname)='_'; end end % plot the parameters for all lines in decresing order %====================================================== out(3, 'plot the parameters for all lines in decresing order'); fig=fig + 1; figure(fig); hold off; orient landscape; subplot(2,2,3); % % plot the total spectro error out(3, 'plot the total erorr'); plot(S_S_tot_param(ind)*xfac, z(ind)/yfac, '--k', 'linewidth', 2) hold on; line_name{1}=['Total: ', name_param{j}]; % ni=2; % plot the individual error out(3, 'plot the individual error'); for ipar=1:npss if ni<=20 %& max(S_Si_param_new{ipar}(ip))/max(S_S_tot_param(ip))>=0.02 plot(S_Si_param_new{ipar}(ip)*xfac, z(ip)/ ... yfac,line_style{ni-1}, 'color', line_color{ni-1},'linewidth', ... line_width{ni-1}, 'markersize', 12); hold on line_name{ni}= kb_names_param1{ni-1}; ni=ni+1; end end %------------- ylim([5,60]); hold on; % set axis %---------- % 1. for the spectro parameters errors ax1 = gca; set(ax1,'XColor','k','YColor','k','fontsize', 14 ); xlabel(['error [',xunit,'] '], 'fontsize', 15, 'fontweight', 'demi'); ylabel(' altitude [km]', 'fontsize', 15, 'fontweight', 'demi'); set(ax1, 'FontSize',14, 'fontweight', 'normal', 'POSITION', ... A_POS(1, :)); % 2. for the retrieval precision ax2 = axes('Position',get(ax1,'Position'), 'XAxisLocation', 'top', ... 'YAxisLocation', 'right', 'Color', 'none', ... 'XColor', 'r', 'YColor', 'k'); % plot retrieval precision plot( s_tot(ip)*xfac, z(ip)/yfac, '*-r','linewidth', 1 ); line_name{length(line_name)+1}=['Retrieval Precision ']; ylim([5,60]); set(ax2, 'XAxisLocation','top', 'YAxisLocation', 'right', 'Color','none',... 'Position', get(ax1,'Position'),... 'XColor','k', 'YColor','k', 'FontSize',14, ... 'fontweight', 'normal', 'YTicklabel', []); title([upper(band), ': ', upper(name_param{j}), ': ' ... ,' error on, ', name], 'fontsize', 16, 'FontWeight', 'bold'); % make legend %============= yp=.99; subplot(2,2,4) %------------- tt=line_name{1}; h=plot([0.05 0.18],[yp yp], '--k', 'linewidth', 2); hold on; text(.23,yp,tt,'FontSize',15, 'color','k') % for ni=2:length(line_name)-1 tt=line_name{ni}; yp=yp-.048; h=plot([0.05 0.18],[yp yp],... line_style{ni-1}, 'color', line_color{ni-1}, ... 'linewidth', line_width{ni-1}, 'markersize', 12); text(.23,yp,tt,'FontSize',15, 'color',line_color{ni-1}); end yp=yp-.048; tt=line_name{length(line_name)}; h=plot([0.05 0.18],[yp yp], '*-r', 'linewidth', 1); text(.23,yp,tt,'FontSize',15, 'color', 'r'); set(gca, 'fontsize', 15); % axis invisible axis([0 1 0 1]) set(gca,'Box','off'); set(gca,'YTickLabel',[]); set(gca,'XTickLabel',[]); set(gca,'YColor',[1 1 1]); set(gca,'XColor',[1 1 1]); set(gca,'FontSize',15, 'POSITION', A_POS(2,:)) axis off; %save figures %============== strsave=[ 'print -depsc ', Q.SAVE,'-', deblank(name), '-'... deblank(name_param{j}),'.eps']; eval(strsave); line_name={''}; % plot the total spectro error for each species fig=fig + 1; figure(fig); hold off; orient landscape; % plot the total spectro error out(3, 'plot the total erorr'); plot(S_tot(ind)*xfac, z(ind)/yfac, '--k', 'linewidth', 2) hold on; line_name{1}=['Total: ', name_param{j}]; % for ni =1:nr_ss plot(S_spec_tot{ni}(ind)*xfac, z(ind)/yfac, line_style{ni}, 'color', line_color{ni}, ... 'linewidth', line_width{ni}, 'markersize', 12); hold on; line_name{ni+1}=['Total - ', name_spec_real{ni}, ': ', num2str(correlation(ni)), ' corr']; end ax1 = gca; set(ax1,'XColor','k','YColor','k','fontsize', 14 ); xlabel(['error [',xunit,'] '], 'fontsize', 15, 'fontweight', 'demi'); ylabel(' altitude [km]', 'fontsize', 15, 'fontweight', 'demi'); set(ax1, 'FontSize',14, 'fontweight', 'normal', 'POSITION', ... A_POS(1, :)); % make legend %============= yp=.99; subplot(2,2,4) %------------- %------------- tt=line_name{1}; h=plot([0.05 0.18],[yp yp], '--k', 'linewidth', 2); hold on; text(.23,yp,tt,'FontSize',15, 'color','k'); for ni=2:length(line_name) tt=line_name{ni}; yp=yp-.05; h=plot([0.05 0.18],[yp yp],... line_style{ni-1}, 'color', line_color{ni-1}, ... 'linewidth', line_width{ni-1}, 'markersize', 12); hold on; text(.23,yp,tt,'FontSize',15, 'color',line_color{ni-1}); end % set axis invisible axis([0 1 0 1]) set(gca,'Box','off'); set(gca,'YTickLabel',[]); set(gca,'XTickLabel',[]); set(gca,'YColor',[1 1 1]); set(gca,'XColor',[1 1 1]); set(gca,'FontSize',15, 'POSITION', A_POS(2,:)) axis off; %save figures %============== strsave=[ 'print -depsc ', Q.SAVE,'-', deblank(name), '-', deblank(name_param{j}), ... '_pwtgss.eps']; eval(strsave); line_name={''}; end end end end %============================================================================== %=== %=== Sub-functions %=== %============================================================================== % converts the deviation in zenith angle into deviation in tangent % altitudes function dztan=za2ztan(Q, dza) global EARTH_RADIUS EARTH_RADIUS = 6378e3 ztan1 = 1e3; za1 = geomtan2za(EARTH_RADIUS, Q.PLATFORM_ALTITUDE, ztan1); za2 = za1+ dza; ztan2 = za2geomtan(EARTH_RADIUS, Q.PLATFORM_ALTITUDE, za2); dztan=abs(ztan2 - ztan1); function A = axes_def() %02; A.BOX = 'on'; % A.FONTSIZE = 15; A.FONTWEIGHT = 'normal'; % A.XLABEL_SIZE = 15; A.XLABEL_WEIGHT = 'bold'; % A.YLABEL_SIZE = 15; A.YLABEL_WEIGHT = 'bold'; %