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

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

    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