Simulating Gaussian Bessel Beam Propagation in Matlab

Question

I should preface this by saying I am not much of a programmer or physicist, so please bear with me.

I am attempting to simulate the propagation of gaussian Bessel beam (as produced by an axicon) using Fourier propagation (Fresnel propagation?) but I am having some problems.

My code is below, as well as the function I found to do the propagation. The first part of the code does a simple approximation of a gaussian Bessel beam and works fine but lacks some of the characteristic features (diverging into a ring beam, reconstruction following occlusions, etc). For this reason, I wanted to take a cross section of the intensity profile (figure 2) and propagate that forward using the fourier method.

enter image description here

Figure three shows the result of propagation of only ten microns (very short compared to the expected length of the Bessel region (figure 1). For such a distance I would have expected almost no change to the intensity cross section.

I don't know what I am doing wrong, perhaps this method is simply unsuitable to this problem, though it seemed to me like a generic method for beam propagation in free space.

If anyone has some guidance, I would be most appreciative. Thanks.

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

I tried increasing the resolution and/or the side length of the simulation window, on the off chance it was from some interactions with the simulation boundary, but didn't see a qualitative difference for short propagation distances (for large z values it does seem to be interacting with the boundaries, but its already deviated from the expected output long before then).