Bragg_reflect

%Notes:
clear all
n_1=1;                              %Refractive index of air
n_2=2.3;                            %Refractve index of zinc sulphite
n_3=1.35;                           %Refractive index of cryolite
n_l=1.5;                            %Refractive index of glass
wav_0=500;                          %Bragg wavelength
wavs=[390:1:700];                   %Visible spectrum that will examined
N=[2 5 10 20];                      %Stacks that will be examined
R=zeros(length(N),length(wavs));    %Matrix where reflectivity will be stored


h_2=wav_0/(4*n_2);                  %Thickness for a layer of zinc sulphite
h_3=wav_0/(4*n_3);                  %Thickness for a layer of cyrolite
h=h_2+h_3;                          %Thickness of a stack

ang_n1=deg2rad(90);                 %Angle of incidence in degrees; 0 is normal incidence
ang_n2=asin((sin(ang_n1)*n_1)/n_2); %Angle after refraction through zinc sulphite
ang_n3=asin((sin(ang_n2)*n_2)/n_3); %Angle after refraction through cryolite
ang_nl=asin((sin(ang_n3)*n_3)/n_l); %Angle after refraction through glass

p_1=n_1*cos(ang_n1);
p_2=n_2*cos(ang_n2);
p_3=n_3*cos(ang_n3);
p_l=n_l*cos(ang_nl);

for j=1:4                           %For loop which tests with each no. of stacks
for k=1:311                         %For loop which finds reflectivity for a given wavelength and no. of stacks

    b_2=((2*pi)/wavs(k))*p_2*h_2;
    b_3=((2*pi)/wavs(k))*p_3*h_3;

    %a=cos(b_2)*cos(b_3)-0.5*((p_2/p_3)+(p_3/p_2))*sin(b_2)*sin(b_3)

%     u=zeros(1,N(j));
% for o=1:N(j);
%       u(1,o)=(sin((o-1)*acos(a)))/sqrt(1-a^2);
% end

%     b_2=((2*pi)/wavs(k))*p_2*h_2;
%     b_3=((2*pi)/wavs(k))*p_3*h_3;

    M_2=[cos(b_2) -i/p_2*sin(b_2);-i*p_2*sin(b_2) cos(b_2)];
    % ^Characteristic matrix for a single layer of zinc sulphite
    M_3=[cos(b_3) -i/p_3*sin(b_3);-i*p_3*sin(b_3) cos(b_3)];
    % ^Characteristic matrix for a single layer of cyrolite
    M_h=M_2*M_3;
    % ^Characteristic matrix for zinc sulphite and cyrolite layers together
    % (a "stack")
    M_2N=M_h^N(j);
    % ^Characteristic matrix for N pairs of zinc sulphite and cyrolite
    % layers, where N is the number of stacks

%     M_2N(1,1)=M_h(1,1)*u(1,N)-u(1,N-1);
%     M_2N(1,2)=M_h(1,2)*u(1,N);
%     M_2N(2,1)=M_h(2,1)*u(1,N);
%     M_2H(2,2)=M_h(2,2)*u(1,N)-u(1,N-1);

    r=((M_2N(1,1)+(M_2N(1,2)*p_l))*p_1-(M_2N(2,1)+(M_2N(2,2)*p_l)))/((M_2N(1,1)+(M_2N(1,2)*p_l))*p_1+(M_2N(2,1)+(M_2N(2,2)*p_l)));
    % ^Reflection coefficient for N stacks
    R(j,k)=abs(r^2);                %Reflectivity for the no. of stacks and wavelength
end
end

figure; hold on
plot(wavs,R(1,:),'r-')
plot(wavs,R(2,:),'b-')
plot(wavs,R(3,:),'g-')
plot(wavs,R(4,:),'y-')
xlim([390 700])
legend('No. of stacks = 2','No. of stacks = 5','No. of stacks = 10','No. of stacks = 20','location','northeast')
xlabel('Wavelength (nm)')
ylabel('Reflectance (dimensionless)')
title('Reflectance of a Bragg reflector over the visible spectrum (Incidence = \pi/2 rad)')

%poly=[1 2*a 4*a^2-1 8*a^3-4*a 16*a^4-12*a^2+1 32*a^5-32*a^3+6*a];
% u=zeros(1,2);
% for i=1:N_2;
%       u(1,i)=(sin(i*acos(a)))/sqrt(1-a^2);
% end