使用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

我使用 Python 完成了 Angular Spectrum Method 的完整实现： https://github.com/rafael-fuente/Diffraction-Simulations--Angular-Spectrum-Method

要将它用于您的示例，您需要克隆存储库并将您的双缝图像放在路径./apertures/double_slit.png

然后，从存储库文件夹中运行以下脚本：

from diffractsim import MonochromaticField, mm, nm, cm


F = MonochromaticField(
    wavelength=632.8 * nm, extent_x=5.6 * mm, extent_y=5.6 * mm, Nx=600, Ny=600
)

F.add_aperture_from_image(
   "./apertures/double_slit.jpg", Nx=1500, Ny=1500
)

rgb = F.compute_colors_at(20*cm)
F.plot_colors(rgb, xlim=[-4.5, 4.5], ylim=[-4.5, 4.5])

这是不同波长的近场（20 厘米）脚本的结果：

正如我们在图中看到的，光的波长越高，干涉条纹的长度越宽。

您还可以使用PolychromaticField class 计算具有广谱（例如白光）的衍射图案：

from diffractsim import PolychromaticField, cf, mm, cm


F = PolychromaticField(
    spectrum=2 * cf.illuminant_d65, extent_x=5.6 * mm, extent_y=5.6 * mm, Nx=400, Ny=500
)

F.add_aperture_from_image(
    "./apertures/double_slit.png", Nx=1400, Ny=1400
)

rgb = F.compute_colors_at(20*cm, spectrum_divisions=40)
F.plot_colors(rgb, xlim=[-4.5, 4.5], ylim=[-4.5, 4.5])

导致：

您的代码的问题是您没有正确使用快速傅里叶变换。

我实现的核心（角谱的传播）在MonochromaticFieldpropagate方法中，在monochromatic_simulator.py中：

def propagate(self, z):
    self.z += z

    # compute angular spectrum
    fft_c = fft2(self.E)
    c = fftshift(fft_c)

    kx = np.linspace(-np.pi * self.Nx // 2 / (self.extent_x / 2), np.pi * self.Nx // 2/ (self.extent_x / 2), self.Nx)
    ky = np.linspace(-np.pi * self.Ny // 2 / (self.extent_y / 2), np.pi * self.Ny // 2 / (self.extent_y / 2), self.Ny)
    kx, ky = np.meshgrid(kx, ky)
    kz = np.sqrt((2 * np.pi / self.λ) ** 2 - kx ** 2 - ky ** 2)

    # propagate the angular spectrum a distance z
    E = ifft2(ifftshift(c * np.exp(1j * kz * z)))

    # compute Field Intensity
    self.I = np.real(E * np.conjugate(E)) 

你可以看到我在执行FFT后使用了fftshift的方法来匹配傅里叶变换的定义。

此外，如果您想按照自己的意愿丢弃渐逝场，使用 numpy.where 而不是循环每个像素是一种更好的方法，因为 Python 循环要慢得多：

    # propagate the angular spectrum a distance z
    mask = (2*np.pi/self.λ)**2 - kx**2 - ky**2 > 0
    A = np.where(mask, c*np.exp(1j*kz * z),  0)
    E = ifft2(ifftshift(A))

此代码应替换 compute_colors_at 方法中的最后几行。 希望这会有所帮助！

Answer 2

要做到这一点可能很棘手。这里是：

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
wavenum_sq = wavenum * wavenum

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，然后对结果进行切片。