高斯贝塞尔光束Matlab傅里叶传播模拟异常问题求助

我并非专业程序员或物理学家，尝试用傅里叶传播（菲涅耳传播）方法在Matlab中模拟锥透镜产生的高斯贝塞尔光束，但遇到问题。

初始近似模型能正常运行，但缺少高斯贝塞尔光束的部分特征（如发散为环形光束、遮挡后重建等）。因此我取该模型的强度截面（图2），用傅里叶传播方法向前传播仅10微米（远短于预期的贝塞尔区域长度），结果强度截面出现了明显变化（图3），而我认为如此短距离下应几乎无变化。我尝试提高分辨率或增大模拟窗口，但短距离传播的结果未出现定性改善。不知问题出在哪，是否该傅里叶传播方法不适用于此场景？恳请指导，谢谢。

主代码 #

clear;
clf;

n = 1.44;     %index of refraction for the axicon
a = 5;        %axicon angle alpha
b = asind(n*sind(a))-a;  %the beta angle
bt = 0;       %bessel order
f1 = 75;
f2 = 8;
mag = f2/f1;

z_cone = tand(a)*3;

w = 3;      %guassian beam diameter in mm
y = 1064;     %wavelength in nm
I_0 = 0.1;    %guassian beam intensity
P = 8;      %power
k = 2*3.14/(y*10^(-6)); %wavenumber in mm

lambda = y/(10^6);

z_r = w/2*cosd(b)/sind(b)*mag^2; %bessel length estimate in mm in air
r_0 = 2.4048/(k*sind(b))*mag; % bessel radius in mm in air
Z_0 = 1/2*z_r; %cross section position for xy plot, taken at some fraction of z_r

imax = 512; %1000 good but slow
jmax = 512; %1000 good but slow
spanr = 0.03;%round(20*r_0*1000)/1000;  %mm
spanz = round(20*z_r)/10;  %mm 0.8
I = zeros(imax,jmax);
ii=[0:1:(imax-1)];
spanxy=spanr; %round(100*spanr/sqrt(2))/100;

xy=I; %defining another matrix to hold an xy field perpendicular to the z axis

%k = 2*3.14/(y*10^(-6)); %wavenumber in mm
%z_r = w/(2*(n-1)*tand(a))*mag^2;
%r_0 = 2.4048/(k*sind(b))*mag;

z = (ii)*spanz/imax; %a vector containing all z values
cof1 = (8*P*k*sind(b)/w)*(z/z_r); %parts equation for approximate bessel beam, not dependent on r so taken out of loop
expon1 = exp(-2*z.^2/z_r^2); %another part

for j = 0:(jmax-1);
    r = spanr/2-(j)*spanr/jmax;

    bessel01 = (besselj(bt,r*2.4048/r_0))^2;

    II(j+1,:)=cof1.*expon1.*bessel01;
end

Z_0index = round(Z_0*imax/spanz);
cof1xy = cof1(Z_0index);
expon1xy = expon1(Z_0index);
test = cof1xy*expon1xy;

for u = 0:(imax-1);
  x = -spanxy/2+u*spanxy/imax;
   for v = 0:(jmax-1);
     y = spanxy/2-v*spanxy/imax;
     r = sqrt(x^2+y^2);

     xy((u+1),(v+1)) = test*((besselj(bt,r*2.4048/r_0))^2);
   end
end

%plot 1, approximation
labelsx = [0:spanz/10:spanz];
labelsy = [(-spanr/2)*-1000:(spanr/10)*-1000:(spanr/2)*-1000];
maxI = max(max(II./(1)));
E = sum(sum(II));
subplot(2,2,1), imshow(II./17.6, hot) %why divide by 17.6? just to reduce over saturation. for visual purposes
axis("tic")
set(gca,'xTick',1:imax/10:imax);
set(gca,'yTick',0:jmax/10:jmax);
set(gca,'xTickLabel',labelsx);
set(gca,'yTickLabel',labelsy);
xlabel("Z-axis (mm)");
ylabel('Radial-axis (\mum)')
title('Spatial Intensity')
H=colorbar;
title(H,'Intensity (%)')
set(H, 'yticklabel', [0,20,40,60,80,100])

%plot 2, intensity cross section (seed for propagation)
labelsx = [(-spanxy/2)*1000:(spanxy/10)*1000:(spanxy/2)*1000];
labelsy = [(-spanxy/2)*-1000:(spanxy/10)*-1000:(spanxy/2)*-1000];
subplot(2,2,2), imshow(xy/17.6, hot)
axis("tic")
set(gca,'xTick',1:imax/10:imax);
set(gca,'yTick',0:imax/10:imax);
set(gca,'xTickLabel',labelsx);
set(gca,'yTickLabel',labelsy);
xlabel('X-axis (\mum)');
ylabel('Y-axis (\mum)')
title('Spatial Intensity')
H=colorbar;
title(H,'Intensity (%)')
set(H, 'yticklabel', [0,20,40,60,80,100])

dz = (spanz-Z_0)/(imax-Z_0index);
zz = [Z_0:(spanz-Z_0)/(imax-Z_0index):spanz]; %actual z positions for the propagated region
zz_adjusted = zz(1,:)-Z_0; %adjusted to set z_0 equal to zero for the purpose of the propagation equation
zpmax = imax-Z_0index;


L = spanxy;
[M,~]=size(xy);
dx=L/M;

% Frequency coordinates sampling
%fx=-1/(2*dx):1/L:1/(2*dx)-1/L;

%invoke fourier propagation
u2 = propTF(xy,spanxy,lambda,0.01);

%plot 3, propagated intensity profile
labelsx = [(-spanxy/2)*1000:(spanxy/10)*1000:(spanxy/2)*1000];
labelsy = [(-spanxy/2)*-1000:(spanxy/10)*-1000:(spanxy/2)*-1000];
%subplot(2,2,3), imshow(h(:,:,100)*10,hot) %imshow(U(:,:,10)/17.6, hot)
subplot(2,2,3), imshow(u2/17.6,hot) %imshow(U(:,:,10)/17.6, hot)
axis("tic")
set(gca,'xTick',1:imax/10:imax);
set(gca,'yTick',0:imax/10:imax);
set(gca,'xTickLabel',labelsx);
set(gca,'yTickLabel',labelsy);
xlabel('X-axis (\mum)');
ylabel('Y-axis (\mum)')
title('Spatial Intensity')
H=colorbar;
title(H,'Intensity (%)')
set(H, 'yticklabel', [0,20,40,60,80,100])

传播函数代码 #

function[u2]=propTF(u1,L,lambda,z)
%Fresnel propagation using the Transfer function
% Based on Computational Fourier Optics by Voelz 
%
% Assuming uniform sampling and presents reflections on the boundaries
% 
%PARAMETERS
% 
% u1      - Complex Amplitude of the beam at the source plane
% L       - Sidelength of the simulation window of the source plane
% lambda  - Wavelength 
% z       - Propagation distance 
% u2      - Complex Amplitude of the beam at the observation plane

% Input array size
[M,~]=size(u1); 

% Sampling inteval size
dx=L/M;   

% Frequency coordinates sampling 
fx=-1/(2*dx):1/L:1/(2*dx)-1/L; 

% Momentum/reciprocal Space 
[FX,FY]=meshgrid(fx,fx);

% Transfer function
H=exp(-1j*pi*lambda*z*(FX.^2+FY.^2)); 

H=fftshift(H);

% Fourier transform of the transfer fucntion
U1=fft2(fftshift(u1));

% Convolution of the system
U2=H.*U1;

% Fourier transform of the convolution to the observation plane
u2=ifftshift(ifft2(U2));  
end

内容的提问来源于stack exchange，提问作者Bishop