D_visualizePerturbedPatterns.m

%{
D_visualizePerturbedPatterns.m

this script computes and overlays patterns from the nominal array
excitation and several monte carlo variations to visualize the relative
deviation in perturbed patterns.  the ensemble mean directivity is also
plotted at each angle.  

see E_generalDirectivityEstimate for detailed comparison with directivity
estimates from the paper.

%}

close all
clear
clc

% load data
% -------------------------------------------------------------------------
tag = 'MeshedHalfWaveDipolesHFSS';      % choose either HFSS or ATOM
load(['data/opdata-',tag,'.mat'])
load(['data/pareto-',tag,'.mat'])
load(['data/montecarlo-',tag,'.mat'])


% constants and settings
% -------------------------------------------------------------------------
eta0 = 376.730313412;   % free space impedance
Nplot = 300;            % number of monte carlo variations to plot

% unpack relevant data
% -------------------------------------------------------------------------
OP = OPsweep{1,1};
F = OP.FthPatt;
R0 = OP.R0;
Iopt = Icell{1,1};
mc = mc{1,1};

% cycle through set of nominal directivity solutions
% -------------------------------------------------------------------------
for k = 1:length(D0List)

    % compute and plot nominal pattern
    % ---------------------------------------------------------------------
    x0 = Iopt(:,idx(k));
    U = abs(F*x0).^2/(2*eta0);
    P = real(x0'*R0*x0)/2;
    D = 4*pi*U/P;
    figure()
    a1 = plot(angleVec,10*log10(D),'linewidth',2);
    hold on

    % compute and plot monte carlo patterns + ensemble mean
    % ---------------------------------------------------------------------
    DS = 0;
    for n = 1:Nplot
        x0 = mc.IMC(:,n,k);
        U = abs(F*x0).^2/(2*eta0);
        P = real(x0'*R0*x0)/2;
        D = 4*pi*U/P;
        a2 = plot(angleVec,10*log10(D),'color',[1,0,0,0.1]);
        DS = DS+D;
    end
    DS = DS/Nplot;
    a3 = plot(angleVec,10*log10(DS),'k','LineWidth',2);
    xlim([min(angleVec),max(angleVec)])
    ylim([-30,30])
    legend([a1,a2,a3],{'D(I_0)','D(I+\deltaI)','\langleD\rangle'},'location','best')
    xlabel('\phi [rad]')
    ylabel('D [dB]')
    title(['D_0 = ',num2str(10*log10(D0List(k))),' dB'])
end