matlab代码：
clcclear allclose all%%N = 200;lambda = 1064e-6;%波长1064nmrow = 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:Nfor b=1:NE2(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))));endaendI2 = 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]. 光子学报, , 34(11):1732-1735.
[2]梁晓晶. 贝塞尔光束的传输变换特性研究[D].
clc; clear; close all;%% 参数设置Nxy = 200; %采样点lambda = 1064e-9;%波长1064nmk = 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:9theta = 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
