hello
I simply added a for loop for the del parameter to see how peak amplitude and peak width of the red curve would evolve vs del (smooth or non smooth)
before that , I made a few mods to speed up the code
- reduced the freq range to what you finally display (instead of computing lot more results that are then throw away)
- modified how you do the peak width computation (here my suggestion will do linear interpolation, so we can use a slightly coarser freq vector without loss of accuracy)
if we use the new code with the complete range for del (from -1 to +1) , we got his result :
It's not evident to find a del value that correspond to a minimum point while having the transmission factor as close as possible to one. In fact it seems there are many possible values that correspond to this request
after some time , we can see that there is an interesting range between 0 and 0.2 that we can further investigate
narrowing down to the del 0 to 0.2 range gives this result
what would be the best del value here ?
seems to me del between 0.04 and 0.08 is a good value
now if you want to selct ONE optimum del value you have to define a criteria that combine peak amplitude and peak width with the correct weighting and minimize this value
I don't know right now how you would define the ponderation here - it's up to you
code :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a program for 1-D Photonic crystal
% For Photonic band structure calculation and Transmission Spectrum
% This program uses Transfer Matrix Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters
P = 0; %%% 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
e1 = e01 - ((e01^2)/2) * (((p11_1/E1) * (v1+1) * P) + ((p12_1/E1) * (3*v1+1) * P));
nA = sqrt(e1);
% Formula for calculating nA
e2 = e02 - ((e02^2)/2) * (((p11_2/E2) * (v2+1) * P) + ((p12_2/E2) * (3*v2+1) * P));
nB = sqrt(e2);
% del_range= -1:0.01:1; %%% Thickness variation parameter (full range search)
del_range= 0:0.001:0.2; %%% Thickness variation parameter (narrow range search)
% Initialization of parameters
dair=1200e-9; % Thickness of air in meters
dA=780e-9; % Thickness of First Layer in meters
dB=830e-9; % Thickness of Second Layer in meters
% ddA=(1+del)*dA; % Thickness of First Layer with thickness variation in meters
% ddB=(1+del)*dB; % Thickness of Second Layer with thickness variation in meters
c=3e8; % Velocity of Light (m/s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency Range
% f=linspace(58,67,16001); % Frequncy in THz
% f=58:0.001:67; % Frequncy in THz
f=60.5:0.01:61.5; % Frequncy in THz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:numel(del_range)
del = del_range(k);
ddA=(1+del)*dA; % Thickness of First Layer with thickness variation in meters
ddB=(1+del)*dB; % Thickness of Second Layer with thickness variation in meters
for loop=1:numel(f)
%%% 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
Da=((w(loop))/c)*(dA*0.5)*nA; %%% wave number of first half layer
Db=((w(loop))/c)*(dB*0.5)*nB; %%% wave number of second half layer
DDA=((w(loop))/c)*ddA*nA; %%% wave number of first layer with thickness variation
DDB=((w(loop))/c)*ddB*nB; %%% wave number of second layer with thickness variation
DDa=((w(loop))/c)*(ddA*0.5)*nA; %%% wave number of first half layer with thickness variation
DDb=((w(loop))/c)*(ddB*0.5)*nB; %%% wave number of second half layer with thickness variation
Dair=((w(loop))/c)*(dair*0.5)*n0; %%% wave number of air layer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of air
a11=cos(Dair); a12=-1i*sin(Dair)/n0; a21=-1i*n0*sin(Dair); a22=cos(Dair);
AA=[a11 a12; a21 a22];
%%% 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];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first half layer
ma11=cos(Da); ma12=-1i*sin(Da)/nA; ma21=-1i*nA*sin(Da); ma22=cos(Da);
ma=[ma11 ma12; ma21 ma22];
%%% Transfer MAtrix elements of Second half layer
lb11=cos(Db); lb12=-1i*sin(Db)/nB; lb21=-1i*nB*sin(Db); lb22=cos(Db);
mb=[lb11 lb12; lb21 lb22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first layer with thickness variation
dm11=cos(DDA); dm12=-1i*sin(DDA)/nA; dm21=-1i*nA*sin(DDA); dm22=cos(DDA);
mA1=[dm11 dm12; dm21 dm22];
%%% Transfer MAtrix elements of Second layer with thickness variation
dl11=cos(DDB); dl12=-1i*sin(DDB)/nB; dl21=-1i*nB*sin(DDB); dl22=cos(DDB);
mB2=[dl11 dl12; dl21 dl22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Transfer Matrix elements of first half layer with thickness variation
dma11=cos(DDa); dma12=-1i*sin(DDa)/nA; dma21=-1i*nA*sin(DDa); dma22=cos(DDa);
ma1=[dma11 dma12; dma21 dma22];
%%% Transfer MAtrix elements of Second half layer with thickness variation
dlb11=cos(DDb); dlb12=-1i*sin(DDb)/nB; dlb21=-1i*nB*sin(DDb); dlb22=cos(DDb);
mb2=[dlb11 dlb12; dlb21 dlb22];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unit cell
m=(mA*mB);
M0=(m^25); %%% A-B structure
M1=(mb*mA*mb)^25; %%% B/2-A-B/2 structure
M2=(ma*mB*ma)^25; %%% A/2-B-A/2 structure
M3=(mb2*mA1*mb2)^25; %%% B/2-A-B/2 structure with thickness variation
M4=(ma1*mB2*ma1)^25; %%% A/2-B-A/2 structure with thickness variation
M = M1*M2;
% M = M2*M1;
% M = M0*AA*M0; %%% defect structure (air)
Mm = M1*M2*M3;
% M = M2*M1*M4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Reflection and Transmission
% R1(loop)=((ns*M(1,1)+ns*n0*M(1,2)-M(2,1)-n0*M(2,2))/(ns*M(1,1)+ns*n0*M(1,2)+M(2,1)+n0*M(2,2)));
% R(loop)=(R1(loop))*conj(R1(loop));
T1(loop)=(2*ns/(ns*M(1,1)+ns*n0*M(1,2)+M(2,1)+n0*M(2,2)));
T(loop)=(n0 / ns) * abs(T1(loop))^2;
Tt1(loop)=(2*ns/(ns*Mm(1,1)+ns*n0*Mm(1,2)+Mm(2,1)+n0*Mm(2,2)));
Tt(loop)=(n0 / ns) * abs(Tt1(loop))^2;
end
% red line data processing
ft = 0.018+f;
[peak_val(k),~] = max(Tt);
half_max = peak_val(k) / 2;
[f_left_edge,f_right_edge] = find_zc(f,T,half_max);
peak_width(k) = f_right_edge - f_left_edge;
end
plot(del_range,peak_val,'b',del_range,peak_width,'r');
xlabel("Del")
legend('peak val','peak width')
return
figure(1)
plot(f,T(:),'k') %%% Topological structure
hold on
plot(ft,Tt(:),'r') %%% Topological coupled structure
hold off
% xlim([60.5,61.5])
% xlim([58,67])
ylim([0,1])
xlabel("Frequency (THz)")
ylabel("Transmission")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)
% put data in rows
x = x(:);
y = y(:);
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
