光纖的模式MATLAB圖模擬

1、光纖模式圖matlab模擬 *摘要:光纖通信是現代化通信的支柱,在光纖通信中,光纖是最重要的部件之一。本文利用電磁波動理論推導了光在光纖中的傳輸模式的本征方程,并使用Matlab軟件繪出不同條件下的模式圖.關鍵詞:光纖模式;電磁波動;Matlab1、 引言對光纖中光的傳播理論的研究,可以有多種方法,比如射線法,標量近似分析法等,但為了更廣泛地描述光纖波導中光的傳播,更詳細地研究光纖的傳輸特性,就必須運用波動光學理論對光纖進行分析.本文從麥克斯韋方程的求解出發(fā)推導光纖的傳播模式本征方程并利用Matlab模擬其模式圖.要對光在光纖中的傳播特性有詳細的理解,必須依靠麥克斯韋方程,結合問題中的邊界條件

2、,求解電磁矢量場.求解的方法一般是:1、先求出亥姆霍茲方程組以及電磁場縱向分量Ez和Hz的具體形式.2、把Ez和Hz有具體形式代入麥克斯韋方程以求取其他電磁場橫向分量、Er、Hr.3、利用界面上電磁場和切向連續(xù)條件，求取模式本征方程1.2、 波動方程由麥克斯韋方程組,我們知道,光纖中電磁場的波動方程可以寫成: (1) 式中參量表示介質的介電常數,表示介質的磁導率.對于在圓柱形光纖中傳播的電磁波.電場和磁場具有如下形式的函數關系： (2)式中為光纖中導波沿z 軸方向的傳播常數,其值由纖芯包層界面處的電磁場邊界條件決定.不同的值對應于不同階的導波模式,它們的場分布也不同.將式(2)代入波動方程式(

3、1)中,可得到矢量亥姆霍茲方程,即 (3)在柱坐標系中,只有沿z軸方向的單位矢量與場點位置無關，所以，在柱坐標系中，只有Ez和Hz才滿足標量亥姆霍茲方程,可得 (4)解得方程(4),可得到電場和磁場的縱向分量Ez和Hz. (5) (6)式中(ra)表示在纖芯內部, 稱為歸一化橫向傳播常數.其大小隨纖芯內場的不同模式而變.(ra)表示在纖芯的外部區(qū)域, 稱為歸一化橫向衰減常數,其大小及符號反映了包層中場的狀態(tài).有了場分量Ez和Hz的表達式,再利用麥氏方程組即可求出場的其它四個分量、Er、Hr.3、 本征方程及模式圖3.1 本征方程光纖中傳播模式及傳輸特性都是由它自身的本征方程確定的. 在光纖的基

4、本參量n1,n2,a ,k0已知的條件下, U,W僅與傳播常數有關.用所導出的各個區(qū)域中電磁場的表達式,再利用電磁場切向分量在纖芯-包層界面上(r=a)連續(xù)的條件,就可以救出模式本征方程,也稱特征方程. (7)3.2 各類模式 根據Jm(u)的振蕩特性,對于一特定的m值,本征方程存在著n個根. 當m=0時，若,相應于 ,只有磁場縱向分量. (8)當m=0時，若 ,相應于 ,只有磁場縱向分量. (9)當m0時，混合模式HE模和EH模 (10)其中:3.3 Matlab模擬流程圖及模式圖 為了分析導波模的傳輸特性,就需要得知各模式傳播常數隨光纖歸一化頻率V的變化情況.這可通過對本征方程(10)求解

5、而得出.其解可寫為 (11)式中 方程(11)是超越方程,在截止和遠離截止的情況下,可以將它簡化成簡單的形式求解,得出各種矢量模式的截止頻率Vc和遠離截止時的u值,從而進行傳輸特性的分析.而在一般情況下(不局限于截止和遠離截止兩種狀態(tài)) ,就需使用計算機對特征方程(7)求數值解.其計算流程圖如圖1所示.圖2分別給出了依照此流程繪制出的TE模、HE11模及EH11模的模式圖. 圖1 計算/ k0-V曲線的程序流程圖(b)(a)圖2/ K0 - V曲線模式圖 (a)EH11模 (b)TE01模 (c)HE11模及(c)四、結束語 光纖中的傳播理論已為人們充分了解,本文直接由麥克斯韋方程組出發(fā),精確

6、地求出電磁場各個分量,根據光纖芯包層界面處電磁場的邊界條件, 在計算機上通過數值求解,得到傳播常數及光纖中的場分布模式,因而對光纖中傳導模的描述完整,結果明確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-

7、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*bes

8、selj(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 longcl

9、earn1=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

10、) 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)

11、; 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)

12、=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

13、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

14、(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=bes