% Experiment of Finding Beta
% This prog is based on Eq. 2.54 of Keiser and Eq. (10.86) of Rao.
% It plots on u-w plane the locus of modes
close all;clear;clc;clf reset
unset = 0.1:0.1:10;wnset = 0.1:0.1:10;
%wnset = sqrt(V^2 - unset.^2);
n1 = 1.50;n2 = 1.47;v = 0;%n1 = 1.47028; n2 = 1.45709;
lamda = 1.55e-6;k1 = 2*pi*n1/lamda;k2 = 2*pi*n2/lamda;
[un,wn]=meshgrid(unset,wnset);
%V = sqrt(un.^2 + wn.^2);
% Case of HE and EH Modes
% mode_eq = (besselj(v+1,un)./(un.*besselj(v,un))   +    besselk(v+1,wn)./(wn.*besselk(v,wn))).* ....
%      (n1^2*besselj(v-1,un)./(un.*besselj(v,un)) - n2^2*besselk(v-1,wn)./(wn.*besselk(v,wn)))+ .... 
%           (besselj(v-1,un)./(un.*besselj(v,un))   -    besselk(v-1,wn)./(wn.*besselk(v,wn))).* ....
%      (n1^2*besselj(v+1,un)./(un.*besselj(v,un)) + n2^2*besselk(v+1,wn)./(wn.*besselk(v,wn))); 
% mode_eq = (besselj(v+1,un)./un   +    besselj(v,un).*besselk(v+1,wn)./(wn.*besselk(v,wn))).* ....
%      (n1^2*besselj(v-1,un)./un - n2^2*besselj(v,un).*besselk(v-1,wn)./(wn.*besselk(v,wn)))+ .... 
%           (besselj(v-1,un)./un   -    besselj(v,un).*besselk(v-1,wn)./(wn.*besselk(v,wn))).* ....
%      (n1^2*besselj(v+1,un)./un + n2^2*besselj(v,un).*besselk(v+1,wn)./(wn.*besselk(v,wn)));  
% Case of TE Modes
% mode_eq = besselj(1,un)./(un.*besselj(0,un)) + besselk(1,wn)./(wn.*besselk(0,wn));% Actual equation, singular
mode_eq = besselj(1,un)./un      + besselj(0,un).*besselk(1,wn)./(wn.*besselk(0,wn));
% Case of TM Modes
% First two are actual equations, but they create singularities 
% mode_eq = (2*pi*n1/lamda)^2*besselj(1,un)./(un*besselj(0,un)) + (2*pi*n2/lamda)^2*besselk(1,wn)./(wn*besselk(0,wn)); 
% mode_eq = n1^2*besselj(1,un)./(un.*besselj(0,un)) + n2^2*besselk(1,wn)./(wn.*besselk(0,wn));  
%mode_eq = n1^2*besselj(1,un)./un + n2^2*besselj(0,un).*besselk(1,wn)./(wn.*besselk(0,wn)); 
mode_eqn = mode_eq / max(max(abs(mode_eq)));
figure(1)
set(gcf,'Renderer','Zbuffer');set(gcf,'Color',[1 1 1]);
meshc(un,wn,mode_eqn);set(gca,'FontSize',13);colormap colorcube;%view([-21.5 60.0])
title('3D Plot of Characteristic Equation for HE and EH Modes','FontSize',12,'FontWeight','bold');
figure(2)
[C,h] = contour(un,wn,mode_eq,'--k');axis square;grid on
set(h,'LineWidth',2);H = clabel(C,h,[0]);set(H,'FontSize',10,'LineWidth',2);
set(gca,'FontSize',13);set(gcf,'Color',[1 1 1]);hold on
V2 = un.^2 + wn.^2;V = sqrt(V2);[C,h] = contour(un,wn,V,'-k');set(h,'LineWidth',2);
H = clabel(C,h);set(H,'FontSize',12,'LineWidth',3);set(H,'LineWidth',2,'BackgroundColor','white','Edgecolor',[.1 .8 .1]);
title('Plot of \itu_n \rm\bf and \itw_n against different V circles','FontSize',12,'FontWeight','bold');
xlabel('\itu_n axis','FontSize',12,'FontWeight','bold');
ylabel('\itw_n axis','FontSize',12,'FontWeight','bold');