%計算纖芯導模的傳播常數與有效折射率
%   equation:T*a1*J(1,T*a1)/J(0,T*a1)=tao*a1*K(1,tao*a1) /K(0,tao*a1)
%
%    a1:radius of the core,a2:radius of the cladding  (um)
%    T:transverse phase shift constant,tao:transverse attenuate constant
%    n1:index of the core, n2:index of the cladding
%    lamda:wavelength in vacuum   (um)
%    V:normalized frequency
%    beta:longitudinal propagation constant
%    k:  wave vector
%    n_eff: effective index of the core
%    x=T*a1, y=tao*a1
format long
a1=4.15;     a2=62.5;
%a1=2.625;     a2=62.5;
% n1=1.458;    n2=1.450;
n1=1.458;    n2=1.45;
lamda=1.55;


k1=2*pi*n1/lamda;
V=2*pi*a1*(n1^2-n2^2)^0.5/lamda
step = 0.000001;
x=0:step:V;
y=(V^2-x.^2).^0.5;
n_co_eff = (k1^2-(x/a1).^2).^0.5*lamda/(2*pi);
% p3=x.*besselj(1,x)./besselj(0,x);
% p4=y.*besselk(1,y)./besselk(0,y);
p3=x.*besselj(1,x).*besselk(0,y);
p4=y.*besselk(1,y).*besselj(0,x);
figure(1);grid on;hold on;
plot(n_co_eff,p3)
plot(n_co_eff,p4,'r')
[mimi,i]=min(abs(p4-p3))
latent_root=i*step/a1;    %橫向傳播常數的特征值T
beta=(k1^2-latent_root^2)^0.5  %縱向傳播常數
co_eff=beta*lamda/(2*pi)