%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SolveHybridModeEVE.m
%
% V - single valued normalized frequency
% n - principle mode number
% l - secondary mode number
% nco - core refractive index
% ncl - cladding refractive index
%
% Solve the hybrid mode eigenvalue equation for a fiber with
% general core/cladding (i.e. not n_c~ n_cl)
% The working here is from Okamoto’s "Fundamenetals of optical
% waveguides" textbook, chapter 3
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function U = SolveHybridModeEVE(V,n,l,nco,ncl)

% Flags
% If you want to plot, turn this flag on
DoPlot = 0;
ShowWarning = 0;

% Define a finite range over which to look for solutions
Us = linspace(0,max(V),20000);

% Define w parameter
Ws = sqrt(V.^2-Us.^2);

% LHS of the eigenvalue equation
LHS1 = (1/2)*(besselj(n-1,Us)-besselj(n+1,Us)) ./ (Us.*besselj(n,Us)) - ...
(1/2)*(besselk(n-1,Ws) + besselk(n+1,Ws)) ./ (Ws.*besselk(n,Ws));
LHS2 = (1/2)*(besselj(n-1,Us)-besselj(n+1,Us)) ./ (Us.*besselj(n,Us)) - ...
(ncl/nco)^2*(1/2)*(besselk(n-1,Ws) + ...
besselk(n+1,Ws)) ./ (Ws.*besselk(n,Ws));
LHS = LHS1.*LHS2;

% RHS of EVE
RHS = n^2*(1./(Us.^2) + 1./(Ws.^2)).*(1./(Us.^2) + (ncl/nco)^2*(1./(Ws.^2)));

% Solve
[indUs,zeroUs] = crossing(LHS-RHS,Us); % Zero finding algorithm


% Sort and extract the lth solution to give the
% (n,l)th mode
if (length(indUs)>=l)
    imin = sort(indUs);
    soli=imin(l);
    U = Us(soli);
else
    % Check solution, issue a warning if nothing found
    soli = -1;
    if ShowWarning
        disp('WARNING: no solution for this m in chosen range')
    end
    U = 0;
end

% For debugging purposes, it an be useful to plot
currentfig= gcf;
if DoPlot
   figure(99)
   hold on
   plot(Us,LHS-RHS);
   line([U,U],[-10,10],'Color','k')
   axis([0,10,-10,10])
   xlabel('U')
   ylabel('LHS-RHS')
   hold off
end
figure(currentfig);
