【光波电子学】MATLAB绘制光纤中线性偏振模式LP之多模光纤的电场分布(光斑)

单模光纤的分析参考本人另一篇博客
【光波电子学】MATLAB绘制光纤中线性偏振模式LP之单模光纤的电场分布(光斑)

详细的理论见另一篇文章
【光波电子学】基于MATLAB的多模光纤模场分布的仿真分析

1 原理

光纤中，在弱导近似下，由$H{E}_{l+1,m}$模和$E{H}_{l-1,m}$模组合的模式是一个线性偏振模式英文是(linear polorized mode)简称为$L{P}_{lm}$。根据资料以下是$L{P}_{lm}$模和导出或构成它的模式、简并度及模式特征方程。

采用高斯模型绘制光纤中的模场强度或电场强度。光纤中光场模式行为是亥姆霍兹方程解的形式，因此模场的精确解在纤芯中用第一类Bessel函数描述，在包层中用第二类变型Bessek函数描述。
根据表3.3 ，可以得到$L{P}_{1m}$模是由于一个$T{E}_{0m}$模、一个$T{M}_{0m}$模以及两个一共4个模式简并而成，其特征方程是
$$\frac{J_1(U)}{U{J}_0(U)}=-\frac{K_1(W)}{W{K}_0(W)}$$

2 Matlab 实现绘制多模光纤模场分布

2.1 第一步

用MATLAB对该特征方程数值求解

clc
clear
close all
tic % 程序运行开始计时
Vmax = 10; %V的最大值为10
N =100; % 共计算(0:0.1:10)中的100个点

for j = 1:N
    V(j) = j/N*Vmax;
    Vtemp = V(j);
    Utemp = NaN;
    i = 0;
    while(isnan(Utemp)&&i<N+1)
%         init = 3.5*(N-1)/N;
        init = Vmax*(N-i)/N;
        try
            Utemp =fzero(@(Utemp)...%求解特征方程的函数
                besselj(1,Utemp)/(Utemp*besselj(0,Utemp))+...
                besselk(1,sqrt(Vtemp^2-Utemp^2))/(sqrt(Vtemp^2-Utemp^2)*...
                besselk(0,sqrt(Vtemp^2-Utemp^2))),init);
        catch
        end
        i = i+1;
    end
    U(j) = Utemp;
end
W = sqrt(V.^2-U.^2);
Ymax = ceil(max([U,W]));

figure
subplot(1,2,1)
plot(V,U)
axis equal
axis([0 Vmax 0 Ymax])
xlabel('V')
ylabel('U')
title('LP_{11}V-U')

subplot(1,2,2)
plot(V,W);
axis equal
axis([0 Vmax 0 Ymax])
xlabel('V')
ylabel('W')
title('LP_{11}V-W')
toc %程序运行结束计时

该图中V-U曲线有三段，分别对应$L{P}_{10}$模、$L{P}_{11}$模和$L{P}_{12}$模的V-U 曲线
根据上面求解得到当V=5时，U=3.1527，W=3.8808

2.2 第二步

绘制二维光场分布

= 5 ;
U =3.1527;
W = 3.8808;

Npoint = 501;
R1 = linspace(-5,5,Npoint);
R2 = linspace(-5,5,Npoint);

X= meshgrid(R1,R2);
Y= meshgrid(R2,R1);
Y = Y';
R = sqrt(X.^2+Y.^2);
Theta = atan(Y./(X+eps));

% 光纤中光芯光场分布
E1 = besselj(1,U*R).*cos(Theta);
I1 = E1.^2;

% 光纤中包层的光场分布
E2 = besselj(1,U).*besselk(1,W.*R).*cos(Theta)/besselk(1,W);
I2 = E2.^2;
I = I1 ;
pos = find(R>=1);
I(pos) = I2(pos);
I  = I/max(max(I1));%对光场分布进行归一化处理

% 作图
imagesc(R1,R2,I,[0 1]);
colormap(gray);
colorbar
xlabel('x')
ylabel('y')
zlabel('z')

2.3 第三步

绘制三维光场分布

V= 5
U =3.1527;
W = 3.8808;

Npoint = 201;
R1 = linspace(0,1,Npoint);
R2 = linspace(1,5,Npoint);
Theta1 = linspace(0,2*pi,Npoint);
Theta2 = linspace(0,2*pi,Npoint);

E1 = zeros(Npoint,Npoint);
E2 = zeros(Npoint,Npoint);
I1 = E1;
I2 = E2;
% 光纤中光芯光场分布
for i = 1:Npoint
   for j=1:Npoint
        E1(i,j) = besselj(1,U*R1(i))*cos(Theta1(j));
        I1(i,j) = E1(i,j).^2;
   end
end
% 光纤中包层的光场分布
for i = 1:Npoint
    for j = 1:Npoint
        E2(i,j) = besselj(1,U).*besselk(1,W.*R2(i))./besselk(1,W)*cos(Theta1(j));
        I2(i,j) = E2(i,j).^2;
    end
end
% 极坐标转换为柱坐标，以便于作图
[Theta1 R1] = meshgrid(Theta1,R1);
[Theta2 R2] = meshgrid(Theta2,R2);
[X1 Y1] = pol2cart(Theta1,R1);
[X2 Y2] = pol2cart(Theta2,R2);

% 三维作图
mesh(X1,Y1,I1)
hold on
mesh(X2,Y2,I2)
colorbar
xlabel('x')
ylabel('y')
zlabel('z')
title(['LP_{11}Mode,V = 5'])