%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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(58,67,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/(pi))*acos((m(1,1)+m(2,2))/2);
end
end
figure()
x = real([-K K]); % data
xtl = [-1 -0.95 -0.9 0.9 0.95 1]; % x-tick labels
xl = [-1 -0.87; 0.87 1]; % x-limits
gap_size = 0.05;
% transform data to (0,1) and split into two parts based on xl
[xt1,xt2] = x_transform(x,xl,gap_size);
% use 5 colors:
c = colororder();
colororder(c(1:size(xt1,2),:))
% plot 10 lines:
plot(xt1,f)
hold on
plot(xt2,f)
ylabel("Frequency (THz)")
xlabel("Wave vector")
ax = gca();
% legend:
lgd = legend('P=0','P=50','P=100','P=150','P=200');
title(lgd,'Hydrostatic Pressure')
lgd.EdgeColor = 'none';
lgd.FontWeight = 'bold';
lgd.Position(1) = ax.Position(1)+(ax.Position(3)-lgd.Position(3))/2;
lgd.Position(2) = ax.Position(2)+ax.Position(4)/2-lgd.Position(4);
% axes appearance:
[xtlt1,xtlt2] = x_transform(xtl,xl,gap_size);
set(ax, ...
'FontWeight','bold', ...
'Box','on', ...
'XLim',[0 1], ...
'XTick',[xtlt1 xtlt2], ...
'XTickLabel',xtl);
% gaps and gap lines:
yl = get(ax,'YLim');
dyl = yl(2)-yl(1);
set(ax,'YLimMode','manual');
patch( ...
'Parent',ax, ...
'XData',0.5+gap_size/2*[-1 -1; 0 0; 1 1; 0 0], ...
'YData',yl+dyl*[1 -1].*[-1; 1; 1; -1]/30, ...
'FaceColor','w', ...
'EdgeColor','none', ...
'Clipping','off', ...
'HandleVisibility','off');
line( ...
'Parent',ax, ...
'XData',reshape(0.5+gap_size/2*[1 -1].*[-1; 0; NaN; 1; 0; NaN],1,[]), ...
'YData',reshape(yl+dyl*[1 -1].*[-1; 1; NaN; 1; -1; NaN]/30,1,[]), ...
'Color','k', ...
'LineWidth',1, ...
'Clipping','off', ...
'HandleVisibility','off');
function [x1,x2] = x_transform(x,xl,gap_size)
[mx,nx] = size(x);
x1 = NaN(mx,nx);
x2 = NaN(mx,nx);
idx1 = x <= xl(1,2);
idx2 = x >= xl(2,1);
w = (1-gap_size)/2;
x1(idx1) = (x(idx1)-xl(1,1))/(xl(1,2)-xl(1,1))*w;
x2(idx2) = (x(idx2)-xl(2,1))/(xl(2,2)-xl(2,1))*w+1-w;
x1(:,all(isnan(x1),1)) = [];
x2(:,all(isnan(x2),1)) = [];
end
