Распространение оптического поля методом углового спектра и проверка ⇐ Python

[Anonymous]

Я написал код для распространения гауссова луча в свободном пространстве методом углового спектра.
Чтобы проверить правильность моделирования или нет, я распространяю гауссов луч в свободном пространстве и проверяю результат ширина луча поля по теоретическим формулам. Ниже приведен мой код, я не понимаю, что в нем не так, ширина луча не соответствует теоретической ширине луча, я был бы очень признателен за любые рекомендации. Спасибо.
Вот мой код:

Код: Выделить всё

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)

# With these parameters, the_w (theoretical beam width) = 0.0707; but ana_w (analytical beam width) = 0.0792

Подробнее здесь: https://stackoverflow.com/questions/793 ... rification