%waveguide with exponential distribution
%This script calculates the analytic values of effective indexes and
%compares them to simulation values

%By default, this script analyzes the first 2 TE modes and 1st 2 TM modes.
%This can be changed in the figure section

%The analytic expressions are found in "Modes in optical waveguides formed
%by diffusion" by E M Conwell

%(There is a slight chance of minor error with the tm expression - see
%below)
function [nTE,nTM]=exp_wg

lam = .633e-6;
d = 2.5e-6; %diffusion length
% e(y) = e0 + eDelta*exp(-|y|/d)
eDelta = 0.3;
e0 = 6.1;
e2 = 1; %air
k0 = 2*pi/lam;
% sqrt(e0 + eDelta) > b/k0 > sqrt(e0)
b0 = linspace( k0*sqrt(e0), k0*sqrt(e0 + eDelta), 1000);
b0 = b0(1:end-1);

te0=TE_eq(b0,eDelta,e0,e2,k0,d,lam);
tm0=TM_eq(b0,eDelta,e0,e2,k0,d,lam);

%TE
intervals=(te0>=0)-(te0<0);
izeros=find(diff(intervals)>0);
X0=[b0(izeros); b0(izeros+1)]';

[nzeros,scrap]=size(X0);

for i=1:nzeros
    nTE(i)=fzero(@(x) TE_eq(x,eDelta,e0,e2,k0,d,lam),X0(i,:))/k0;
end
nTE=nTE(end:-1:1);

%TM
intervals=(tm0>=0)-(tm0<0);
izeros=find(diff(intervals)>0);
X0=[b0(izeros); b0(izeros+1)]';

[nzeros,scrap]=size(X0);

for i=1:nzeros
    nTM(i)=fzero(@(x) TM_eq(x,eDelta,e0,e2,k0,d,lam),X0(i,:))/k0;
end
nTM=nTM(end:-1:1);

function te0=TE_eq(b0,eDelta,e0,e2,k0,d,lam)

zeta =4*pi*sqrt(eDelta)*d/lam;  %argument for bessel function
p0_0 = sqrt( b0.^2 - e0*k0^2 );
p2_0 = sqrt( b0.^2 - e2*k0^2 );

te0 = ( besselj(2*d*p0_0-1,zeta) - besselj(2*d*p0_0+1,zeta) ) ./ besselj(2*d*p0_0,zeta) + p2_0*lam/pi/sqrt(eDelta);

function tm0=TM_eq(b0,eDelta,e0,e2,k0,d,lam)

zeta =4*pi*sqrt(eDelta)*d/lam;  %argument for bessel function
p0_0 = sqrt( b0.^2 - e0*k0^2 );
p2_0 = sqrt( b0.^2 - e2*k0^2 );
%the term J' appears for tm and it is defined (in a 1974 paper) to be the
%derivative with respect to its argument.  By comparing the te equations in
%73 and 74, J' is assumed to be the difference of J divided by 2
%Therefore, ff and fff below are used to reflect this assumption.
ff = 1;   
fff = 2;

tm0 = ( besselj(2*d*p0_0-1,zeta) - besselj(2*d*p0_0+1,zeta) ) ./ besselj(2*d*p0_0,zeta) + fff*lam*sqrt(eDelta)/4/pi/e0/d + ...
        e0/e2*p2_0*lam/ff/pi/sqrt(eDelta) * (1+eDelta/e0);