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 THz\n', 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 THz\n', 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.