How to plot average of logarithm of transmission versus sample length by positional-randomize structure?

Question

GULZAR 2024-3-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2091516-how-to-plot-average-of-logarithm-of-transmission-versus-sample-length-by-positional-randomize-struct

编辑： GULZAR 2024-3-21

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program uses Transfer Matrix Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization of parameters
wavelength0=1.5e-6;
nA = 1.1;  % Refractive index of material A
nB = 3.5;   % Refractive index of material B
dA = wavelength0/(4*(nA));   % Thickness of layer A in meters
dB = wavelength0/(4*(nB));   % Thickness of layer B in meters
% Define parameters
num_samples = 100;  % Number of disorder configurations
L_max = 220e-6; % Maximum sample length in meters
dL = 2e-6;     % Length increment in meters
n0=1;                               % Refractive index of Incident and transmission medium
c=3e8;                              % Velocity of Light (m/s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% wavelength Range
lambda=linspace(1490,1510,10);                    % wavelength in nm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Preallocate for speed
log_transmissions = zeros(num_samples, floor(L_max/dL) + 1);
% Loop over sample lengths
for L = dL:dL:L_max
    num_layers = round(L / (dA + dB)); % Determine number of layers
    % Loop over disorder realizations
    for sample = 1:num_samples
        % Generate random layer sequence
        layer_order = randi([1, 2], 1, num_layers);
        % Initialize transfer matrix
        N = [1 0; 0 1];
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for loop=1:length(lambda)
            f=(c*1e9./lambda);
            % angular frequency
            w=2*pi*f;
            a=dA+dB;                                         % Lattice Constant
            DA=((w(loop))/c)*dA*nA;         % phase angle in first layer
            DB=((w(loop))/c)*dB*nB;         % phase angle in second layer
            g0=n0;
            gA=nA;
            gB=nB;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Transfer Matrix elements of first layer
            m11=cos(DA);
            m12=-1i*sin(DA)/gA;
            m21=-1i*gA*sin(DA);
            m22=cos(DA);
            m1=[m11 m12; m21 m22];
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Transfer MAtrix elements of Second layer
            l11=cos(DB);
            l12=-1i*sin(DB)/gB;
            l21=-1i*gB*sin(DB);
            l22=cos(DB);
            m2=[l11 l12; l21 l22];
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Build transfer matrix
            for j = 1:num_layers
                if layer_order(j) == 1
                    N = N *m1 ;
                else
                    N = N *m2;
                end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Transmission
            T1(loop)=(2*g0/(g0*N(1,1)+g0*g0*N(1,2)+N(2,1)+g0*N(2,2)));            
            T(loop)=(T1(loop)*conj(T1(loop)));
            % Store log transmission
            log_transmissions(sample, round(L/dL) + 1) = log10(T1(loop));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
    end
end
toc
Elapsed time is 10.546770 seconds.
% Calculate ensemble average
    avgT = mean(log_transmissions, 1);
    Ls = linspace(dL, L_max, length(avgT));
figure(1)
plot(Ls*1e6,real(avgT(:)),'o','LineWidth',1.5)
xlim([0,250])