光纤模式分布 matlab,matlab计算单模光纤模式分布(公布源代码及参考文献)

最近在使用matlab计算单模光纤纤芯模及包层模模场分布时，有一些问题一直悬而未决，多次咨询原作者后虽解决了部分问题，但是余下的问题原作者也不理我了，特发此贴以广交学习光纤方面的同学、老师及科研人员，希望大家互相学习，互相交流，共同探讨，共同进步。

首先，求解单模光纤纤芯模特征方程，光纤为Corning SMF-28，具体参数见以下代码，参考文献为：

% [1]. Optical Fibre Waveguide Analysis_Charles Tsao_Oxford 1992 ----% %

%     Page 315--------------% %

% [2]. Eigenvalue and field equations of three-layered uniaxial fibers

% and their applications to the characteristics of long-period fiber

% gratings with applied axial strain ---------% %

% Zi-jia Zhang,J.Opt.Soc.Am. A/Vol.22,No,11/November 2005 --------%

% -------------------------------------------------------------------------

首先光纤模型采用的是三层模型，相应代码如下ecothf_HE.m：(用过matlab计算光纤模场分布的都应该能看懂下面的代码)

% -------------------------------------------------------------------------

% The Eigenvalue Equation of the Core Hybrid HE or EH mode for

% a general Three-layered Single mode Fiber  ------------% %

% % --------------------------------Reference---------------------------% %

% [1]. Optical Fibre Waveguide Analysis_Charles Tsao_Oxford 1992 ----% %

%     Page 315--------------% %

% [2]. Eigenvalue and field equations of three-layered uniaxial fibers

% and their applications to the characteristics of long-period fiber

% gratings with applied axial strain ---------% %

% Zi-jia Zhang,J.Opt.Soc.Am. A/Vol.22,No,11/November 2005 --------%

% -------------------------------------------------------------------------

function f = ecothf_HE(U1)

global lambda n1 n2 n3 ra1 ra2

ns1 = n1*n1;

ns2 = n2*n2;

ns3 = n3*n3;

s21 = ns2/ns1;

s23 = ns2/ns3;

kw = 2*pi/lambda;             %  kw = k0

kr1 = kw*ra1;

krs1 = kr1*kr1;

kr2 = kw*ra2;

krs2 = kr2*kr2;

bvs12 = krs1*(ns1-ns2);         % V12^2 = k0^2*ra1^2*(n1^2-n2^2)

bvs23 = krs2*(ns2-ns3);

bus1 = U1*U1;                  % U1^2

nes = ns1-bus1/krs1;           % neff^2 = n1^2-U1^2/(k0*ra1)^2

ne = sqrt(nes);

u1 = U1/ra1;

w2 = kw*sqrt(nes-ns2);        % w2 = k0*sqrt(n2^2-neff^2)

w3 = kw*sqrt(nes-ns3);        % w3 = k0*sqrt(neff^2-n3^2)

bu1 = u1*ra1;                 % U1 = u1*ra1

bw2 = w2*ra1;                 % W2 = w2*ra1

bws2 = bw2*bw2;

bw3 = w3*ra2;                 % W3 = w3*ra2

bws3 = bw3*bw3;               % W3^2 = k0^2*ra2^2*(neff^2-n3^2)

u11 = bu1;

w21 = w2*ra1;

w22 = w2*ra2;

w32 = w3*ra2;

ar = ra2/ra1;                 % ar = alpha2

ars = ar*ar;

%%

j011 = besselj(0,u11);         % j011 = besselj(0,u1*a1)

j111 = besselj(1,u11);         % j111 = besselj(1,u1*a1)

j211 = besselj(2,u11);         % j211 = besselj(2,u1*a1)

i021 = besseli(0,w21);         % i021 = besseli(0,w2*a1)

i121 = besseli(1,w21);         % i121 = besseli(1,w2*a1)

i221 = besseli(2,w21);         % i221 = besseli(2,w2*a1)

i022 = besseli(0,w22);         % i022 = besseli(0,w2*a2)

i122 = besseli(1,w22);         % i122 = besseli(1,w2*a2)

i222 = besseli(2,w22);         % i222 = besseli(2,w2*a2)

k021 = besselk(0,w21);         % k021 = besselk(0,w2*a1)

k121 = besselk(1,w21);         % k121 = besselk(1,w2*a1)

k221 = besselk(2,w21);         % k221 = besselk(2,w2*a1)

k022 = besselk(0,w22);         % k022 = besselk(0,w2*a2)

k122 = besselk(1,w22);         % k122 = besselk(1,w2*a2)

k222 = besselk(2,w22);         % k222 = besselk(2,w2*a2)

k032 = besselk(0,w32);         % k032 = besselk(0,w3*a2)

k132 = besselk(1,w32);         % k132 = besselk(1,w3*a2)

k232 = besselk(2,w32);         % k232 = besselk(2,w3*a2)

% besselj'(n,x) = 0.5*(besselj(n-1,x)-besselj(n+1,x))

% besseli'(n,x) = 0.5*(besseli(n-1,x)+besseli(n+1,x))

j11 = 0.5*(j011-j211);         % besselj'(1,u1*a1) = 0.5*(besselj(0,u1*a1)-besselj(2,u1*a1))

i21 = 0.5*(i021+i221);         % besseli'(1,w2*a1) = 0.5*(besseli(0,w2*a1)+besseli(2,w2*a1))

i22 = 0.5*(i022+i222);         % besseli'(1,w2*a2) = 0.5*(besseli(0,w2*a2)+besseli(2,w2*a2))

k21 = -0.5*(k021+k221);         % besselk'(1,w2*a1) = -0.5*(besselk(0,w2*a1)+besselk(2,w2*a1))

k22 = -0.5*(k022+k222);         % besselk'(1,w2*a2) = -0.5*(besselk(0,w2*a2)+besselk(2,w2*a2))

k32 = -0.5*(k032+k232);         % besselk'(1,w3*a2) = -0.5*(besselk(0,w3*a2)+besselk(2,w3*a2))

p1 = i122*k121-i121*k122;

q1 = i122*k21-i21*k122;

r1 = i22*k121-i121*k22;

s1 = i22*k21-i21*k22;

kb = k32/(bw3*k132);           % Kb = besselk'(1,w3*a2)/(W3*besselk(1,w3*a2))

jus = bus1*j111*j111;          % jus = U1^2*(besselj(2,u1*a1))^2

ju = sqrt(jus);

aw = 1/(ar*bw2);               % au = 1/(alpha2*U2)

aws1 = ars*bws2;               % aus1 = alpha2^2*U2^2

qw1 = q1/bw2;

rw1 = r1*aw;

sw1 = s1*aw;

was2 = 4/(pi*pi*aws1*bws2);

jp = j11*p1;                   % jp = besselj'(1,u1*a1)*p1

kp = kb*p1;                    % kp = k3*p2 = K*p1

kq = kb*q1;                    % kq = k3*q2 = K*q1

psq = jp-s21*qw1*ju;

psr = kp+s23*rw1;

pq = jp-qw1*ju;

pr = kp+rw1;

deltas = nes;                      % deltas = delta^2 = (neff)^2

delta = sqrt(deltas);

x1 = n1*bus1*bws2/(delta*bvs12);

x2 = n3*ars*bws2*bws3/(delta*bvs23);

xs1 = x1*x1;

xs2 = x2*x2;

f1 = (p1*p1-2*ns2*x1*x2/(aws1*bws2))*jus;

f2 = xs1*xs2*(j11*(kp-rw1)+ju*(kq-sw1)/bw2);

f3 = (j11*(kp-s23*rw1)+ju*s21*(kq-s23*sw1)/bw2);

f4 = -xs1*psq*pq;

f5 = -xs2*psr*pr*jus;

f= f1+f2*f3+f4+f5;

另外编写了一个ecothp_HE.m文件求解上面的方程，首先绘函数f的图形，找到f=0点的大致区间，再求根，具体代码如下：

%-------------------------------------------------------------------------

% coremode_HE11:eigenvalue of core modes of three layer fiber ,ploting f

%--------------------------------------------------------------------------

clear all

close all

clc

format long

global lambda n1 n2 n3 ra1 ra2 U_HE11

n1 = 1.4681;

n2 = 1.4628;

n3 = 1;

ra1 = 4.15e-6;

ra2 = 62.5e-6;

lambda = 1.55e-6;

ns1 = n1*n1;

ns2 = n2*n2;

ns3 = n3*n3;

kw = 2*pi/lambda;

kws = kw*kw;

kr1 = kw*ra1;

krs1 = kr1*kr1;

bvs1 = krs1*(ns1-ns2);      % v1^2 = k0^2*ra1^2*(n1^2-n2^2)

bv1 = sqrt(bvs1);

N = 1000;             % the total caculation counts

Umax = bv1;           % For the three-layer SMF,the effective index(ne) of

% the core modes should lie between n1 and n2,

% i.e. n1 > ne > n2.So,U1 lies between bv1 and

% bv2,i.e. bv1 > U1 > 0.

dU = Umax/N;

U = 0;                % initial value of U1

for k = 1:N

U = U+dU;         % U

xu(k) = U;

fcoHE = ecothf_HE(U);

fco(k) = fcoHE;

end

figure()

plot(xu,fco)

xlabel('U')

ylabel('F(U)')

title('Core-mode HE11 of three-layer fiber','fontweight','bold')

grid on

figure()

fplot(@ecothf_HE,[0,1])

title('HE11 mode between [0,1]','fontweight','bold')

xlabel('U')

ylabel('F(U)')

grid on;

U_HE11 = fzero(@ecothf_HE,[0.3,0.4])

U = U_HE11;          % The first root of eigenvalue equation corresponds to the HE11 mode

parameters(U)

a = fzero(@ecothf_HE,[1.4902,1.4904])

parameters(a)

复制代码

以上代码中调用了parameters.m文件，这是求光纤模场其他的一些参数，如neff，beta等，很简单的代码，如下：

function F = parameters(U1)

format long

global lambda n1 n2 n3 ra1 ra2

n1 = 1.4681;

n2 = 1.4628;

n3 = 1;

ns1 = n1*n1;

ns2 = n2*n2;

ns3 = n3*n3;

ra1 = 4.15e-6;

ras1 = ra1*ra1;

ra2 = 62.5e-6;

lambda = 1.55e-6;

kw = 2*pi/lambda;

kws = kw*kw;

kr1 = kw*ra1;

kr2 = kw*ra2;

krs1 = kr1*kr1;

krs2 = kr2*kr2;

vs1 = krs1*(ns1-ns2);      % v1^2 = k0^2*ra1^2*(n1^2-n2^2)

vs2 = krs1*(ns2-ns3);      % v2^2 = k0^2*ra1^2*(n2^2-n3^2)

v1 = sqrt(vs1);

v2 = sqrt(vs2);

bu = U1;

bus = bu*bu;

u1 = bu/ra1

neff = sqrt(ns1-bus/krs1)

beta = kw*neff

运行ecothp_HE.m后，

从函数f的图形中可以看出，f=0有两个点，也就是f=0有两个根，Corning SMF单模光纤应该只有一个纤芯模，即f=0应该只有一个根，该选择哪一个根作为纤芯模的参数呢？？？

然后，假设选择第一个根作为光纤参数，即U_HE11 =  1.561149962979787，再计算纤芯模模场分布归一化系数Eco，参考文献为：(见2楼)

[Last edited by syulzh on 2012-11-22 at 09:00]

以上问题均已解决，若有光纤方面的同仁，请加本人qq交流，共同进步。

[Last edited by syulzh on 2013-9-30 at 19:08]