基于matlab光纤的模式图模拟.doc

. . . . .基于matlab的光纤模式图模拟摘要:光纤通信是现代化通信的支柱,在光纤通信中,光纤是最重要的部件之一。本文利用电磁波动理论推导了光在光纤中的传输模式的本征方程,并使用Matlab软件绘出不同条件下的模式图.关键词:光纤模式;电磁波动;Matlab1、 引言对光纤中光的传播理论的研究,可以有多种方法,比如射线法,标量近似分析法等,但为了更广泛地描述光纤波导中光的传播,更详细地研究光纤的传输特性,就必须运用波动光学理论对光纤进行分析.本文从麦克斯韦方程的求解出发推导光纤的传播模式本征方程并利用Matlab模拟其模式图.要对光在光纤中的传播特性有详细的理解,必须依靠麦克斯韦方程,结合问题中的边界条件,求解电磁矢量场.求解的方法一般是:1、先求出亥姆霍兹方程组以及电磁场纵向分量Ez和Hz的具体形式.2、把Ez和Hz有具体形式代入麦克斯韦方程以求取其他电磁场横向分量、Er、Hr.3、利用界面上电磁场和切向连续条件，求取模式本征方程1.2、 波动方程由麦克斯韦方程组,我们知道,光纤中电磁场的波动方程可以写成: (1) 式中参量表示介质的介电常数,表示介质的磁导率.对于在圆柱形光纤中传播的电磁波.电场和磁场具有如下形式的函数关系： (2)式中为光纤中导波沿z 轴方向的传播常数,其值由纤芯包层界面处的电磁场边界条件决定.不同的值对应于不同阶的导波模式,它们的场分布也不同.将式(2)代入波动方程式(1)中,可得到矢量亥姆霍兹方程,即 (3)在柱坐标系中,只有沿z轴方向的单位矢量与场点位置无关，所以，在柱坐标系中，只有Ez和Hz才满足标量亥姆霍兹方程,可得 (4)解得方程(4),可得到电场和磁场的纵向分量Ez和Hz. (5) (6)式中(ra)表示在纤芯内部, 称为归一化横向传播常数.其大小随纤芯内场的不同模式而变.(ra)表示在纤芯的外部区域, 称为归一化横向衰减常数,其大小及符号反映了包层中场的状态.有了场分量Ez和Hz的表达式,再利用麦氏方程组即可求出场的其它四个分量、Er、Hr.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模 (10)其中:3.3 Matlab模拟流程图及模式图 为了分析导波模的传输特性,就需要得知各模式传播常数随光纤归一化频率V的变化情况.这可通过对本征方程(10)求解而得出.其解可写为 (11)式中 方程(11)是超越方程,在截止和远离截止的情况下,可以将它简化成简单的形式求解,得出各种矢量模式的截止频率Vc和远离截止时的u值,从而进行传输特性的分析.而在一般情况下(不局限于截止和远离截止两种状态) ,就需使用计算机对特征方程(7)求数值解.其计算流程图如图1所示.图2分别给出了依照此流程绘制出的TE模、HE11模及EH11模的模式图. 图1 计算/ k0-V曲线的程序流程图(b)(a)图2/ K0 - V曲线模式图 (a)EH11模 (b)TE01模 (c)HE11模及(c)四、结束语 光纤中的传播理论已为人们充分了解,本文直接由麦克斯韦方程组出发,精确地求出电磁场各个分量,根据光纤芯包层界面处电磁场的边界条件, 在计算机上通过数值求解,得到传播常数及光纤中的场分布模式,因而对光纤中传导模的描述完整,结果明确4.五、参考文献：1佘守宪等.导波光学物理基础M.北京:北方交通大学出版社,2002.2陈军等.光学电磁理论M.北京:科学出版社,2005.3陈抗生等.微波与光导波技术教程M.浙江:浙江大学出版社,2000.4薛苏云等.阶跃折射率光纤的电磁场模式研究J.河海大学常州分校学 报.2000.14(4):16-20.附录一: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:6 k02=V2/k; k01=sqrt(k02); for BeiTa=n:0.00001:n1 %n nn=n+0.00001 U2=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) else z1=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; end end end endplot(x,y);axis(0 6.5 1.447 1.45);附录二:Matlab实现HE11模程序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:1 k02=V2/k; k01=sqrt(k02); for BeiTa=n:0.00001:n1 %n nn=n+0.000001 U2=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) else z1=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; end end end endfor V=1:0.01:6 k02=V2/k; k01=sqrt(k02); for BeiTa=n:0.00001:n1 %n U2=u1*k02-(a2)*k02*(BeiTa2); U=sqrt(U2); W2=(a2)*k02*(BeiTa2)-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/(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; end end end end plot(x,y);axis(0 6.5 1.447 1.45);附录三:Matlab实现EH11模程序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:1 k02=V2/k; k01=sqrt(k02); for BeiTa=n:0.00001:n1 %n nn=n+0.000001 U2=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) else z1=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; end end end endfor V=1:0.01:6 k02=V2/k; k01=sqrt(k02); for BeiTa=n:0.00001:n1 %n U2=u1*k02-(a2)*k02*(BeiTa2); U=sqrt(U2); W2=(a2)*k02*(BeiTa2)-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/(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=Bei