I=sqrt(-1);
Mode_Type = menu('Which type of mode do you want to study?','TE Mode','TM Mode' );
N=input('Number of pairs of layers : '); %Define the number of pairs of layers for the stack
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Structure_Parameter=zeros(4*N+1,3);
Structure_Parameter([1:2:4*N+1],1)=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Incident_Index=input('Incident Index : '); %Define the refractive index of the incident medium
Index_1=input('Type 1 layer index : '); %Define the refractive index of the type 1 layer in the stack
Index_2=input('Type 2 layer index : '); %Define the refractive index of the type 2 layer in the stack
Substrate_Index=input('Substrate Index : '); %Define the refractive index of the substrate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Structure_Parameter(1,2)=Incident_Index;
Structure_Parameter(1,3)=Index_1;
Structure_Parameter([3:4:4*N-1],2)=Index_1;
Structure_Parameter([3:4:4*N-1],3)=Index_2;
Structure_Parameter([5:4:4*N-3],2)=Index_2;
Structure_Parameter([5:4:4*N-3],3)=Index_1;
Structure_Parameter(4*N+1,2)=Index_2;
Structure_Parameter(4*N+1,3)=Substrate_Index;
Structure_Parameter([2:4:4*N-2],2)=Index_1;
Structure_Parameter([4:4:4*N],2)=Index_2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thickness_1=input('Type 1 layer thickness (nm) : '); %Define the thickness of the type 1 layer in the stack
thickness_1=thickness_1*1e-9;
thickness_2=input('Type 2 layer thickness (nm) : '); %Define the thickness of the type 2 layer in the stack
thickness_2=thickness_2*1e-9;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Structure_Parameter([2:4:4*N-2],3)=thickness_1;
Structure_Parameter([4:4:4*N],3)=thickness_2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lowest_Incident_Angle=input('Lowest Incident Angle (deg) : '); %Define lowest incident angle to be calculated
Highest_Incident_Angle=input('Highest Incident Angle (deg) : '); %Define highest incident angle to be calculated
Angle_step=input('Angle step (deg) : '); %Define angle change step
Lowest_Wavelength=input('Lowest wavelength (nm) : '); %Define the lowest wavelenght to be calculated
Lowest_Wavelength=Lowest_Wavelength*1e-9;
Highest_Wavelength=input('Highest wavelength (nm) : '); %Define the highest wavelenght to be calculated
Highest_Wavelength=Highest_Wavelength*1e-9;
Wavelength_step=input('Wavelength step (nm) : '); %Define the wavelenght change step
Wavelength_step=Wavelength_step*1e-9;
x_layer = menu('Do you want to include an x type layer?','Yes','No' ); %Ask the user to define if he wishes to include an x type layer
while x_layer == 1
Index_x=input('Type x layer index : '); %Define the refractive index of the type x layer in the stack
thickness_x=input('Type x layer thickness (nm) : '); %Define the thickness of the type x layer in the stack
thickness_x = thickness_x*1E-9;
m = input('Number of pair of layers that the replacement take place : ');
replace = menu('Which type of layer does it replace?' , 'Type 1' , 'Type 2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Structure_Parameter(4*m+2*replace-5,3)=Index_x;
Structure_Parameter(4*m+2*replace-4,2)=Index_x;
Structure_Parameter(4*m+2*replace-4,3)=thickness_x;
Structure_Parameter(4*m+2*replace-3,2)=Index_x;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_layer = menu('Do you want to include an x type layer?','Yes','No' ); %Ask the user to define if he wishes to include an x type layer
end
%Converting the range of incident angles to radians
Lowest_Incident_Angle = Lowest_Incident_Angle * (pi / 180);
Highest_Incident_Angle = Highest_Incident_Angle * (pi / 180);
Angle_step = Angle_step * (pi / 180 );
%Initialize array index to be loop
Reflectance = [];
Transmittance = [];
Angle = [];
Wavelength = [];
PhaseChange = [];
%Begin outer loop for calculating incident angle for the range specified before
row = 1; %This is the angle
for ang = Lowest_Incident_Angle : Angle_step : Highest_Incident_Angle
k_s=asin(Incident_Index*sin(ang)/Substrate_Index);
%Working out the optical admittances of the layers
Y = 2.6544*1e-3;
Incidence_admittance = (Incident_Index)*Y*cos(ang);
Substrate_admittance = (Substrate_Index)*Y*cos(k_s);
%Begin inner loop for calculating the range of wavelengths specified
col = 1; %This is the wave length number
for lambda = Lowest_Wavelength : Wavelength_step : Highest_Wavelength
M=[1 0;0 1];
if Mode_Type==1
for i=1:4*N+1
if Structure_Parameter(i,1)==1
n1=Structure_Parameter(i,2);
n2=Structure_Parameter(i,3);
k1=asin(Incident_Index*sin(ang)/n1);
k2=asin(Incident_Index*sin(ang)/n2);
A=n1*cos(k1)/(2*n2*cos(k2));
Q=[0.5+A 0.5-A;0.5-A 0.5+A];
elseif Structure_Parameter(i,1)==0
k=asin(Incident_Index*sin(ang)/Structure_Parameter(i,2));
fi = ((2*pi*Structure_Parameter(i,2))/(lambda))*Structure_Parameter(i,3)*cos(k);
Q=[exp(I*fi) 0 ; 0 exp(-I*fi)];
end
M=Q*M;
end
elseif Mode_Type==2
for i=1:4*N+1
if Structure_Parameter(i,1)==1
n1=Structure_Parameter(i,2);
n2=Structure_Parameter(i,3);
k1=asin(Incident_Index*sin(ang)/n1);
k2=asin(Incident_Index*sin(ang)/n2);
A=n1*cos(k1)/(2*n2*cos(k2));
B=cos(k1-k2);
C=cos(k1+k2);
Q=[B*(A+0.5) C*(A-0.5);C*(A-0.5) B*(A+0.5)];
elseif Structure_Parameter(i,1)==0
k=asin(Incident_Index*sin(ang)/Structure_Parameter(i,2));
fi = ((2*pi*Structure_Parameter(i,2))/(lambda))*Structure_Parameter(i,3)*cos(k);
Q=[exp(I*fi) 0 ; 0 exp(-I*fi)];
end
M=Q*M;
end
end
r=-M(2,1)/M(2,2); %Define reflection coefficient
t=M(1,1)+M(1,2)*r; %Define transmission coefficient
R1 = abs(r).^2; %Define reflectance
R2 = R1;
k_s = asin(Incident_Index*sin(ang)/Substrate_Index);
T1 = ((Substrate_Index*cos(k_s))*(abs(t).^2))/(Incident_Index*cos(ang)); %Define Transmittance
T2 = T1;
Reflectance(col,row) = R2;
Transmittance(col,row) = T2;
%Calculating elements of the characteristic matrix
B = M(1,1) + M(1,2)*Substrate_admittance;
C = M(2,1) + M(2,2)*Substrate_admittance;
%Define and calculate phase change
PC = abs((Incidence_admittance*B - C)/(Incidence_admittance*B + C));
T3 = PC;
PhaseChange(col,row) = T3;
Wavelength(col) = lambda;