fourier transform of gaussian beam to simulate far-field beam

Hi,

i have been working on some simulation of optical stuffs
such as the simulation of light beam in free space.

i want to use the fraunhofer approximation method as such
to fourier transform my gaussian beam by using fftshift(fft
(gaussian)). this is simulate the propagated beam at far
field.
f(x,y)= exp(-(x.^2+y.^2)/Wo)

g(x,y)= ho F(x/λd, y/λd) ,where ho=(j/λd)exp(-jkd),
d is the far-field dist,
λ is the wavelength,
g(x,y) is the far-
field beam function,
F is the fourier transform
of f(x,y) which is the
initial beam,
Wo is the initial beam
width.

However i cannot get the correct result in matlab.
The beam-width of the far-field beam seems to be affected
by the spacing of the array. And this makes the simulated
propagated beam different from the expected result. Varying
the spacing of the array(dx and dy), sometimes i will get a
small beam width at far-field. In fact i am using a very
small beam width (eg. 0.001m) and the expected far field
beam waist should be bigger than the initial beam width.

Also, how should i normalize the F(vx,vy)? is it divide by
N^2 ?
By using the Parseual’s Theorem, the power of the initial
and farfield beam should be the same. but in my program,
the power of the initial beam and the farfield beam, they
are not equal.

below is my coding:

N = 2^9;%2^9; %//size of iterations for x and y directions
w_0 = 1.84D-3; %//Beam waist in m
wavelength = 1.064D-6;
dx = 0.05D-3;%//unit length of iterations; per m
dy = dx;
x = [1:N];
y = x;
z =0; % starting position of z
k_0 = 2*pi/wavelength;
dk = pi/N;

Rayleigh= w_0.^2*pi/wavelength

[X,Y] = meshgrid((x-1-N/2),(x-1-N/2));
radius = sqrt((X*dx).^2 + (Y*dx).^2);

u_0 = A.*exp(-(radius/w_0).^2); %//gaussian

P_0= sum(sum(abs(u_0).^2));

u_F= 1/N/N*fft2(u_0);
u_F= fftshift(u_F);

z_d =0.5985; % far-field dist
h_0= j/wavelength/(z_d).*exp(-j*k_0.*z_d);

g_F=h_0.*u_F0;
P_1= sum(sum(abs(g_F).^2));

figure(1)
surf(abs(u_0),'facecolor','flat','edgecolor','none')
camlight left;
lighting phong
figure(2)
surf(abs(g_F),'facecolor','flat','edgecolor','none')
camlight left;
lighting phong


Can anyone help me with this problem ? Thanks...
.