光纖光柵理論基礎(chǔ)與傳感技術(shù)-求解三層光纖包層模本征方程

1、精選優(yōu)質(zhì)文檔-傾情為你奉上%-% 求解三層光纖包層模本征方程,畫(huà)曲線f(U1),% eclthp:eigenvalue of cladding modes of three layer fiber ,ploting f % 光纖光柵理論基礎(chǔ)與傳感技術(shù) % zhzijia %- clear all,clcformat longglobal lambda n1 n2 n3 ra1 ra2 %更改方便,引入全局量n1 = 1.4681;n2 = 1.4628;n3 = 1;ns1 = n1*n1;ns2 = n2*n2;ns3 = n2*n3; ra1 = 4.15e-6;ra2 = 62.5e-6

2、; lambda = 1550e-9; kw = 2*pi/lambda; kws = kw*kw;kr1 = kw*ra1;kr2 = kw*ra2;krs1 = kr1*kr1;krs2 = kr2*kr2;vs1 = krs1*(ns1-ns2);vs2 = krs1*(ns2-ns3);v1 = sqrt(vs1);v2 = sqrt(vs2);m0 = 1000; %計(jì)算總次數(shù) U10 = 10;%U1終值，dU1 = U10/m0;U1 = 1.01*v1;%U1起始值for k = 1:m0 k U1 = U1+dU1;% U1 xu(k) = U1; fcl1 = eclthf

3、(U1); fcl(k) = fcl1;endfigure()plot(xu,fcl)xlabel('U1')ylabel('F')title('三層光纖包層模本征方程')grid on %-%光纖光柵理論基礎(chǔ)與傳感技術(shù) %求解三層光纖包層模本征方程,畫(huà)曲線f(U1),%eclthf:eigenvalue of cladding modes of three layer,function % zhzijia %-function f = eclthf(U1) %eigenvalue of cladding modes of three layer

4、 fiber,function global lambda n1 n2 n3 ra1 ra2 ns1=n1*n1;ns2 = n2*n2;ns3 = n2*n3;s21 = ns2/ns1;s23 = ns2/ns3;ras1 = ra1*ra1;ras2 = ra2*ra2;kw = 2*pi/lambda; %kkws = kw*kw;kr1 = kw*ra1;kr2 = kw*ra2;krs1 = kr1*kr1;krs2 = kr2*kr2;v1 = kr1*sqrt(ns1-ns2);v2 = kr1*sqrt(ns1-ns3);vw = kr2*sqrt(ns2-ns3);vs12

5、 = krs1*(ns1-ns2);vs23 = krs2*(ns2-ns3); ar = ra2/ra1; ars = ar*ar;bus = U1*U1;nes = ns1-bus/krs1;u1 = U1/ra1;u2 = kw*sqrt(ns2-nes);w3 = kw*sqrt(nes-ns3);u11 = U1;u21 = kr1*sqrt(ns2-nes);u22 = kr2*sqrt(ns2-nes); w32 = kr2*sqrt(nes-ns3);bu1 = u1*ra1;bu2 = u2*ra1;bw3 = w3*ra2;bus1 = bu1*bu1;bus2 = krs

6、1*(ns2-nes);bws3 = krs2*(nes-ns3); sigma2 = nes;sigma = sqrt(sigma2);xs1 = (n1*bus1*bus2/(sigma*vs12)2;xs2 = ns3*(ars*bus2*bws3/(sigma*vs23)2; j011 = besselj(0,u11);j111 = besselj(1,u11);j211 = besselj(2,u11); j021 = besselj(0,u21);j221 = besselj(2,u21);j121 = besselj(1,u21); j022 = besselj(0,u22);j

7、122 = besselj(1,u22);j222 = besselj(2,u22); y021 = bessely(0,u21);y221 = bessely(2,u21);y121 = bessely(1,u21); y022 = bessely(0,u22);y122 = bessely(1,u22);y222 = bessely(2,u22); k032 = besselk(0,w32);k132 = besselk(1,w32);k232 = besselk(2,w32); j11 = 0.5*(j011-j211);% j1'j21 = 0.5*(j021-j221);j2

8、2 = 0.5*(j022-j222);y21 = 0.5*(y021-y221);y22 = 0.5*(y022-y222);%k22 = -0.5*(k032+k232); %k'%kb = k22/(w3*k132); % k'/wkjb = u1*j111; %ujk3 = -0.5*(k032+k232)/(bw3*k132); % k'/wk p1 = j122*y121-j121*y122;q1 = j122*y21-j21*y122;r1 = j22*y121-j121*y22;s1 = j22*y21-j21*y22; ju12 = (bu1/bu2)

9、*j111;jus = bus1*j111*j111;jp = j11*p1;kp = k3*p1;kq = k3*q1;au = 1/(ar*bu2);aus = ars*bus2*bus2;rau = r1*au;x12 = 8/(pi*pi)*bus1*bws3*ns2/(sigma2*vs12*vs23);f01 = (p1*p1+x12)*jus;f031 = xs1*xs2*(j11*(kp+rau)-ju12*(kq+au*s1);f032 = (j11*(kp+s23*rau)-s21*ju12*(kq+s23*au*s1);f03 = f031*f032;f04 = -xs1

10、*(jp-q1*ju12)*(jp-s21*q1*ju12);f05 = -xs2*jus*(kp+rau)*(kp+s23*rau);f = f01+f03+f04+f05; %-% 求解三層光纖包層模本征方程,包層模傳播常數(shù)% eclthb:eigenvalue of cladding modes of three layer,beta % 光纖光柵理論基礎(chǔ)與傳感技術(shù) % zhzijia %-clear all,clcformat longglobal lambda n1 n2 n3 ra1 ra2 %更改方便,引入全局量n1 = 1.4681;n2 = 1.4628;n3 = 1;ns1

11、 = n1*n1;ns2 = n2*n2;ns3 = n2*n3; ra1 = 4.15e-6;ra2 = 62.5e-6; lambda = 1550e-9; kw = 2*pi/lambda; kws = kw*kw;kr1 = kw*ra1;kr2 = kw*ra2;krs1 = kr1*kr1;krs2 = kr2*kr2;vs1 = krs1*(ns1-ns2);vs2 = krs1*(ns2-ns3);v1 = sqrt(vs1);v2 = sqrt(vs2);m0 = 1000; %計(jì)算總次數(shù) U10 = 7;%U1終值，dU1 = U10/m0;U1 = 1.01*v1;%U1

12、起始值bum = U1;m = 0;m1 = 50;% 求出beta的最大數(shù)目for k = 1:m0 bum = bum+dU1; xu(k) = bum; ya = eclthf(bum); bua = bum; y(k) = ya; yaa = ya; if k = 1 yab = ya; end if yaa*yab<0 m = m+1 beta(m) = eclthbf(bua,bub); if m>= m1 break end end bub = bum; yab = ya;end plot(xu,y)xlabel('U_1')ylabel('f

13、 (U_1)')grid on beta' %-% 求解三層光纖包層模本征方程,包層模傳播常數(shù),判斷零點(diǎn)，% eclthbf:eigenvalue of cladding modes of three layer,beta ,function% 光纖光柵理論基礎(chǔ)與傳感技術(shù) % zhzijia %-function beta = eclthbf(bua,bub) global lambda n1 n2 n3 ra1 ra2 ns1 = n1*n1;ns2 = n2*n2;ns3 = n2*n3;kw = 2*pi/lambda; %k kws = kw*kw;kr1 = kw*ra1;kr2 = kw*ra2;krs1 = kr1*kr1;krs2 = kr2*kr2;ee = 1e-10;c = 1;m1 = 0;while c m1 = m1+1; if m1 = 100 break end bu2 = 0.5*(bua+bub); z = eclthf(bu2); z1 =