function [neff, Lambda, n_core, n_clad] = fiber_mode_calc(a, core, clad)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is used to calculate the mode of a step index optical fiber.
% 
% Inputs 
%   a: core radius
%   core: number of the material used for the core
%   clad: number of the material used for the cladding
%
%Outputs
%   neff: vector of effective mode indices
%   Lambda: vector of wavelength (in m)
%   n_core: vector of core mode indices
%   n_clad: vector of cladding indices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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];

% Fiber parameters
%a=1.8e-6;   %core radius in meters
%core=2; %number for the core material
%clad=1; %number for the cladding material

c=3e8;
lambda=linspace(1.2, 1.6, 201); % the wavelength in micrometers

%check to see if these materials are valid
n_core = sqrt(1+(A1(core).*lambda.^2)./(lambda.^2-L1(core).^2)+(A2(core).*lambda.^2)./(lambda.^2-L2(core).^2)+(A3(core).*lambda.^2)./(lambda.^2-L3(core).^2));
n_clad = sqrt(1+(A1(clad).*lambda.^2)./(lambda.^2-L1(clad).^2)+(A2(clad).*lambda.^2)./(lambda.^2-L2(clad).^2)+(A3(clad).*lambda.^2)./(lambda.^2-L3(clad).^2));
I=find(n_core<n_clad);

if length(I)>0
    'NOT VALID MATERIALS'
else

    %use SI units now
    Lambda=lambda*1e-6;

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Part 2 Calculate the effective index of refraction
%
%   (a) Loop through each of the wavelengths in the vector Lambda
%   (b) Solve the equation J1(pa)/J0(pa) = q/a*K1(qa)/K0(qa)
%       This calculation needs to be very exact in order to numerical
%       derivatives
%   (c) Save the calculated neff to a vector

clear Bz FMIN_Error b k0_mat neff

    for lp_lam = 1:length(Lambda)   %loop through wavelengths
        [lp_lam, length(Lambda)]
        ll=lambda(lp_lam);
        ncore = n_core(lp_lam);
        nclad = n_clad(lp_lam);
    
        k0 = 2*pi/Lambda(lp_lam);
        B1 =k0.*nclad;
        B2 =k0.*ncore;

        %the range of possible propagation constants
        B=linspace(B1,B2,500);
        B=B(2:end);

        p=((ncore*k0)^2-B.^2).^.5;
        q=(B.^2-(nclad*k0)^2).^.5;
  
    %the two sides of the mode equation
        LHS = besselj(1,p*a)./(besselj(0,p*a));
        RHS = (q./p).*besselk(1,q*a)./besselk(0,q*a);
        
        J1=find(isinf(LHS));
        J2=find(isinf(RHS));
        
        LHS(J1)=1e10*ones(size(J1));
        RHS(J2)=1e10*ones(size(J2));
        
%find the beta where the two sides are the closest together
        Diff1 = 1e20*abs(LHS-RHS);
        [C,I]= min(Diff1);

%The calculation of beta needs to be very exact otherwise the numerical 
%derivates will not give the correct answer
        lp_no=0;
        ret=0;
        while C>1e-10 & ret==0
            if lp_no>30
                ret=1;
            end
            lp_no=lp_no+1;
%        [lp_a, lp_no]

            if I==1
                BI1=B2;
                BI2=B(I+1);
            else
                BI1=B(I-1);
                BI2=B(I+1);
            end
    
            B = linspace(BI1,BI2,100);

            p=((ncore*k0)^2-B.^2).^.5;
            q=(B.^2-(nclad*k0)^2).^.5;
  
    %the two sides of the mode equation
            LHS = besselj(1,p*a)./(besselj(0,p*a));
            RHS = (q./p).*besselk(1,q*a)./besselk(0,q*a);
        
            J1=find(isinf(LHS));
            J2=find(isinf(RHS));
        
            LHS(J1)=1e10*ones(size(J1));
            RHS(J2)=1e10*ones(size(J2));
        
%find the beta where the two sides are the closest together
            Diff1 = 1e20*abs(LHS-RHS);
            [C,I]= min(Diff1);
        end
        
%Save the calculated mode        

        Btest = B(I);   
        Bz(lp_lam)=Btest;
    
        FMIN_Error(lp_lam) = C; %the error
    
        b(lp_lam)=((Btest/k0)^2-nclad^2)/(ncore^2-nclad^2);
        k0_mat(lp_lam) = k0;
        neff(lp_lam) = Btest/k0;
    end
end