光纤模式分布 matlab,matlab计算单模光纤模式分布(公布源代码及参考文献)
最近在使用matlab计算单模光纤纤芯模及包层模模场分布时，有一些问题一直悬而未决，多次咨询原作者后虽解决了部分问题，但是余下的问题原作者也不理我了，特发此贴以广交学习光纤方面的同学、老师及科研人员，希望大家互相学习，互相交流，共同探讨，共同进步。
首先，求解单模光纤纤芯模特征方程，光纤为Corning SMF-28，具体参数见以下代码，参考文献为：
% [1]. Optical Fibre Waveguide Analysis_Charles Tsao_Oxford 1992 ----% %
% Page 315--------------% %
% [2]. Eigenvalue and field equations of three-layered uniaxial fibers
% and their applications to the characteristics of long-period fiber
% gratings with applied axial strain ---------% %
% Zi-jia Zhang,J.Opt.Soc.Am. A/Vol.22,No,11/November 2005 --------%
% -------------------------------------------------------------------------
首先光纤模型采用的是三层模型，相应代码如下ecothf_HE.m：(用过matlab计算光纤模场分布的都应该能看懂下面的代码)
% -------------------------------------------------------------------------
% The Eigenvalue Equation of the Core Hybrid HE or EH mode for
% a general Three-layered Single mode Fiber ------------% %
% % --------------------------------Reference---------------------------% %
% [1]. Optical Fibre Waveguide Analysis_Charles Tsao_Oxford 1992 ----% %
% Page 315--------------% %
% [2]. Eigenvalue and field equations of three-layered uniaxial fibers
% and their applications to the characteristics of long-period fiber
% gratings with applied axial strain ---------% %
% Zi-jia Zhang,J.Opt.Soc.Am. A/Vol.22,No,11/November 2005 --------%
% -------------------------------------------------------------------------
function f = ecothf_HE(U1)
global lambda n1 n2 n3 ra1 ra2
ns1 = n1*n1;
ns2 = n2*n2;
ns3 = n3*n3;
s21 = ns2/ns1;
s23 = ns2/ns3;
kw = 2*pi/lambda; % kw = k0
kr1 = kw*ra1;
krs1 = kr1*kr1;
kr2 = kw*ra2;
krs2 = kr2*kr2;
bvs12 = krs1*(ns1-ns2); % V12^2 = k0^2*ra1^2*(n1^2-n2^2)
bvs23 = krs2*(ns2-ns3);
bus1 = U1*U1; % U1^2
nes = ns1-bus1/krs1; % neff^2 = n1^2-U1^2/(k0*ra1)^2
ne = sqrt(nes);
u1 = U1/ra1;
w2 = kw*sqrt(nes-ns2); % w2 = k0*sqrt(n2^2-neff^2)
w3 = kw*sqrt(nes-ns3); % w3 = k0*sqrt(neff^2-n3^2)
bu1 = u1*ra1; % U1 = u1*ra1
bw2 = w2*ra1; % W2 = w2*ra1
bws2 = bw2*bw2;
bw3 = w3*ra2; % W3 = w3*ra2
bws3 = bw3*bw3; % W3^2 = k0^2*ra2^2*(neff^2-n3^2)
u11 = bu1;
w21 = w2*ra1;
w22 = w2*ra2;
w32 = w3*ra2;
ar = ra2/ra1; % ar = alpha2
ars = ar*ar;
%%
j011 = besselj(0,u11); % j011 = besselj(0,u1*a1)
j111 = besselj(1,u11); % j111 = besselj(1,u1*a1)
j211 = besselj(2,u11); % j211 = besselj(2,u1*a1)
i021 = besseli(0,w21); % i021 = besseli(0,w2*a1)
i121 = besseli(1,w21); % i121 = besseli(1,w2*a1)
i221 = besseli(2,w21); % i221 = besseli(2,w2*a1)
i022 = besseli(0,w22); % i022 = besseli(0,w2*a2)
i122 = besseli(1,w22); % i122 = besseli(1,w2*a2)
i222 = besseli(2,w22); % i222 = besseli(2,w2*a2)
k021 = besselk(0,w21); % k021 = besselk(0,w2*a1)
k121 = besselk(1,w21); % k121 = besselk(1,w2*a1)
k221 = besselk(2,w21); % k221 = besselk(2,w2*a1)
k022 = besselk(0,w22); % k022 = besselk(0,w2*a2)
k122 = besselk(1,w22); % k122 = besselk(1,w2*a2)
k222 = besselk(2,w22); % k222 = besselk(2,w2*a2)
k032 = besselk(0,w32); % k032 = besselk(0,w3*a2)
k132 = besselk(1,w32); % k132 = besselk(1,w3*a2)
k232 = besselk(2,w32); % k232 = besselk(2,w3*a2)
% besselj'(n,x) = 0.5*(besselj(n-1,x)-besselj(n+1,x))
% besseli'(n,x) = 0.5*(besseli(n-1,x)+besseli(n+1,x))
j11 = 0.5*(j011-j211); % besselj'(1,u1*a1) = 0.5*(besselj(0,u1*a1)-besselj(2,u1*a1))
i21 = 0.5*(i021+i221); % besseli'(1,w2*a1) = 0.5*(besseli(0,w2*a1)+besseli(2,w2*a1))
i22 = 0.5*(i022+i222); % besseli'(1,w2*a2) = 0.5*(besseli(0,w2*a2)+besseli(2,w2*a2))
k21 = -0.5*(k021+k221); % besselk'(1,w2*a1) = -0.5*(besselk(0,w2*a1)+besselk(2,w2*a1))
k22 = -0.5*(k022+k222); % besselk'(1,w2*a2) = -0.5*(besselk(0,w2*a2)+besselk(2,w2*a2))
k32 = -0.5*(k032+k232); % besselk'(1,w3*a2) = -0.5*(besselk(0,w3*a2)+besselk(2,w3*a2))
p1 = i122*k121-i121*k122;
q1 = i122*k21-i21*k122;
r1 = i22*k121-i121*k22;
s1 = i22*k21-i21*k22;
kb = k32/(bw3*k132); % Kb = besselk'(1,w3*a2)/(W3*besselk(1,w3*a2))
jus = bus1*j111*j111; % jus = U1^2*(besselj(2,u1*a1))^2
ju = sqrt(jus);
aw = 1/(ar*bw2); % au = 1/(alpha2*U2)
aws1 = ars*bws2; % aus1 = alpha2^2*U2^2
qw1 = q1/bw2;
rw1 = r1*aw;
sw1 = s1*aw;
was2 = 4/(pi*pi*aws1*bws2);
jp = j11*p1; % jp = besselj'(1,u1*a1)*p1
kp = kb*p1; % kp = k3*p2 = K*p1
kq = kb*q1; % kq = k3*q2 = K*q1
psq = jp-s21*qw1*ju;
psr = kp+s23*rw1;
pq = jp-qw1*ju;
pr = kp+rw1;
deltas = nes; % deltas = delta^2 = (neff)^2
delta = sqrt(deltas);
x1 = n1*bus1*bws2/(delta*bvs12);
x2 = n3*ars*bws2*bws3/(delta*bvs23);
xs1 = x1*x1;
xs2 = x2*x2;
f1 = (p1*p1-2*ns2*x1*x2/(aws1*bws2))*jus;
f2 = xs1*xs2*(j11*(kp-rw1)+ju*(kq-sw1)/bw2);
f3 = (j11*(kp-s23*rw1)+ju*s21*(kq-s23*sw1)/bw2);
f4 = -xs1*psq*pq;
f5 = -xs2*psr*pr*jus;
f= f1+f2*f3+f4+f5;
另外编写了一个ecothp_HE.m文件求解上面的方程，首先绘函数f的图形，找到f=0点的大致区间，再求根，具体代码如下：
%-------------------------------------------------------------------------
% coremode_HE11:eigenvalue of core modes of three layer fiber ,ploting f
%--------------------------------------------------------------------------
clear all
close all
clc
format long
global lambda n1 n2 n3 ra1 ra2 U_HE11
n1 = 1.4681;
n2 = 1.4628;
n3 = 1;
ra1 = 4.15e-6;
ra2 = 62.5e-6;
lambda = 1.55e-6;
ns1 = n1*n1;
ns2 = n2*n2;
ns3 = n3*n3;
kw = 2*pi/lambda;
kws = kw*kw;
kr1 = kw*ra1;
krs1 = kr1*kr1;
bvs1 = krs1*(ns1-ns2); % v1^2 = k0^2*ra1^2*(n1^2-n2^2)
bv1 = sqrt(bvs1);
N = 1000; % the total caculation counts
Umax = bv1; % For the three-layer SMF,the effective index(ne) of
% the core modes should lie between n1 and n2,
% i.e. n1 > ne > n2.So,U1 lies between bv1 and
% bv2,i.e. bv1 > U1 > 0.
dU = Umax/N;
U = 0; % initial value of U1
for k = 1:N
U = U+dU; % U
xu(k) = U;
fcoHE = ecothf_HE(U);
fco(k) = fcoHE;
end
figure()
plot(xu,fco)
xlabel('U')
ylabel('F(U)')
title('Core-mode HE11 of three-layer fiber','fontweight','bold')
grid on
figure()
fplot(@ecothf_HE,[0,1])
title('HE11 mode between [0,1]','fontweight','bold')
xlabel('U')
ylabel('F(U)')
grid on;
U_HE11 = fzero(@ecothf_HE,[0.3,0.4])
U = U_HE11; % The first root of eigenvalue equation corresponds to the HE11 mode
parameters(U)
a = fzero(@ecothf_HE,[1.4902,1.4904])
parameters(a)
复制代码
以上代码中调用了parameters.m文件，这是求光纤模场其他的一些参数，如neff，beta等，很简单的代码，如下：
function F = parameters(U1)
format long
global lambda n1 n2 n3 ra1 ra2
n1 = 1.4681;
n2 = 1.4628;
n3 = 1;
ns1 = n1*n1;
ns2 = n2*n2;
ns3 = n3*n3;
ra1 = 4.15e-6;
ras1 = ra1*ra1;
ra2 = 62.5e-6;
lambda = 1.55e-6;
kw = 2*pi/lambda;
kws = kw*kw;
kr1 = kw*ra1;
kr2 = kw*ra2;
krs1 = kr1*kr1;
krs2 = kr2*kr2;
vs1 = krs1*(ns1-ns2); % v1^2 = k0^2*ra1^2*(n1^2-n2^2)
vs2 = krs1*(ns2-ns3); % v2^2 = k0^2*ra1^2*(n2^2-n3^2)
v1 = sqrt(vs1);
v2 = sqrt(vs2);
bu = U1;
bus = bu*bu;
u1 = bu/ra1
neff = sqrt(ns1-bus/krs1)
beta = kw*neff
运行ecothp_HE.m后，
从函数f的图形中可以看出，f=0有两个点，也就是f=0有两个根，Corning SMF单模光纤应该只有一个纤芯模，即f=0应该只有一个根，该选择哪一个根作为纤芯模的参数呢？？？
然后，假设选择第一个根作为光纤参数，即U_HE11 = 1.561149962979787，再计算纤芯模模场分布归一化系数Eco，参考文献为：(见2楼)
[Last edited by syulzh on 2012-11-22 at 09:00]
以上问题均已解决，若有光纤方面的同仁，请加本人qq交流，共同进步。
[Last edited by syulzh on 2013-9-30 at 19:08]