Matlab：贝塞尔高斯光束自由传输matlab仿真

matlab代码：

clc
clear all
close all
%%
N = 200;
lambda = 1064e-6;               %波长1064nm
row = linspace(-1,1,N); col = linspace(-1,1,N);
[x,y] = meshgrid(row,col);
[theta,r] = cart2pol(x,y);
w = 3;                          %高斯光束束腰宽度
k = 2*pi/lambda;                %波数
k_r = 20;                       %径向波矢 - 常量
k_z = sqrt(k^2-k_r^2);          %轴向波矢
z = 0;
for n = 0 : 3                    %贝塞尔函数阶数n = 0,1,2,3等等
    E = besselj(n,k_r*r).*exp(-r.^2/w^2).*exp(1i*k_z*z).*exp(1i*n*theta);
    I = E.*conj(E);
    I = I/max(max(I));          %归一化
    figure;mesh(x,y,I)
    set(gca,'fontname','times new roman','fontsize',16);
    title([num2str(n),'阶贝塞尔-高斯光束光强分布'],'fontname','华文中宋','fontsize',16);
    xlabel('x/mm','fontname','times new roman','fontsize',16);
    ylabel('y/mm','fontname','times new roman','fontsize',16);
    zlabel('归一化强度','fontname','华文中宋','fontsize',16);
end
%% 自由传输部分(菲涅尔衍射积分)
n = 1;
E1 = besselj(n,k_r*r).*exp(-r.^2/w^2).*exp(1i*k_z*z).*exp(1i*n*theta);
Z = 100;   %传输距离
x00 = linspace(-0.5,0.5,N); y00 = linspace(-0.5,0.5,N);
[x2,y2] = meshgrid(x00,y00);
for a=1:N
    for b=1:N
        E2(a,b) = -1i/lambda/Z*exp(1i*k*Z)*sum(sum(E1.*exp(1i*k/2/Z.*((x00(a)-x).^2+(y00(b)-y).^2))));
    end
    a
end
I2 = E2.*conj(E2);      I2 = I2/max(max(I2));
figure;mesh(x2,y2,I2)
set(gca,'fontname','times new roman','fontsize',16);
title([num2str(n),'阶贝塞尔-高斯光束自由传输后光强分布'],'fontname','华文中宋','fontsize',16);
xlabel('x/mm','fontname','times new roman','fontsize',16);
ylabel('y/mm','fontname','times new roman','fontsize',16);
zlabel('归一化强度','fontname','华文中宋','fontsize',16);

结果如下：





这里用的是菲涅尔衍射积分来计算传出100mm后的光场，可以改变Z的大小来计算传输更远距离传输距离的光场，不过如果传输距离过大，菲涅尔衍射积分公式将不再适用，请慎重使用。
参考文献：
[1]Gao Zenghui, Lü Baida. 非傍轴贝塞尔-高斯光束在自由空间的传输[J]. 光子学报, 2005, 34(11):1732-1735.
[2]梁晓晶. 贝塞尔光束的传输变换特性研究[D].

clc;    clear;  close all;
%% 参数设置
Nxy = 200;                      %采样点
lambda = 1064e-9;               %波长1064nm
k = 2*pi/lambda;                %波数
w = 3e-3;                       %高斯光束束腰宽度
[x,y] = meshgrid(linspace(-1.5*w,1.5*w,Nxy));
[theta,r] = cart2pol(x,y);
k_r = 7000;                       %径向波矢 - 常量
k_z = sqrt(k^2-k_r^2);          %轴向波矢
z = 0;
c = 3e8;
T = lambda/c;   %周期
m = 0;          %拓扑荷数 
omega = 2*pi*c/lambda;
for jt = 1:9
    theta = imrotate(theta, 90*(jt-1));      %每个周期相位旋转90° 本来是想旋转任意角度的 可惜相位旋转之后矩阵维度不匹配了
    E = besselj(m,k_r*r).*exp(-r.^2/w^2).*exp(1i*k_z*z).*exp(1i*m*theta).*exp(1i*omega*T*(jt-1));
    I = E.*conj(E);     I = I/max(max(I));
    m = m + 1;          %因为涡旋相位是旋转对称的 相位旋转也看不出来 所以改变拓扑荷数看着更直观一些
    figure;mesh(x*1e3,y*1e3,I);
    view(2);
    set(gca,'fontname','times new roman','fontsize',16);
    title([num2str(jt-1),'T时的光强分布'],'fontname','华文中宋','fontsize',16);
    xlabel('x/mm','fontname','times new roman','fontsize',16);
    ylabel('y/mm','fontname','times new roman','fontsize',16);
    zlabel('归一化强度','fontname','华文中宋','fontsize',16);
    %set( gca, 'XTick', [], 'YTick', [] );
    F1 = getframe;
    imwrite(F1.cdata,[ 'C:\Users\Administrator\Desktop\贝塞尔\t = ',num2str(jt-1),'T.png']);  %保存图片
end