%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a program for 1-D Photonic crystal
% For Photonic band structure calculation
% This program uses Transfer Matrix Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters
P = [0 50 100 150 200]; %%% Hydrostatic pressure in MPa
n0 = 1; % Refractive index of air
ns = 1.46; % Refractive index of subtrate
n01 = 1.578; % Refractive index of PS
n02 = 1.484; % Refractive index of PMMA
% Calculate Epsilon for n01, n02
e01 = n01^2; e02 = n02^2;
% Constants
p11_1 = 0.32; p12_1 = 0.31; p11_2 = 0.3; p12_2 = 0.297;
v1 = 0.35;
E1 = 3.3e3; % Convert to scientific notation (3.3 x 10^3)
v2 = 0.37;
E2 = 3.0303e3; % Convert to scientific notation (3.0303 x 10^3)
% Formula for calculating nA
% Initialization of parameters
dA=780e-9; % Thickness of First Layer in meters
dB=830e-9; % Thickness of Second Layer in meters
c=3e8; % Velocity of Light (m/s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency Range
f=linspace(55,70,16001); % Frequncy in THz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:length(P)
for loop=1:length(f)
e1 = e01 - ((e01^2)/2) * (((p11_1/E1) * (v1+1) * P(k)) + ((p12_1/E1) * (3*v1+1) * P(k)));
nA = sqrt(e1);
% Formula for calculating nA
e2 = e02 - ((e02^2)/2) * (((p11_2/E2) * (v2+1) * P(k)) + ((p12_2/E2) * (3*v2+1) * P(k)));
nB = sqrt(e2);
%%% angular frequency
w=2*pi*f*1e12;
a=dA+dB; %%% Lattice Constant
DA=((w(loop))/c)*dA*nA; %%% wave number of first layer
DB=((w(loop))/c)*dB*nB; %%% wave number of second layer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first layer
m11=cos(DA); m12=-1i*sin(DA)/nA; m21=-1i*nA*sin(DA); m22=cos(DA);
mA=[m11 m12; m21 m22];
%%% Transfer MAtrix elements of Second layer
l11=cos(DB); l12=-1i*sin(DB)/nB; l21=-1i*nB*sin(DB); l22=cos(DB);
mB=[l11 l12; l21 l22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Unit cell
m=(mA*mB);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Photonic Band Structure
K(loop,k)=(1/(2*pi))*acos((m(1,1)+m(2,2))/2);
end
end
figure(1)
plot(f,real(K))
legend('P=0','P=50','P=100','P=150','P=200')
xlabel("Frequency (THz)")
ylabel("Wave vector")
