%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code calculates the modes of a symmetric slab waveguide
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% variables:
% lr: wavelength
% n1: cladding index
% n2: core index
% lr: wavlength (in micrometers)
% d: vector with film thicknesses
% pol: polarization state (1: TE, 2:TM)
%
% All length units are in micrometers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global V
%waveguide parameters
%required parameters: n1, n2, d, lr
n1=1.5; %top region (the cover)
n2=1.6; %waveguide film indices
d=6; %waveguide film thicknesses
lr=1.0; %wavelength
ko=2*pi/lr;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Equations:
%
%Dispersion equations
% h=sqrt((n2*ko)^2-beta^2)
% q=sqrt(beta^2-(n1*ko)^2)
%combine these 2 equations
% h^2+q^2 = (n2^2-n1^2)*ko^2
%
%Mode equations
% h*tan(h*d/2) = q (symmetric TE modes)
% -h*cot(h*d/2)= q (antisymmetric TE modes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multiply both sides of the equations by d
% call x=hd, y=qd
%
% h*tan(h*d/2) = q becomes x*tan(x/2)=y2
% -h*cot(h*d/2) = q becomes -x*cot(x/2)=y3
%
% h^2+q^2 = (n2^2-n1^2)*ko^2 becomes x^2+y1^2 = d^2*ko^2*(n2^2-n1^2) = V
%
%normalized frequency
V=d*pi/lr*sqrt(n2^2-n1^2)
%break x up into 2 pi intervals
intervals=ceil(1.01*V)
%(1) plot the three equations
% Dispersion equation
% Symmetric modes
% Antisymmetric modes
figure(1)
clf
hold 'on'
if intervals==1
x=linspace(0, 1.5*V,101);
x=sort([x,V]);
I=find(x<=V);
y1=sqrt(V^2-x.^2); %dispersion
y2=x.*tan(x/2); %symmetric
y3=-x.*cot(x/2); %antisymetric
plot(x(I),y1(I),'b')
plot(x,y2,'r')
plot(x,y3,'g')
grid
axis image
axis([0, max(x), 0 max(x)])
xlabel('h d')
ylabel('q d')
hold 'off'
else
for lp=1:intervals
%only plot in pi/2 intervals to eliminate positive/negative changes
x=linspace(pi/2*(lp-1),pi/2*lp,101);
x=x(2:end);
I=find(x0
for lp_asym=1:n_asym
tmp=fminbnd('asym_modes',pi/2+(lp_asym-1)*pi,pi+(lp_asym-1)*pi)
x_asym=[x_asym,tmp]
end
end
hd=2*sort([x_sym,x_asym]);
beta=sqrt((n2*ko)^2-(hd/d).^2)
neff=beta/ko
%(4) plot the mode fields
figure(2)
x1=linspace(-d/2,d/2,101);
x2=linspace(d/2,2*d,101);
x3=linspace(-2*d,-d/2,101);
nn=0;
while nn