使用python的角谱方法

Question

我正在尝试使用角谱方法在数值上传播给定的（电场）。为此，我遵循“傅立叶光学的原理和应用”（罗伯特·泰森），第3章，第2页

我尝试使用以下代码重新创建数学

import numpy as np
import imageio

U = imageio.imread("ap.png")[:,:, 0] # load circular aperture
_lambda = 800e-9

def propagate2(self,z):
    A = np.fft.fft2(U, norm="ortho") # 2D FFT 
    alpha = np.fft.fftfreq(U.shape[0])*_lambda # direction cosine in x direction
    beta = np.fft.fftfreq(U.shape[1])*_lambda # direction cosine in x direction
    gamma = np.zeros([alpha.shape[0], beta.shape[0]])
    k = 2*np.pi/_lambda # wavevector

    for i,j in itertools.product(range(alpha.shape[0]), range(beta.shape[0])): # determine phase value for each (i,j)
        if alpha[i]**2+beta[j]**2 < 1:
            gamma[i,j] = np.sqrt(1-alpha[i]**2-beta[j]**2)
        else:
            gamma[i,j] = 1j*np.sqrt(np.abs(1-alpha[i]**2-beta[j]**2))
    phi = np.exp(1j*k*z*gamma)
    field = np.fft.ifft2(A*phi, norm="ortho") # 2D IFFT
    return field

此代码应该产生通常的双缝衍射图案，但是（如下所示）根本不会产生衍射。

我相当确定我的alpha和beta值存在问题，但是我似乎找不到它。非常感谢您的帮助。

ap.png：

Answer 1

要正确解决这个问题可能很棘手。这是：

u = ... # this is your 2D complex field that you wish to propagate
z = ... # this is the distance by which you wish to propagate

dx, dy = 1, 1 # or whatever

wavelen = 1 # or whatever
wavenum = 2 * np.pi / wavelen

kx = np.fft.fftfreq(u.shape[0], dx / (2 * np.pi))
ky = np.fft.fftfreq(u.shape[1], dy / (2 * np.pi))

# this is just for broadcasting, the "indexing" argument prevents NumPy
# from doing a transpose (default meshgrid behaviour)
kx, ky = np.meshgrid(kx, ky, indexing = 'ij', sparse = True)

kz_sq = kx * kx + ky * ky

# we zero out anything for which this is not true, see many textbooks on
# optics for an explanation
mask = wavenum * wavenum > kz_sq

g = np.zeros((len(kx), len(ky)), dtype = np.complex_)
g[mask] = np.exp(1j * np.sqrt(wavenum_sq - kz_sq[mask]) * z)

res = np.fft.ifft2(g * np.fft.fft2(u)) # this is the result

您可能想要垫住u以防止回绕。在这种情况下，只需计算形状加倍的g，然后将结果切成薄片。