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;   
        col = col + 1;   
           
    end   
       
    Angle(row) = ang;   
    row = row + 1;   
end   
   
   
   
   
figure(1); clf;   
   
plot (Wavelength,Reflectance);   
title('Reflectivity vs Wavelength');   
xlabel('Wavelength (nm)');   
ylabel('Reflectivity');   
hold on;   
grid on;   
zoom;   
   
   
figure(2); clf;   
   
plot (Wavelength,Transmittance);   
title('Transmittance vs Wavelength');   
xlabel('Wavelength (nm)');   
ylabel('Transmittance');   
hold on;   
grid on;   
zoom;