Bessel光束经透镜的matlab模拟传输_光束传输方法的matlab模拟

哇咔咔，试着写了一个贝塞尔光束经过透镜传输的matlab模拟代码。

看了篇文献，感觉贝塞尔光束经过透镜蛮神奇的，就按照里面的参数模拟了一下，结果还是蛮一致的。

文献：梁晓晶. 贝塞尔光束的传输变换特性研究[D].天津大学,2016.

用了一个网格划分+ABCD系统的科林斯积分公式捏。

为了运算速度，使用了矩阵算法，大大加速了速度。参考了一下这个文献（matlab矩阵库yyds）

龚海龙, 李国俊, 陈琪, et al. 一种快速实现光场传输的算法研究[J]. 光学学报, 2016

然后顺便写了一节用于跑出视频并保存的代码，去掉注释即可食用（别忘了去掉上面一节的getframe那一句哦~）龚海龙

ps：当然，如果只是需要几个点的截图只要把循环去掉即可食用~

%% 贝塞尔光束 经ABCD系统传输
clc
clear all
close all
%%
N = 2000; %光源处取格点数
N1=2000;  %观测处取格点数
lambda = 663e-6;               %波长663nm
front=2;%生成区间
hind=2;%观测取间
Z =100;           %传输距离z mm
f=50;  %透镜焦距（mm）
A=1-Z/f; 
B=Z;
C=-1/f;
D=1;%矩阵光学总矩阵
 
hindx = linspace(-front,front,N); hindy = linspace(-front,front,N);
[x1,y1] = meshgrid(hindx,hindy);
[theta,r] = cart2pol(x1,y1);
k = 2*pi/lambda;                %波矢
k_r = 119.175;                       %径向波矢 
k_z = sqrt(k^2-k_r^2);          %轴向波矢
z = 0;
w=1;%束腰
n =  1;                  %贝塞尔函数阶数n = 0,1,2,3。。。
E1 = besselj(n,k_r*r).*exp(1i*k_z*z).*exp(1i*n*theta);%.*exp(-r.*r ./(w^2)) ;% 就问这束腰你要不要吧
    I = E1.*conj(E1);
    I = I/max(max(I));          %归一化
    figure
    plot(x1(1,:),I(N/2,:));
    title([num2str(n),'阶贝塞尔光束光强分布'],'fontname','华文中宋','fontsize',16);
    xlabel('x/mm','fontname','times new roman','fontsize',16);
    zlabel('归一化强度','fontname','华文中宋','fontsize',16);
    
    figure;mesh(x1,y1,I)
    set(gca,'fontname','times new roman','fontsize',16);
    colorbar;
    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);
    view(0,90) %调整立体图的观察角度
%% 矩阵 启动！
figure
for i=1:90
Z =i*0.5+40;           %传输距离z mm
f=50;
A=1-Z/f; 
B=Z;
C=-1/f;
D=1;%矩阵光学总矩阵
    gedianx = linspace(-hind,hind,N1); gediany = linspace(-hind,hind,N1);
    [x2,y2] = meshgrid(gedianx,gediany);
    M1=E1.*exp(1i*k*A/2/B *(x1.^2+y1.^2));
    M2=exp(1i*k*D/2/B *(x2.^2+y2.^2));
    Mx=exp(-1i*k/B *x2(1,:)'*x1(1,:));
    My=exp(-1i*k/B *y1(:,1)*y2(:,1)');
    E2=-1i/lambda/B *exp(1i*k*Z).*M2.*(Mx*M1*My); %.*(hind/N1).*(hind/N1)  咱不归一就得写这个;
    I1 = E2.*conj(E2);   I1 = I1/max(max(I1));
 %这个是画立体图
%     mesh(x2,y2,I1)
%      colorbar;
%     set(gca,'fontname','times new roman','fontsize',16);
%     title([num2str(n),'阶贝塞尔光束经透镜传输',num2str(Z),'mm后光强分布'],'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);
%     view(0,90)
%     
 
%这个是画截面图
    plot(y2(:,1),I1(:,N1/2));
    title([num2str(n),'阶贝塞尔光束经透镜传输',num2str(Z),'mm后光强分布'],'fontname','华文中宋','fontsize',16);
    xlabel('y/mm','fontname','times new roman','fontsize',16);
    zlabel('归一化强度','fontname','华文中宋','fontsize',16);
    drawnow; 
 %   pic1(i)=getframe(gcf);%获得当前帧（得到视频）
end
%% 出视频
% now=VideoWriter("test1.avi"); %生成avi视频的名字
% now.FrameRate=15; %帧数
% open(now);
%  writeVideo(now,pic1);
%  close(now);

 跑完差不多是这么个样子：

 和文献上的图还是蛮吻合的捏就是说。