%------------------------------------------------------------------------------ % NAME: qpopt_fmono % % This function optimises the monochromatic frequency grid based % on the polynomial order selected for the spectra (F_ORDER). % % The optimisation is performed in this way: % 1. Spectra are calculated for the set-up atmospheres for % a very fine F_MONO. % 2. The end points of each frequency range are included in % F_MONO. % 3. The spectra are interpolated from the present F_MONO % to the fine test grid. Interpolation applied follows F_ORDER. % 4. The maximum difference between interpolated and calculated % spectra (for the fine grid) is determined. % 5. If the difference exceeds the given precision, the frequency % where the maximum difference was found is put into F_MONO % 6. Points 3-5 are repeated the wanted precision is obtained. % % The procedure is repeated for each frequency range specified by % FRANGES. The precision must be set for each frequency range. % A typical example would be that FRANGES specifies two ranges, % where one range is the primary band and the other range is the % image band. % % The bands can be specified in any order, but the final F_MONO is % sorted in ascending order. % % If the function is called without output argument, F_MONO is stored % to the corresponding file specified by Q. % % FORMAT: f_mono = qpopt_fmono(Q,franges,prec,za_pencil % [, fminspace, do_plot ]) % % OUT: f_mono Selected monochromatic frequency grid. % IN: Q Setting structure. % franges 2-column marix, where column 1 gives the lower % end point of each frequency range and column 2 % gives the upper end point. % prec Precision for each frequency range. This is a % vector with the same length as number of rows % of FRANGES. % za_pencil Vector with zenith angles to be used in the % function. % OPTIONAL: fminspace Spacing of the test frequency grids. Default is % 20% of the Doppler width for a very heavy molecule. % The default is used if FMINSPACE=[]. % do_plot Flag to show example on position of found frequency % grid and interpolation error. Default is 0 (meaning % no plotting). %------------------------------------------------------------------------------ % HISTORY: 2001.11.20 Created by Patrick Eriksson. function f_mono = qpopt_fmono(Q,franges,prec,za_pencil,fminspace,do_plot) %= Check F_ORDER % if Q.F_ORDER~=1 & Q.F_ORDER~=3 error('The Q-filed F_ORDER must be 1 or 3.'); end %= Number of frequency ranges % nranges = size( franges, 1 ); % if length( prec) ~= nranges error('The size of FRANGES and the length of PREC do not match'); end %= Determine of minimum frequency spacing, if not given %= The spacing is set to 1/5 of the Doppler width for a molecule with %= molecular weight of 100 (at a temperature of 200 K). %= (For calculation of the Doppler width, see my thesis, page 221). % if ~exist('fminspace') | isempty(fminspace) fmean = (min(franges(:,1))+max(franges(:,2))) / 2; fminspace = 7.16e-7 * fmean * sqrt(2) / 5; end out(1,1); out(1,'Optimising F_MONO with respect to interpolation error.'); if Q.F_ORDER == 1 out(1,'Polynomial representation: linear'); else out(1,'Polynomial representation: cubic'); end if exist('do_plot') & do_plot out(1,'Plotting will be performed.'); end out(1,sprintf('Spacing of test grid: %.0f kHz',fminspace/1e3)); out(2,0); %=== Loop ranges and determine f_mono for the range % f_mono = []; % nza = length( za_pencil ); % for i = 1:nranges out(2,sprintf('Doing frequency range %d/%d',i,nranges)); out(2,sprintf('Required precision %.3f K',prec(i))); out(3,'Calculating reference spectra ...'); %= Create test frequency grid % f1 = franges( i, 1 ); f2 = franges( i, 2 ); nf = max([ ceil((f2 - f1)./fminspace+0.0001), 2 ]); % ftest = linspace( f1, f2, nf )'; %= Calculate test spectra for the set-up atmospheres % nsu = sstring_length( Q.SETUP_VMR ); % Q2 = Q; % Y = zeros( nf, nza, nsu ); % for j = 1:nsu % out(3,sprintf('Reference atmosphere %d/%d',j,nsu)); % Q2.APRIORI_VMR = sstring_get_i( Q.SETUP_VMR, j ); Q2.APRIORI_PTZ = sstring_get_i( Q.SETUP_PTZ, j ); % y = qp_ympb( Q2, ftest, za_pencil ); % Y(:,:,j) = reshape( y, nf, nza ); plot(Y) pause % end % if Q.F_ORDER == 1 ind = 1; else ind = [1;round(nf/2)]; end % imax = nf; % emax = 1e99; % while emax > prec(i) % ind = [ind;imax]; ind = sort( ind ); % emax = 0; % for j = 1:nsu if Q.F_ORDER == 1 Y2 = interp1( ftest(ind), Y(ind,:,j), ftest, 'linear' ); else Y2 = interp1( ftest(ind), Y(ind,:,j), ftest, 'cubic' ); end [ej,ij] = max( max( abs(Y2-Y(:,:,j)), [], 2 ) ); % if ej > emax emax = ej; imax = ij; end end % if emax > prec(i) out(3,sprintf('Error=%.3f, adding %.6f GHz, nr=%d',... emax,ftest(imax)/1e9,length(ind)+1)); else out(3,sprintf('Error=%.3f',emax)); end end %= Show results if DO_PLOT % if exist('do_plot') & do_plot % i1 = ceil(nza/2); % figure(1), clf plot( ftest/1e9, Y(:,i1,nsu), '-', ftest(ind)/1e9, Y(ind,i1,nsu), '.' ) xlabel('Frequency [GHz]'); ylabel('Tb [K]') title( sprintf( ... 'Last test atmosphere, zenith angle = %.2f deg',za_pencil(i1)) ); % figure(2), clf plot( ftest/1e9, (Y2(:,i1) - Y(:,i1,nsu))*1e3, '-' ) xlabel('Frequency [GHz]'); ylabel('Interpolation error [mK]') title( sprintf( ... 'Last test atmosphere, zenith angle = %.2f deg',za_pencil(i1)) ); % out(1,'Press a key to continue.'); pause; end f_mono = [ f_mono; ftest(ind) ]; end %=== Sort and remove duplicates of frequencies % f_mono = unique( f_mono ); %=== Save if no output argument % if ~nargout % ss = 'Optimized by QPOPT_FMONO_BY_INTERP'; ss = str2mat(ss,''); for i = 1:nranges s = sprintf('Range %d',i); ss = str2mat(ss,s); s = sprintf('%.6f-%.6f GHz, precision %.3f',... franges(i,1)/1e9,franges(i,2)/1e9,prec(i)); ss = str2mat(ss,s); end % write_datafile( fullfile(Q.CALCGRIDS_DIR,Q.F_MONO), f_mono, 'VECTOR',10, ss); % end out(1,-1);