fresnelSynthNuMaxVarPublic
Compares Institute Fresnel's data set with synthetic data generated by the exact contrast for several approximation orders nuMax. This is useful to find an appropriate approximation order nuMax for the incident field's matching.
Warning: This function clears variables and creates a directory. Be careful when using.
Contents
Syntax
Start this routine in the folder code by
init; fresnelSynthPaper
Description
This routine is similar to fresnelSynthPaper. Therefore we refer to fresnelSynthPaper.html for a description.
In contrast to the routine fresnelSynthPaper computations with various values for the approximation order nuMax are supported. Therefore the parameter nuMaxMax can be set.
In addition, resulting figures are saved in output/fresnelSynthNuMaxVarPublic.
Note that a suitable choice of nuMax is essentially for a suitable matching of the incident field.
More About
This file is essentially a copy of the intern routine fresnelSynthNuMaxVar with adaptions and documentation. It was used to create Fig. 3.9 in [1]. (The routine fresnelSynthNuMaxVar originated from fresnelSynth. It is a small part, that was adapted.)
References
- [1] Florian Bürgel. Effective and Efficient Reconstruction Schemes for the Inverse Medium Problem in Scattering. PhD thesis, Universität Bremen, 2019.
See Also
Code
Initialization
close all; % clearvars -except inseti htext; clear all; disp(' ') disp('-----') disp('fresnelSynth') disp(' ') % -- initialize tic; init; calcyExact = 0; % Make dir to save figures dirOutput = 'output/fresnelSynthNuMaxVarPublic/'; if exist(dirOutput,'dir')~=7 fprintf('Make dir: %s.\n',dirOutput) mkdir(dirOutput); end % -- set parameters % Two possibilities: % * Use "-- load paperFresnelSynth" from start to end. % or % * Use "-- Set parameters from paperFresnelSynth (is a copy)" from start to end. % -- load paperFresnelSynth: start -- % if ~exist('inseti','var') % %inseti = 'setiContrastFresnelSynth'; % synthetic fresnel data % inseti = 'paperFresnelSynth'; % synthetic fresnel data % end % seti.inseti = inseti; % eval(inseti) % -- load paperFresnelSynth: end -- % -- Set parameters from paperFresnelSynth (is a copy): start -- % -- setGrid seti.dim = 2; seti.rCD = 0.2; seti.nCD = 256; % original % -- fresnel data (loadData.m) seti.expData = 'fresnel'; seti.fresnelFreq = 5*1E9; % standard 3 GHz (4 GHz) seti.fresnelFile = 'inexpdata/fresnel_opus_1/dielTM_dec8f.exp'; % 1--8 GHz, test: single cylinder, 4 GHz %seti.fresnelFile = 'data/fresnel_opus_1/twodielTM_8f.exp'; % 1--8 GHz filename = seti.fresnelFile; seti.nuMax = 10; seti.ampCalc = 1; % default (B1) % -- setKernel cLight = 2.99792458E8; % light velocity seti.k = 2*pi*seti.fresnelFreq/cLight; seti.model = 'helmholtz'; % -- setMeas (measurements) % simulate fresnel_opus_1: seti.incType = 'pointSource'; % used in op. 1 and 2 (correct? should be...) seti.measType = 'nearField'; % used in op. 1 and 2 seti.incNb = 36; seti.measNb = 72; seti.radSrc = 0.72; seti.radMeas = 0.76; % -- inversion process (regularization) %seti.invNo = 3; seti.invNo = 6; seti.delta = 0; % to test %seti = setiAuxRec(seti,usefindpdatau,usefindalpha); seti.alpha = 1E8; % old: 2E4 corresponding 0.05 (in old code) seti.beta = 0; % because we expect no sharp edges... seti.physBounds = [-500,500,0,inf]; seti.nOut = 30; % -- minimization with pda seti.pdaN = 50; % -- Set parameters from paperFresnelSynth (is a copy): end --
Load Data from Institute Fresnel
seti = checkConsisExpData(seti); for freq = 1:8 % INPUT: frequency in GHz
seti.fresnelFreq = freq*1E9; nuMaxMax = 20; % INPUT: nuMaxMax % nuMaxMax = 100; % INPUT: nuMaxMax errCmean = zeros(nuMaxMax,1); errRe = zeros(nuMaxMax,1); errIm = zeros(nuMaxMax,1); for nuMax = 1:nuMaxMax
disp('----------------------------------') fprintf('nuMax = %g\n',nuMax); disp('----------------------------------') seti.nuMax = nuMax; disp(['loading fresnel data from "' seti.fresnelFile '" for frequency = ' num2str(seti.fresnelFreq/1E9) ' GHz']) ; seti = loadData(seti); %seti = evalc('loadData(seti)'); % setGeomSim is called in loadData %evalc: computes but no disp output... %seti.errC is inside... yFresnel = seti.FmeasDelta; % is loaded Fresnel data %disp('- yExactInc calculation using incField from Fresnel data') %disp(' but incField is NOT exact, so yExactInc might not be Fresnel data') FmeasExactInc = mimo(seti, seti.qROIexact, 'simo'); FmeasExactInc = fresnelMiss(FmeasExactInc,yFresnel);
Compare yFresnel and yExactInc
Compare real-world data from Institute Fresnel, yFresnel, with exact computed data, yExactInc.
disp('Compare yFresnel and yExactInc') errCmean(nuMax) = mean(seti.errC); errRe(nuMax) = norm(real(yFresnel-FmeasExactInc))/norm(real(FmeasExactInc)); errIm(nuMax) = norm(imag(yFresnel-FmeasExactInc))/norm(imag(FmeasExactInc)); %fprintf('nuMax = %g | ',nuMax); %fprintf('rel. err (real part) = %g | ',errRe); %fprintf('rel. err (imag part) = %g\n',errIm);
end
Notation of Values
- Cmean: rel. err. of
(mean because min = max in our experience)
- errRe: rel. err of data y (real part)
- errIm: rel. err of data y (imag part)
disp('----------') disp(' nuMax | Cmean in % | errRe in % | errIm in %'); % with *100 below... then %... for nuMax = 1:nuMaxMax fprintf('%d & %4.1f & %4.1f & %4.1f\n',nuMax,errCmean(nuMax)*100, errRe(nuMax)*100, errIm(nuMax)*100); end disp('----------')
Figures
figFontSize = 20; % default in ipscatt publish: 20 (pubFontSize) (does influence the result; tested with 40) % lw = 2; % LineWidth (default for printing; is very thin in PhD thesis) lw = 3; % LineWidth (default for talks) figure(101); % h = plot(t,step,'b-',t,stepinv,'r:',t,t,'k-.'); % h = plot(1:nuMaxMax, errCmean, 1:nuMaxMax, errRe, 1:nuMaxMax, errIm); h = plot(1:nuMaxMax,errCmean,'k-.',... 1:nuMaxMax,errRe,'r-',... 1:nuMaxMax, errIm,'b:'); % -. dash-dot line; - solid line, : dotted line set(h(1),'LineWidth',lw); set(h(2),'LineWidth',lw); set(h(3),'LineWidth',lw); axis square; % xlim([-5 5]); ylim([0 2]); set(gca, 'FontSize', figFontSize) % set(gca, 'XTick', [a 0 b]); % set(gca, 'YTick', [a 0 b]); % set(gca, 'XTickLabelMode', 'manual', 'XTickLabel', {'a','0','b'}); % set(gca, 'YTickLabelMode', 'manual', 'YTickLabel', {'a','0','b'}); % legend('mean(errC)','data errRe','data errIm'); % xlabel('nuMax'); % title('suitable amplitudes: variation of nuMax')
Save Figures as png, eps and fig
seti.dirDatetime = datestr(now, 'yyyymmddTHHMMSS'); % additional information [~,name,~] = fileparts(seti.fresnelFile); filePng = sprintf('%s%s_%d_GHz_%s.png',dirOutput,name,seti.fresnelFreq/1E9,seti.dirDatetime); figNo = 101; figNostr = sprintf('-f%d',figNo); print(figNostr,'-dpng',filePng); fileEpsc = strrep(filePng,'.png','.eps'); print(figNostr,'-depsc',fileEpsc); fileFig = strrep(filePng,'.png','.fig'); savefig(figNo,fileFig);
end