TMM Multilayer structure Transmission

% Define constants
c0 = 3e8; % Speed of light in vacuum (m/s)
dA = 30e-6; % Thickness of SiO2 layer in meters
dB = 7e-6; % Thickness of Si layer in meters
dD1 = 2e-6; % Thickness of D1 layer in meters
dD0 = 12e-6; % Thickness of D0 layer in meters
dD2 = 2e-6; % Thickness of D2 layer in meters
N = 5; % Number of SiO2 and Si periods
nA = 2.1; % Refractive index of SiO2
nB = 3.4; % Refractive index of Si

% Define frequency range
frequencies = linspace(1, 2.5, 1000); % Frequency range from 1 THz to 2.5 THz

% Define refractive index (n) and absorption coefficient (alpha) for 0% PG (water)
n_water = -0.07042 * frequencies.^3 + 0.4293 * frequencies.^2 – 0.9245 * frequencies + 2.744;
alpha_water = -17.0559 * frequencies.^2 + 151.258 * frequencies + 81.2348;

% Compute complex refractive index for defect layer (D0)
n_complex_water = n_water + 1i * (c0 * alpha_water) ./ (4 * pi * frequencies * 1e12);

% Define other refractive indices for layers
nD1 = nA; % Refractive index of D1 (SiO2)
nD2 = nA; % Refractive index of D2 (SiO2)
nS = nA; % Refractive index of substrate (SiO2)
n0 = 1; % Refractive index of ambient medium
theta0 = 0; % Angle of incidence in air (normal incidence)

% Calculate f_PBG using Equation 26
f_PBG = (c0 * (nA + nB)) / (4 * (dA + dB) * (nA * nB));

% Display f_PBG result
fprintf(‘The photonic bandgap frequency (f_PBG) is %.2f THzn’, f_PBG / 1e12);

% Calculate f_R using Equation 27
m = 1; % Order
f_R = m * c0 / (2 * dD0 * mean(real(n_complex_water)));

% Display f_R result
fprintf(‘The resonance frequency (f_R) is %.2f THzn’, f_R / 1e12);

% Initialize results
transmittance = zeros(length(frequencies), 1);

% Loop over each frequency
for j = 1:length(frequencies)
f = frequencies(j) * 1e12; % Convert frequency to Hz
omega = 2 * pi * f;
lambda = c0 / f;

% Calculate incident angles using Snell’s Law
thetaA = asin(n0 / nA * sin(theta0)); % Angle in SiO2 layer
thetaB = asin(n0 / nB * sin(theta0)); % Angle in Si layer
thetaD1 = asin(n0 / nD1 * sin(theta0)); % Angle in D1 layer
thetaD0 = asin(n0 / real(n_complex_water(j)) * sin(theta0)); % Angle in D0 layer (using real part)
thetaD2 = asin(n0 / nD2 * sin(theta0)); % Angle in D2 layer

% Calculate Sigma using Equation 20
sigmaA = (2 * pi * dA * nA * cos(thetaA)) / lambda;
sigmaB = (2 * pi * dB * nB * cos(thetaB)) / lambda;
sigmaD1 = (2 * pi * dD1 * nD1 * cos(thetaD1)) / lambda;
sigmaD0 = (2 * pi * dD0 * real(n_complex_water(j)) * cos(thetaD0)) / lambda; % Using real part of refractive index
sigmaD2 = (2 * pi * dD2 * nD2 * cos(thetaD2)) / lambda;

phiA = nA * cos(thetaA);
phiB = nB * cos(thetaB);
phiD1 = nD1 * cos(thetaD1);
phiD0 = real(n_complex_water(j)) * cos(thetaD0); % Using real part of refractive index
phiD2 = nD2 * cos(thetaD2);

% Transfer matrices
aA = [cos(sigmaA), -1i / phiA * sin(sigmaA); -1i * phiA * sin(sigmaA), cos(sigmaA)];
aB = [cos(sigmaB), -1i / phiB * sin(sigmaB); -1i * phiB * sin(sigmaB), cos(sigmaB)];
aD1 = [cos(sigmaD1), -1i / phiD1 * sin(sigmaD1); -1i * phiD1 * sin(sigmaD1), cos(sigmaD1)];
aD0 = [cos(sigmaD0), -1i / phiD0 * sin(sigmaD0); -1i * phiD0 * sin(sigmaD0), cos(sigmaD0)];
aD2 = [cos(sigmaD2), -1i / phiD2 * sin(sigmaD2); -1i * phiD2 * sin(sigmaD2), cos(sigmaD2)];

% Total matrix
M_total = (aA * aB)^N * aD1 * aD0 * aD2 * (aA * aB)^N;

% Calculate admittance values
phi0 = n0 * cos(theta0);
phis = nS * cos(thetaA);

% Calculate transmittance using the provided formula
t = (2 * phi0) / ((M_total(1,1) + M_total(1,2) * phis) * phi0 + (M_total(2,1) + M_total(2,2) * phis));
T = phis / phi0 * abs(t)^2 * 100;
transmittance(j) = T;
end

% Plot results
figure;
plot(frequencies, transmittance); % Already in percentage
xlabel(‘Frequency (THz)’);
ylabel(‘Transmittance (%)’);
title(‘Transmittance of THz PG Sensor vs. Frequency’);

The output should look like the uploaded image. Check the code for mistakes. The amplitude of T (in %) at the cental frequency doesn’t match with the image.% Define constants
c0 = 3e8; % Speed of light in vacuum (m/s)
dA = 30e-6; % Thickness of SiO2 layer in meters
dB = 7e-6; % Thickness of Si layer in meters
dD1 = 2e-6; % Thickness of D1 layer in meters
dD0 = 12e-6; % Thickness of D0 layer in meters
dD2 = 2e-6; % Thickness of D2 layer in meters
N = 5; % Number of SiO2 and Si periods
nA = 2.1; % Refractive index of SiO2
nB = 3.4; % Refractive index of Si