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;