光纤的模式MATLAB图模拟

光纤模式图matlab 模拟
**************
摘要:光纤通信是现代化通信的支柱,在光纤通信中,光纤是最重要的部件之一。

本文利用电磁波动理论推导了光在光纤中的传输模式的本征方程,并使用Matlab 软件绘出不同条件下的模式图.
关键词:光纤模式;电磁波动;Matlab
一、引言
对光纤中光的传播理论的研究,可以有多种方法,比如射线法,标量近似分析法等,但为了更广泛地描述光纤波导中光的传播,更详细地研究光纤的传输特性,就必须运用波动光学理论对光纤进行分析.本文从麦克斯韦方程的求解出发推导光纤的传播模式本征方程并利用Matlab 模拟其模式图.
要对光在光纤中的传播特性有详细的理解,必须依靠麦克斯韦方程,结合问题中的边界条件,求解电磁矢量场.求解的方法一般是:1、先求出亥姆霍兹方程组以及电磁场纵向分量E z 和H z 的具体形式.2、把E z 和H z 有具体形式代入麦克斯韦方程以求取其他电磁场横向分量θE 、E r 、θH 、H r .3、利用界面上电磁场θE 和θH 切
向连续条件，求取模式本征方程[1].
二、波动方程
由麦克斯韦方程组,我们知道,光纤中电磁场的波动方程可以写成: (1)
式中参量ε表示介质的介电常数,μ表示介质的磁导率.
对于在圆柱形光纤中传播的电磁波.电场和磁场具有如下形式的函数关系：
(2)
式中β为光纤中导波沿z 轴方向的传播常数,其值由纤芯———包层界面处的电磁场边界条件决定.不同的β值对应于不同阶的导波模式,它们的场分布也不同.
将式(2)代入波动方程式(1)中,可得到矢量亥姆霍兹方程,即 (3)
在柱坐标系中,只有沿z 轴方向的单位矢量与场点位置无关，所以，在柱坐标系中，
E
t
E 222∂∂=∇εμH t H 222∂∂=∇εμ)(),(),(t z j e r E y x E ωβθ-=)(),(),(t z j e r H
y x H ωβθ-=02022=+∇E k n E 02022=+∇H k n H
只有Ez 和Hz 才满足标量亥姆霍兹方程,可得 (4)
解得方程(4),可得到电场和磁场的纵向分量Ez 和Hz.
(5)
(6)
式中(r≤a)表示在纤芯内部, 称为归一化横向传播常数.其大小随纤
芯内场的不同模式而变.(r>a)表示在纤芯的外部区域, 称为归一
化横向衰减常数,其大小及符号反映了包层中场的状态.有了场分量Ez 和Hz 的表达式,再利用麦氏方程组即可求出场的其它四个分量θE 、E r 、θH 、H r.
三、本征方程及模式图 3.1 本征方程
光纤中传播模式及传输特性都是由它自身的本征方程确定的. 在光纤的基本参量n1,n2,a ,k0已知的条件下, U,W 仅与传播常数β有关.用所导出的各个区域中电磁场的表达式,再利用电磁场切向分量在纤芯-包层界面上(r=a)连续的条件,就可以救出模式本征方程,也称特征方程.
(7)
3.2 各类模式
根据Jm(u)的振荡特性,对于一特定的m 值,本征方程存在着n 个根.
当m=0时，若 ,相应于 ,只有磁场纵向分量. (8)
当m=0时，若 ,相应于 ,只有磁场纵向分量. (9)
当m>0时，混合模式HE 模和EH 模
0))(()(1)(1)(2
220
22222=-+∂∂+∂∂+∂∂z z z
z z z z z H E n k H E r H E r r H E r βθ⎪⎪⎩⎪⎪
⎨⎧>≤=);
()()
();()()(a r e r a W K W K A a r e r a U J U J A E im m m im m m z
θθ⎪⎪⎩⎪⎪⎨⎧>≤=);
()()();()()(a r e r a W K W K B a r e r a U J U J B H im m m im m m
z θθ2
2120β-=n k a U 22220β
-=n k a W 222204202221)11()()()(])
()()()(][)()()()([W U k m UW V k m W WK W K n U UJ U J n W WK W K U UJ U J m m
m m m m m m +=='+''+'ββ0,0≠=B A 0,0≠=z z H E )
(0)()(1)()(10
000
模TE W K W K W U J U J U ='+'0,0=≠B A 0
,0=≠z z H E )(0)()()()(002
20021模TM W K W K W n U J U J U n ='+')]
11)(11(44[)1(21)1(21222222122
221212*********W
U W n n U m K n n K n n K n n J ++--+±+-=
(10)
其中:
3.3 Matlab 模拟流程图及模式图
为了分析导波模的传输特性,就需要得知各模式传播常数β随光纤归一化频率V 的变化情况.这可通过对本征方程(10)求解而得出.其解可写为
(11)
式中
方程(11)是超越方程,在截止和远离截止的情况下,可以将它简化成简单的形
式求解,得出各种矢量模式的截止频率Vc 和远离截止时的u 值,从而进行传输特性的分析.而在一般情况下(不局限于截止和远离截止两种状态) ,就需使用计算机对特征方程(7)求数值解.其计算流程图如图1所示.图2分别给出了依照此流程绘制出的TE 模、HE 11模及EH 11模的模式图.
图1 计算β/ k0-V 曲线的程序流程图
)
()(,)()(W WK W K K U UJ U J J m
m
m m '='= 0])1[(2121
22
=±++F K n n J )]
11)(11(44[)1(2
1222222122
221
21222121W
U W
n n U
m K n n K n n F ++--+±=
四、结束语
光纤中的传播理论已为人们充分了解,本文直接由麦克斯韦方程组出发,精确地求出电磁场各个分量,根据光纤芯———包层界面处电磁场的边界条件, 在计算机上通过数值求解,得到传播常数及光纤中的场分布模式,因而对光纤中传导模的描述完整,结果明确[4].
五、参考文献：
[1]佘守宪等.导波光学物理基础[M].北京:北方交通大学出版社,2002. [2]陈军等.光学电磁理论[M].北京:科学出版社,2005.
[3]陈抗生等.微波与光导波技术教程[M].浙江:浙江大学出版社,2000. [4]薛苏云等.阶跃折射率光纤的电磁场模式研究[J].河海大学常州分校学 报.2000.14(4):16-20.
11
EH 0/k ββ=2
2210n n a k V -=(a)
01
TE 0
/k ββ=2
2
210n n a k V -=(b)
11
HE 0/k ββ=2
2210n n a k V -=(c)
图2 β/ K0 - V 曲线模式图 (a)EH 11模 (b)TE 01模 (c)HE 11模及
附录一:Matlab实现TE模程序
format long
clear
n1=1.45
n2=1.447
a=4.5
k=(a^2)*(n1^2-n2^2);
u1=(a^2)*(n1^2);
w1=(a^2)*(n2^2);
delta1=(n1^2+n2^2)/(2*n1^2);
delta2=(n1^2-n2^2)/(2*n1^2);
i=1;
n=n2;
for V=0:0.01:6
k02=V^2/k;
k01=sqrt(k02);
for BeiTa=n:0.00001:n1
%n
nn=n+0.00001
U2=u1*k02-(a^2)*k02*(BeiTa^2);
U=sqrt(U2);
W2=(a^2)*k02*(BeiTa^2)-w1*k02;
W=sqrt(W2);
if(U==0 || W==0)
break;
%disp('sss')
else
z1=n1^2*besselj(1,U)/(U*besselj(0,U));
z2=n2^2*besselk(1,W)/(W*besselk(0,W));
z3=1/(U^2)+delta1*(1/W^2-z2)-sqrt(delta2^2*(1/W^2-z2)^2+(BeiTa/n1)^2*(V/(U* W))^4);
if(abs(z1+z2)<0.01)
x(i)=V;
y(i)=BeiTa;
i=i+1;
z1;
z2;
%disp('nnn')
n=BeiTa;
break;
end
end
end
end
plot(x,y);
axis([0 6.5 1.447 1.45]);
附录二:Matlab实现HE11模程序
format long
clear
n1=1.45
n2=1.447
a=4.5
k=(a^2)*(n1^2-n2^2);
u1=(a^2)*(n1^2);
w1=(a^2)*(n2^2);
delta1=(n1^2+n2^2)/(2*n1^2);
delta2=(n1^2-n2^2)/(2*n1^2);
i=1;
n=n2;
for V=0:0.001:1
k02=V^2/k;
k01=sqrt(k02);
for BeiTa=n:0.00001:n1
%n
nn=n+0.000001
U2=u1*k02-(a^2)*k02*(BeiTa^2);
U=sqrt(U2)
W2=(a^2)*k02*(BeiTa^2)-w1*k02;
W=sqrt(W2)
if(U==0 || W==0)
break;
%disp('sss')
else
z1=besselj(0,U)/(U*besselj(1,U));
z2=besselk(0,W)/(W*besselk(1,W));
z3=1/(U^2)+delta1*(1/W^2+z2)+sqrt(delta2^2*(1/W^2+z2)^2+(BeiTa/n1)^2*(V/(U* W))^4);
if(abs(z1-z3)<0.01)
x(i)=V;
y(i)=BeiTa;
i=i+1;
z1;
z2;
%disp('nnn')
n=BeiTa;
break;
end
end
end
end
for V=1:0.01:6
k02=V^2/k;
k01=sqrt(k02);
for BeiTa=n:0.00001:n1
%n
U2=u1*k02-(a^2)*k02*(BeiTa^2);
U=sqrt(U2);
W2=(a^2)*k02*(BeiTa^2)-w1*k02;
W=sqrt(W2);
if(U==0 || W==0)
break;
else
z1=besselj(0,U)/(U*besselj(1,U));
z2=besselk(0,W)/(W*besselk(1,W));
z3=1/(U^2)+delta1*(1/W^2+z2)+sqrt(delta2^2*(1/W^2+z2)^2+(BeiTa/n1)^2*(V/(U* W))^4);
if(abs(z1-z3)<0.01)
x(i)=V;
y(i)=BeiTa;
i=i+1;
z1;
z2;
n=BeiTa;
break;
end
end
end
end
plot(x,y);
axis([0 6.5 1.447 1.45]);
附录三:Matlab实现EH11模程序
format long
clear
n1=1.45
n2=1.447
a=4.5
k=(a^2)*(n1^2-n2^2);
u1=(a^2)*(n1^2);
w1=(a^2)*(n2^2);
delta1=(n1^2+n2^2)/(2*n1^2);
delta2=(n1^2-n2^2)/(2*n1^2);
i=1;
n=n2;
for V=0:0.001:1
k02=V^2/k;
k01=sqrt(k02);
for BeiTa=n:0.00001:n1
%n
nn=n+0.000001
U2=u1*k02-(a^2)*k02*(BeiTa^2);
U=sqrt(U2)
W2=(a^2)*k02*(BeiTa^2)-w1*k02;
W=sqrt(W2)
if(U==0 || W==0)
break;
%disp('sss')
else
z1=besselj(0,U)/(U*besselj(1,U));
z2=besselk(0,W)/(W*besselk(1,W));
z3=1/(U^2)+delta1*(1/W^2+z2)+sqrt(delta2^2*(1/W^2+z2)^2+(BeiTa/n1)^2*(V/(U* W))^4);
if(abs(z1-z3)<0.01)
x(i)=V;
y(i)=BeiTa;
i=i+1;
z1;
z2;
%disp('nnn')
n=BeiTa;
break;
end
end
end
end
for V=1:0.01:6
k02=V^2/k;
k01=sqrt(k02);
for BeiTa=n:0.00001:n1
%n
U2=u1*k02-(a^2)*k02*(BeiTa^2);
U=sqrt(U2);
W2=(a^2)*k02*(BeiTa^2)-w1*k02;
W=sqrt(W2);
if(U==0 || W==0)
break;
else
z1=besselj(0,U)/(U*besselj(1,U));
z2=besselk(0,W)/(W*besselk(1,W));
z3=1/(U^2)+delta1*(1/W^2+z2)+sqrt(delta2^2*(1/W^2+z2)^2+(BeiTa/n1)^2*(V/(U* W))^4);
if(abs(z1-z3)<0.01)
x(i)=V;
y(i)=BeiTa;
i=i+1;
z1;
z2;
n=BeiTa;
break;
end
end
end
end
plot(x,y);
axis([0 6.5 1.447 1.45]);。