%Chap4Exercise5.m
%VCSEL mirror design
clear;
clf;

% '1' is GaAs layer, '2' is AlAs layer
wavelength=980e-9;  % center wavelength in meters
n1=3.5;             % 3.5 is refractive index of GaAs
n2=3.0;             % 3.0 is refractive index of AlAsd1=69.5e-9;         
d1=wavelength/n1/4 % thickness of dielectric with refractive n1 (GaAs layer)
d2=wavelength/n2/4 % thickness of dielectric with refractive n2 (AlAs layer)
d3=wavelength/n1   % thickness of dielectric with refractive n1 in center of pass-band filter

istep=0;
resolution=0.00005; % resolution in nm
range=[-100:1:-2 -1:resolution:1 1:1:100]*1e-9+wavelength; %wavelength scanning range

for wavelength=range
istep=istep+1;
k0=2*pi/wavelength; %input from air
k1=2*pi*n1/wavelength;
k2=2*pi*n2/wavelength;

% in EM, exp(-ikz) is often considered positive propagation +z direction
% Assuming light is propagating perpendicular to the layers moving
% from below (-z) to above (+z)
% A*exp(-ik1*z)+B*exp(ik1*z)=C*exp(-ik2*z)+D*exp(-ik2*z)
% D0,D1,D2 are interface matrices
% [D0]*[input from air]=[D1]*[output to n1],
% [D1]*[input from n1]=[D2]*[output to n2]

D0=[1 1
    k0 -k0];
D1=[1 1
    k1 -k1];
D2=[1 1
    k2 -k2];
    
% P1,P2,P3 are propagation matrices
% input from the below, output above material, [input from below]=[P]*[output above],
P1=[exp(i*k1*d1) 0
    0 exp(-i*k1*d1)];        
P2=[exp(i*k2*d2) 0
    0 exp(-i*k2*d2)];
P3=[exp(i*k1*d3) 0           % for dielectric with refractive index n1 in center of
    0 exp(-i*k1*d3)];        % passband filter
 
% Solution = 
% [Input A;B]=inv(D0)*{D1*P1*inv(D1)*D2*P2*inv(D2)}*{2nd period}*{3rd}
% *..*{defect,wider n1}*...*{D1*P1*inv(D1)*D2*P2*inv(D2)}*D0*[output C;D]
G1=D1*P1*inv(D1)*D2*P2*inv(D2);   	% for period of lower mirror
G2=D2*P2*inv(D2)*D1*P1*inv(D1);   	% for period of upper mirror
lower=inv(D0)*(G1^25)*D1;         	% propagation through lower 25 layer-pairs
upper=inv(D1)*(G2^25)*D0;         	% propagation through upper 25 layer-pairs
Total=lower*P3*upper;

r(istep)=Total(2,1)/Total(1,1);    	% reflected electromagnetic field
R(istep)=r(istep)*conj(r(istep));  	% total reflection intensity is r x r*

lower_r(istep)=lower(2,1)/lower(1,1);
lower_R(istep)=lower_r(istep)*conj(lower_r(istep));  % total reflection

upper_r(istep)=upper(2,1)/upper(1,1);
upper_R(istep)=upper_r(istep)*conj(upper_r(istep));  % total reflection

end    %%% loop to next value of istep

figure(1);
plot(range,lower_R);
title(['Chapt4Exercise5a, Bragg mirror']);ylabel('Reflection');xlabel('Wavelength, \lambda (m)');

%find FWHM of transmission
T=1-R;
N=length(T)-202

[tmp,centerindex]=max(T); % find index of peak center
Tmax=max(T);
i = 101;
while sign(T(i)-(Tmax/2)) == sign(T(i-1)-(Tmax/2))
    i = i+1;
end                              %first crossing is between (i-1) & (i)
interp = ((Tmax/2)-T(i-1)) / (T(i)-T(i-1));
llead = range(i-1) + interp*(range(i)-range(i-1));
i = centerindex+1;               %start search for next crossing at center
while ((sign(T(i)-(Tmax/2)) == sign(T(i-1)-(Tmax/2))) & (i <= N-1))
    i = i+1;
end

    interp = ((Tmax/2)-T(i-1)) / (T(i)-T(i-1));
    ltrail = range(i-1) + interp*(range(i)-range(i-1));
    FWHM = (ltrail - llead); %(m)

%FWHM=size(find(R<max(R)/2),2)*resolution*1e-9  %accurate only when one peak in range


figure(2);
plot(range,R);
title(['Chapt4Exercise5b, VCSEL mirror' ]);ylabel('Reflection');xlabel('Wavelength, \lambda (m)');
temp=[' FWHM = '  num2str(FWHM) '  m'];
text(1550e-9, 0.5, temp);

