1、光纤模式图 matlab 模拟*摘要:光纤通信是现代化通信的支柱,在光纤通信中 ,光纤是最重要的部件之一。本文利用电磁波动理论推导了光在光纤中的传输模式的本征方程,并使用 Matlab软件绘出不同条件下的模式图.关键词:光纤模式;电磁波动;Matlab1、引言对光纤中光的传播理论的研究,可以有多种方法,比如射线法,标量近似分析法等,但为了更广泛地描述光纤波导中光的传播,更详细地研究光纤的传输特性,就必须运用波动光学理论对光纤进行分析.本文从麦克斯韦方程的求解出发推导光纤的传播模式本征方程并利用Matlab模拟其模式图.要对光在光纤中的传播特性有详细的理解,必须依靠麦克斯韦方程,结合问题中的边界
2、条件,求解电磁矢量场.求解的方法一般是:1 、先求出亥姆霍兹方程组以及电磁场纵向分量E z和H z的具体形式.2、把E z和H z有具体形式代入麦克斯韦方程以求取其他电磁场横向分量 、E r、 、H r.3、利用界面上电磁场 和 切 EH向连续条件，求取模式本征方程 1.2、波动方程由麦克斯韦方程组,我们知道,光纤中电磁场的波动方程可以写成:(1) 式中参量表示介质的介电常数,表示介质的磁导率 .对于在圆柱形光纤中传播的电磁波.电场和磁场具有如下形式的函数关系：(2)式中为光纤中导波沿 z 轴方向的传播常数,其值由纤芯包层界面处的电磁场边界条件决定.不同的值对应于不同阶的导波模式,它们的场分布
3、也不同.将式(2)代入波动方程式 (1)中,可得到矢量亥姆霍兹方程,即(3)在柱坐标系中,只有沿z 轴方向的单位矢量与场点位置无关，所以，在柱坐标系Et22Ht22)()(),( tzjeryxE )(),(),( tzjeryx022Ekn 022Hkn中，只有Ez和Hz 才满足标量亥姆霍兹方程 ,可得(4)解得方程(4), 可得到电场和磁场的纵向分量Ez 和Hz.(5)(6)式中(ra)表示在纤芯内部, 称为归一化横向传播常数.其大小随纤芯内场的不同模式而变.(ra)表示在纤芯的外部区域, 称为归一化横向衰减常数,其大小及符号反映了包层中场的状态.有了场分量Ez和Hz的表达式,再利用麦氏方
4、程组即可求出场的其它四个分量 、E r、 、H r.3、本征方程及模式图3.1 本征方程光纤中传播模式及传输特性都是由它自身的本征方程确定的. 在光纤的基本参量n1,n2,a ,k0已知的条件下, U,W仅与传播常数有关.用所导出的各个区域中电磁场的表达式,再利用电磁场切向分量在纤芯-包层界面上(r=a)连续的条件,就可以救出模式本征方程,也称特征方程.(7)3.2 各类模式根据Jm(u) 的振荡特性,对于一特定的m值,本征方程存在着 n个根.当m=0 时，若 ,相应于 ,只有磁场纵向分量.(8)当m=0 时，若 ,相应于 ,只有磁场纵向分量.(9)当m0 时，混合模式 HE模和 EH模0)(
5、)1)(1)( 22022 zzzz HEnkHErErHr );()() areraWKAUJEimmiz );()() areraKBJUHimmiz 210nk 20nkW220420 221)1()( )()()()( WUkmWVk KnJnKUJ mm BA,zzHE)(0)(1)(00 模TEKUJ,zz )()()(02021 模TMWnJ )1)(1(4)1(2)1(222212 WUnmKnKnJ (10)其中:3.3 Matlab模拟流程图及模式图为了分析导波模的传输特性,就需要得知各模式传播常数随光纤归一化频率V的变化情况.这可通过对本征方程(10)求解而得出 .其解可
6、写为(11)式中方程(11) 是超越方程,在截止和远离截止的情况下,可以将它简化成简单的形式求解,得出各种矢量模式的截止频率Vc和远离截止时的u值,从而进行传输特性的分析.而在一般情况下(不局限于截止和远离截止两种状态) ,就需使用计算机对特征方程(7)求数值解 .其计算流程图如图1所示.图2分别给出了依照此流程绘制出的TE模、HE 11模及EH 11模的模式图.图1 计算/ k0-V 曲线的程序流程图)(,)(WKUJmm0)1(22FnJ )1)(1(422211 WUnmKF 四、结束语光纤中的传播理论已为人们充分了解,本文直接由麦克斯韦方程组出发,精确地求出电磁场各个分量,根据光纤芯包
7、层界面处电磁场的边界条件, 在计算机上通过数值求解,得到传播常数及光纤中的场分布模式,因而对光纤中传导模的描述完整,结果明确 4.五、参考文献：1佘守宪等. 导波光学物理基础M.北京:北方交通大学出版社,2002.2陈军等.光学电磁理论 M.北京:科学出版社,2005.3陈抗生等. 微波与光导波技术教程M.浙江:浙江大学出版社,2000.4薛苏云等 .阶跃折射率光纤的电磁场模式研究J.河海大学常州分校学报.2000.14(4):16-20.1EH0/k210nakV(a)01TE0/k210nakV(b)1HE0/210nakV(c)图 2 / K0 - V 曲线模式图(a)EH11 模(b)
8、TE01 模(c)HE11 模及附录一:Matlab实现TE模程序format longclearn1=1.45n2=1.447a=4.5k=(a2)*(n12-n22);u1=(a2)*(n12);w1=(a2)*(n22);delta1=(n12+n22)/(2*n12);delta2=(n12-n22)/(2*n12);i=1;n=n2;for V=0:0.01:6k02=V2/k;k01=sqrt(k02);for BeiTa=n:0.00001:n1%nnn=n+0.00001U2=u1*k02-(a2)*k02*(BeiTa2);U=sqrt(U2);W2=(a2)*k02*(Be
9、iTa2)-w1*k02;W=sqrt(W2);if(U=0 | W=0)break;%disp(sss)elsez1=n12*besselj(1,U)/(U*besselj(0,U);z2=n22*besselk(1,W)/(W*besselk(0,W);z3=1/(U2)+delta1*(1/W2-z2)-sqrt(delta22*(1/W2-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;endendend endplot(x,y);ax
10、is(0 6.5 1.447 1.45);附录二:Matlab实现HE 11模程序format longclearn1=1.45n2=1.447a=4.5k=(a2)*(n12-n22);u1=(a2)*(n12);w1=(a2)*(n22);delta1=(n12+n22)/(2*n12);delta2=(n12-n22)/(2*n12);i=1;n=n2;for V=0:0.001:1k02=V2/k;k01=sqrt(k02);for BeiTa=n:0.00001:n1%nnn=n+0.000001U2=u1*k02-(a2)*k02*(BeiTa2);U=sqrt(U2)W2=(a2
11、)*k02*(BeiTa2)-w1*k02;W=sqrt(W2)if(U=0 | W=0)break;%disp(sss)elsez1=besselj(0,U)/(U*besselj(1,U);z2=besselk(0,W)/(W*besselk(1,W);z3=1/(U2)+delta1*(1/W2+z2)+sqrt(delta22*(1/W2+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;endendend endfor V=1:0.01
12、:6k02=V2/k;k01=sqrt(k02);for BeiTa=n:0.00001:n1%nU2=u1*k02-(a2)*k02*(BeiTa2);U=sqrt(U2);W2=(a2)*k02*(BeiTa2)-w1*k02;W=sqrt(W2);if(U=0 | W=0)break;elsez1=besselj(0,U)/(U*besselj(1,U);z2=besselk(0,W)/(W*besselk(1,W);z3=1/(U2)+delta1*(1/W2+z2)+sqrt(delta22*(1/W2+z2)2+(BeiTa/n1)2*(V/(U*W)4);if(abs(z1-z3
13、)0.01)x(i)=V;y(i)=BeiTa;i=i+1;z1;z2;n=BeiTa;break;endendend endplot(x,y);axis(0 6.5 1.447 1.45);附录三:Matlab 实现EH 11模程序format longclearn1=1.45n2=1.447a=4.5k=(a2)*(n12-n22);u1=(a2)*(n12);w1=(a2)*(n22);delta1=(n12+n22)/(2*n12);delta2=(n12-n22)/(2*n12);i=1;n=n2;for V=0:0.001:1k02=V2/k;k01=sqrt(k02);for B
14、eiTa=n:0.00001:n1%nnn=n+0.000001U2=u1*k02-(a2)*k02*(BeiTa2);U=sqrt(U2)W2=(a2)*k02*(BeiTa2)-w1*k02;W=sqrt(W2)if(U=0 | W=0)break;%disp(sss)elsez1=besselj(0,U)/(U*besselj(1,U);z2=besselk(0,W)/(W*besselk(1,W);z3=1/(U2)+delta1*(1/W2+z2)+sqrt(delta22*(1/W2+z2)2+(BeiTa/n1)2*(V/(U*W)4);if(abs(z1-z3)0.01)x(i
15、)=V;y(i)=BeiTa;i=i+1;z1;z2;%disp(nnn)n=BeiTa;break;endendend endfor V=1:0.01:6k02=V2/k;k01=sqrt(k02);for BeiTa=n:0.00001:n1%nU2=u1*k02-(a2)*k02*(BeiTa2);U=sqrt(U2);W2=(a2)*k02*(BeiTa2)-w1*k02;W=sqrt(W2);if(U=0 | W=0)break;elsez1=besselj(0,U)/(U*besselj(1,U);z2=besselk(0,W)/(W*besselk(1,W);z3=1/(U2)+delta1*(1/W2+z2)+sqrt(delta22*(1/W2+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;endendend endplot(x,y);axis(0 6.5 1.447 1.45);
