%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is used to 
%   (1) Calculate the index of refraction of various glass materials 
%       that could be used in optical fiber
%   (2) Calculate the core radius for a particular cut-off wavelength
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Part 1 The glass materials
%
% Use the glas materials in the reference
%   J. W. Fleming, "Material dispersion in lightguide glasses,"
%   Electr. Lett., vol. 14, p. 326-328, May 1978.
%
% The indices of refraction are described using the Sellmeier equations
% given by
%       n^2 - 1 = sum( (Ai * lambda^2)/( lambda^2 - li))
%       The wavelength, lambda, is in units of micrometers
%
% These are the coefficients for the 6 materials
A1 = [ .696750; .711040; .695790; .690618; .691116; .796468];
A2 = [ .408218; .451885; .452497; .401996; .399166; .497614];
A3 = [ .890815; .704048; .712513; .898817; .890423; .358924];
L1 = [ .069066; .064270; .061568; .061900; .068227; .094359];
L2 = [ .115662; .129408; .119921; .123662; .116460; .093386];
L3 = [9.900559;9.425478;8.656641;9.098960;9.993707;5.999652];

% basic parameters
c=3e8;

% plot of the varous indices of refraction
lambda=linspace(1.2, 1.6, 201); % the wavelength in micrometers

NN=[];
NG=[];
for lp_mat=1:6
    a1=A1(lp_mat);
    a2=A2(lp_mat);
    a3=A3(lp_mat);
    l1=L1(lp_mat);
    l2=L2(lp_mat);
    l3=L3(lp_mat);
    
    %calculate the refractive indices
    n = sqrt(1+(a1.*lambda.^2)./(lambda.^2-l1.^2)+(a2.*lambda.^2)./(lambda.^2-l2.^2)+(a3.*lambda.^2)./(lambda.^2-l3.^2));

    %make a matrix of refractive indices
    NN=[NN;n];


    %calculate the group refractive indices
    w=2*pi*c./lambda;
    dw=diff(w);
    ng=n(2:end)+w(2:end).*diff(n)./diff(w);
   
    %make a matrix of group refractive indices
    NG=[NG;ng];

end

%the indices of refraction as a function of wavelength
subplot(2,1,1)
plot(lambda,NN)

subplot(2,1,2)
plot(lambda(2:end),NG)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Part 2 Cut-off
%
%   Calculate the waveguide radius for a a particular cut-off wavelength
%   Use the equation 2.405 = 2*pi/lambda*a*sqrt(ncore^2-nclad^2)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_cutoff=1.2;   %cut-off wavelength
core=2; %number for the core material
clad=1; %number for the cladding material

%check to see if these materials are valid
I=find(NN(core,:)<NN(clad,:));
if length(I)>0
    'NOT VALID MATERIALS'
else
    
    %calculate the core radius for the particular cut-off wavelength
    n_core = sqrt(1+(A1(core).*l_cutoff^2)./(l_cutoff^2-L1(core).^2)+(A2(core).*l_cutoff^2)./(l_cutoff^2-L2(core).^2)+(A3(core).*l_cutoff^2)./(l_cutoff^2-L3(core).^2))
    n_clad = sqrt(1+(A1(clad).*l_cutoff^2)./(l_cutoff^2-L1(clad).^2)+(A2(clad).*l_cutoff^2)./(l_cutoff^2-L2(clad).^2)+(A3(clad).*l_cutoff^2)./(l_cutoff^2-L3(clad).^2))

    a=2.405*l_cutoff/(2*pi*sqrt(n_core^2-n_clad^2))
        
end