% Test hybrid mode EVE solver by plotting beta vs V dispersion 
% curves

% Clean up workspace
clear all
close all

% Set parameters
nlist = 1:4; % List of principle mode numbers
l = 1; % Secondary mode number
nco = 1.45; % Core index
ncl = 1; % Cladding index
lambda = 800e-9; % Wavelength
k=2*pi/lambda; % Wave number
a = (linspace(200,2000,200)/2)*1e-9; % Here, I chose to vary V by changing the fiber radius. Changing the wavelength is also okay.
V = k*a*sqrt(nco^2-ncl^2); 

% Set up figure
figure(1)
clf

for kk = 1:length(nlist)

    for ll = 1:length(V)
        U(ll) = SolveHybridModeEVE(V(ll),nlist(kk),l,nco,ncl); 
    end
    
    % Calculate beta
    betas = sqrt(k^2*nco^2 - (U./a).^2);
    betas(find(U==0)) = k*ncl; % Cases where no solution is found correspond to free space beta
    
    % Calculate other quantities of interest from beta
    bs = ((betas/k).^2 - ncl^2) / (nco^2-ncl^2);
    neffs = betas/k;
    figure(1)
    hold on
    plot(V,bs,'b-')
    hold off


end

% Format plot
set(gca,'FontSize',14)
xlabel('V','FontSize',16)
ylabel('\beta','FontSize',16)
axis tight