作者B站ID“大洋咩咩咩咩”,投稿有对该代码的讲解视频(角谱传播),如果有那个地方有疑问可以去视频中寻找答案~
文章及代码解决三个问题:
夫琅禾费衍射(理想凸透镜成像)做FFT2后,仿真光斑的实际尺寸确定?
光场在自由空间传播一定距离z后的光斑什么样,怎么用FFT2做?
如何仿真泽尼克多项式光斑,例如慧差光斑,以及实际尺寸?
clc,clear
N_obj_sam=512; %采样频率
D=2.5;%% 圆孔径mm
plotW=4;%% 绘图宽度mm
z=600; %% 衍射距离 mm
lambda=632.8e-6; %% 波长mm
N=6;%% 补零数
%============================================
%物光场
plotW2=plotW/2;
R=D/2;
Obj_x=linspace(-plotW2,plotW2,N_obj_sam);
[Obj_x,Obj_y]=meshgrid(Obj_x,Obj_x);
Obj_A2=zeros(N_obj_sam);Obj_P=Obj_A2;
Obj_A2((Obj_x.^2+Obj_y.^2)<=R^2)=1;
Obj_E=Obj_A2.*exp(1i*Obj_P);
%============================================
%物光场补零
Obj_Ezf=F_ZerosFilling(Obj_E,N_obj_sam,N*N_obj_sam);
A=fftshift(fft2(fftshift(Obj_Ezf)));
[m,n]=size(A);%% 角谱的规格
dtheta=lambda/(N*plotW);%% 角谱角分辨率
A_x=(-m/2:m/2-1)*dtheta;A_y=(-n/2:n/2-1)*dtheta;% 角谱的角分布
[A_x,A_y]=meshgrid(A_x,A_y);%% 构造角谱
%角谱传播
k=2*pi/lambda;%波矢量
H=exp(1i*k*z*sqrt(1-(A_x.^2+A_y.^2)));%自由空间光学传递函数
B=A.*H;
%逆变换成像传播
Img_Ezf=fftshift(ifft2(fftshift(B)));
Img_E=F_UnZerosFilling(Img_Ezf,N_obj_sam*2);
subplot(231)
F_2DPlot(Obj_A2,plotW2,'mm','mm',strcat('物透过率 R=',num2str(R),'mm'))
subplot(232)
F_2DPlot(Obj_P,plotW2,'mm','mm',strcat('物相位 \lambda=',num2str(lambda*10^6),'nm'))
subplot(233)
x_plot_obj=linspace(-plotW2,plotW2,N_obj_sam);
plot(x_plot_obj,Obj_A2(N_obj_sam/2,:),'color',[0 0 0],'LineWidth',3)
xlabel('mm');ylabel('I');xlim([-plotW2 plotW2]);
title(strcat('通光孔径 D=',num2str(D),'mm'))
subplot(234)
Img_A2=F_NorOne(abs(Img_E).^2);
F_2DPlot(Img_A2,plotW,'mm','mm',strcat('像面光斑 z=',num2str(z),'mm'))
subplot(235)
F_2DPlot(angle(Img_E),plotW,'mm','mm','像面相位')
subplot(236)
x_plot_img=linspace(-plotW,plotW,N_obj_sam*2);
plot(x_plot_img,Img_A2(N_obj_sam,:),'color',[0 0 0],'LineWidth',3)
xlabel('mm');ylabel('I');xlim([-plotW plotW]);
title(strcat('艾里斑半径理论值 R_A_i_r_y= 1.22\lambdaz/D=',num2str(1.22*lambda*z/D),'mm'))
function F_2DPlot(z,k,xl,yl,tl)
%本函数目的:绘制3维曲面图,去边框,带标题和label
%输入变量:k为放缩坐标轴的因子
image([-1*k 1*k],[-1*k 1*k],z,'CDataMapping','scaled')
xlabel(xl);
ylabel(yl);
title(tl);
axis xy
axis equal
xlim([-1*k 1*k]);
ylim([-1*k 1*k]);
view([0 90]);
end
function MFill=F_ZerosFilling(MOnes,aOnes_F,NumPixcel)
%子函数序号#7
MFill=zeros(NumPixcel);
[aOnes_T,~]=size(MOnes);%实际输入的波前函数,可以是366或367
StartP=(NumPixcel-aOnes_F)/2;%获取规格差距
MFill(StartP+1:StartP+aOnes_T,StartP+1:StartP+aOnes_T)=MOnes;
end
function O_A2_Plot=F_UnZerosFilling(O_A2,Samp_p)
[a,b]=size(O_A2);
O_A2_Plot=O_A2(a/2-Samp_p/2+1:a/2+Samp_p/2,b/2-Samp_p/2+1:b/2+Samp_p/2);
end
function A_n=F_NorOne(A)
mA=min(A(:));
A=A-mA;
MA=max(A(:));
A_n=A/MA;
end