利用 FFT 模拟菲涅尔衍射积分

利用 FFT 模拟菲涅尔衍射积分

一束光线穿过一个孔径为 $S^{\prime }$ 的平面，在距离平面为 $L$ 的时候，其波函数可以由菲涅尔积分定义：

$$\Psi (\mathbf{r},t)=C{\int }_{S^{\prime }}\frac{e^{ik|r-r^{\prime }|}}{|r-r^{\prime }|}\cos (\theta )d^2{r}^{\prime },withC=\frac{k{\Psi }_0{e}^{-itw}}{2\pi i}$$

基于菲涅尔近似，在角度接近 0 度的时候，上式可以简化为：

$$|r-r^{\prime }|\approx z+\frac{(x-x^{\prime }{)}^2+(y-y^{\prime }{)}^2}{2z}$$

并且可以假设：

$$|r-r^{\prime }|\approx L$$

菲涅尔积分可以变成：

$$\Psi (\mathbf{r},t)=R{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x^{\prime },y^{\prime })e^{\frac{ik}{2z}(x^{\prime 2}+y^{\prime 2})}{e}^{-\frac{ikx}{z}{x}^{\prime }-\frac{iky}{z}{y}^{\prime }}d{x}^{\prime }d{y}^{\prime }$$

$$withR=\frac{k{\Psi }_0{e}^{i(kz-tw)}}{2\pi iz}{e}^{ik\frac{x^2+y^2}{2z}}$$

$$f(x^{\prime },y^{\prime })=\{\begin{array}{ll}1& \text{ if }(x^{\prime },y^{\prime })\in S^{\prime }\\ 0& \text{ if }(x^{\prime },y^{\prime })\notin S^{\prime }\end{array}$$

上式可以转换成傅里叶变换的形式：

$$\Psi (\mathbf{r},t)=R\cdot \mathcal{F}[f(x^{\prime },y^{\prime })e^{\frac{ik}{2z}(x^{\prime 2}+y^{\prime 2})}]$$

$$\mathcal{F}[f(x^{\prime },y^{\prime })e^{\frac{ik}{2z}(x^{\prime 2}+y^{\prime 2})}]={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x^{\prime },y^{\prime })e^{\frac{ik}{2z}(x^{\prime 2}+y^{\prime 2})}{e}^{-\frac{ikx}{z}{x}^{\prime }-\frac{iky}{z}{y}^{\prime }}d{x}^{\prime }d{y}^{\prime }$$

为了便于计算，我们可以将连续域内的 $(x^{\prime },y^{\prime })$ 进行离散化，比如，可以将 $(x^{\prime },y^{\prime })\in [-L_x,L_x]\times [-L_y,L_y]$ 划分成等间距的 $N_x,N_y$ 份：

$$x_{n_x}^{\prime }=\left\{-L_x+n_x\frac{2L_x}{N_x},0\le n_x\le N_x-1\right\}$$

$$y_{n_y}^{\prime }=\left\{-L_y+n_y\frac{2L_y}{N_y},0\le n_y\le N_y-1\right\}$$

接下来，就可以定义离散傅里叶变换：

$$G(u_x,u_y)=\sum _{n_x=0}^{N_x-1}\sum _{n_y=0}^{N_y-1}g(x_{n_x}^{\prime },y_{n_y}^{\prime })e^{-i2\pi u_x\frac{n_x}{N_x}-i2\pi u_y\frac{n_y}{N_y}}$$

在距离为 L 的平面上，

$$x_{u_x}=\frac{(u_x-N_x/2)L\lambda }{2L_x}$$

$$y_{u_x}=\frac{(u_y-N_y/2)L\lambda }{2L_y}$$

### 代码转载自：
##  https://rafael-fuente.github.io/solving-the-diffraction-integral-with-the-fast-fourier-transform-fft-and-python.html

import numpy as np
import matplotlib.pyplot as plt

class Sheet():
    def __init__(self,extentX, extentY, Nx, Ny):
        self.x = np.linspace(extentX[0],extentX[1],Nx)
        self.y = np.linspace(extentY[0],extentY[1],Ny)
        self.xx,self.yy = np.meshgrid(self.x, self.y)

        self.Nx = int(Nx)
        self.Ny = int(Ny)
        self.f = np.zeros((int(self.Ny), int(self.Nx)))

    def rectangular_slit(self,x0, y0, lx, ly):
        """
        Creates a slit centered at the point (x0, y0) with width lx and height ly
        """
        self.f += np.select( [((self.xx > (x0 - lx/2) ) & (self.xx < (x0 + lx/2) )) & ((self.yy > (y0 - ly/2) ) & (self.yy < (y0 + ly/2) )),  True], [1, 0])
        
Lx = 1.4
Ly = 0.4

Nx= 2500
Ny= 1500

sheet = Sheet(extentX = [-Lx, Lx] , extentY = [-Ly, Ly], Nx= Nx, Ny= Ny)

#slit separation 
mm = 1e-3
D = 500 * mm

sheet.rectangular_slit(x0 = -D/2, y0 = 0, lx = 22 * mm , ly = 88 * mm)
sheet.rectangular_slit(x0 = +D/2, y0 = 0, lx = 22 * mm , ly = 88 * mm)


# distance from slit to the screen (mm)
z = 5000

# wavelength (mm)
λ = 18.5*1e-7
k = 2*np.pi/λ

fft_c = np.fft.fft2(sheet.f * np.exp(1j * k/(2*z) *(sheet.xx**2 + sheet.yy**2)))
c = np.fft.fftshift(fft_c)

fig = plt.figure(figsize=(6, 6))
ax1 = fig.add_subplot(2,1,1)  
ax2 = fig.add_subplot(2,1,2,sharex=ax1, yticklabels=[])

abs_c = np.absolute(c)

#screen size mm
Lx_screen = Nx*z*λ/(4*Lx)
Ly_screen = Ny*z*λ/(4*Ly)

x_max = (np.pi/Lx * (Nx//2 - 1))*z*λ/(2*np.pi)
y_max = (np.pi/Ly * (Ny//2 - 1))*z*λ/(2*np.pi)

ax1.imshow(abs_c, extent = [-Lx_screen, Lx_screen, -Ly_screen,Ly_screen], cmap ='gray', interpolation = "bilinear", aspect = 'auto')

ax2.plot(np.linspace(-Lx_screen,Lx_screen, len(abs_c[0])), abs_c[len(abs_c)//2]**2)

ax1.set_ylabel("y (mm)")
ax2.set_xlabel("x (mm)")
ax2.set_ylabel("Probability Density $|\psi|^{2}$")
ax1.set_xlim([-2,2])
ax2.set_xlim([-2,2])
ax1.set_ylim([-1,1])
plt.setp(ax1.get_xticklabels(), visible=False)
plt.show()

参考文献：
https://rafael-fuente.github.io/solving-the-diffraction-integral-with-the-fast-fourier-transform-fft-and-python.html