角谱法高斯光束自由空间传播仿真的束腰宽度验证问题求助

大家好，我用角谱法写了一段自由空间高斯光束传播的仿真代码，现在想通过对比理论束腰宽度来验证仿真是否正确，但结果和理论值对不上，实在找不到问题出在哪，希望能得到大家的指导，谢谢！

以下是我的代码：

import numpy as np

wvl   = 1550e-9             # wavelength
wz    = 0.05                # beam waist
gdim  = 1                   # spatial extent of the grid;
resol = 512                 # grid resolution
I     = 1                   # intensity amplitude
zr    = (np.pi * wz**2)/wvl # rayleigh range
k     = 2*np.pi/wvl

# Theoretical beam width at reciever (Lz km) # Gaussian 
wzt     = lambda Lz, wg, lmda: wg * np.sqrt(1 + ((Lz * lmda)/(np.pi * wg**2))**2)

wz_     = wzt(0, wz, wvl)
z       = zr

# define grid
dx      =  np.sqrt((wvl*z)/resol) #(wvl * np.sqrt(z**2 + gdim**2))/gdim
gdim    =  dx * resol
x, y    =  np.meshgrid(np.arange(-gdim/2, gdim/2, dx),
                     np.arange(-gdim/2, gdim/2, dx))

r    = np.sqrt(x**2 + y**2)

def gaussian(wz, r, I):
    """
    wz: Beam width at z=0;
    r : Radial coordinates
    I : Intensity distribution
    """
    Fin = I*np.exp(-2*(r/wz)**2)
    return Fin
gbeam = gaussian(wz, r, I)
def ft2(g, delta=None):
    """
    Computes the 2D Fourier transform of g with scaling.
    
    Parameters:
    g : numpy.ndarray
        Input 2D array to transform.
    delta : float
        Spacing in the spatial domain.
    
    Returns:
    numpy.ndarray
        The scaled 2D Fourier transform of g.
    """
    G = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(g)))
    return G

def ang_spec_prop(e, z, resol, dx, wvl):
    k = (2*np.pi/wvl)
    
    fxx = np.fft.fftshift(np.fft.fftfreq(resol, dx))
    fxx, fyy = np.meshgrid(fxx, fxx)
    
    alfa = k ** 2 - (4 * (np.pi ** 2)) * (fxx ** 2 + fyy ** 2)
    
    tmp = np.sqrt(np.abs(alfa))
    kz  = np.where(alfa >= 0, tmp, 1j*tmp)
     
    h1   = np.exp(z * kz * 1j)
    e_ = ft2(ft2(e)*h1)    
    return e_
prop_as = ang_spec_prop(gbeam, z, gbeam.shape[0], dx, wvl)

# Calculate the beam width
def anacal_gbw(ub, r):
    # Find the intensity where it falls to 1/e^2 of its maximum
    ub       = np.abs(ub)**2
    r_flat   = r.flatten()
    imax     = ub.max()/(np.exp(1)**2)
    in_flat  = ub.flatten()
    wg_calc  = r_flat[np.abs(in_flat-imax).argmin()]
    return wg_calc
ana_w = anacal_gbw(prop_as, r) 
the_w = wzt(z, wz, wvl)

用这些参数计算时，理论束腰宽度the_w为0.0707，但仿真后得到的束腰宽度ana_w是0.0792，两者不匹配。

备注：内容来源于stack exchange，提问作者harshit tiwari