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

一束光线穿过一个孔径为 S′S'S′ 的平面，在距离平面为 LLL 的时候，其波函数可以由菲涅尔积分定义：

Ψ(r,t)=C∫S′eik∣r−r′∣∣r−r′∣cos⁡(θ)d2r′,withC=kΨ0e−itw2πi\Psi(\mathbf{r}, t) = C \int_{S'} \frac{e^{ik|r-r'|}}{|r-r'|} \cos(\theta)d^2r', \quad with \quad C = \frac{k \Psi_{0} e^{-itw}}{2 \pi i} Ψ(r,t)=C∫S′∣r−r′∣eik∣r−r′∣cos(θ)d2r′,withC=2πikΨ0e−itw

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

∣r−r′∣≈z+(x−x′)2+(y−y′)22z|r - r'| \approx z + \frac{(x - x')^2 + (y - y')^2}{2z} ∣r−r′∣≈z+2z(x−x′)2+(y−y′)2

并且可以假设：

∣r−r′∣≈L|r - r'| \approx L ∣r−r′∣≈L

菲涅尔积分可以变成：

Ψ(r,t)=R∫−∞∞∫−∞∞f(x′,y′)eik2z(x′2+y′2)e−ikxzx′−ikyzy′dx′dy′\Psi(\mathbf{r}, t) = R \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x', y') e^{\frac{ik}{2z}(x'^2 + y'^2)} e^{-\frac{ikx}{z}x'-\frac{iky}{z}y'} dx'dy'Ψ(r,t)=R∫−∞∞∫−∞∞f(x′,y′)e2zik(x′2+y′2)e−zikxx′−zikyy′dx′dy′

withR=kΨ0ei(kz−tw)2πizeikx2+y22zwith \quad R = \frac{k \Psi_{0} e^{i(kz - tw)}}{2 \pi i z} e^{ik \frac{x^2 + y^2}{2z}} withR=2πizkΨ0ei(kz−tw)eik2zx2+y2

f(x′,y′)={1if (x′,y′)∈S′0if (x′,y′)∉S′f(x', y') = \begin{cases} 1 & \text{ if } (x', y')\in S' \\ 0 & \text{ if } (x', y')\notin S' \end{cases} f(x′,y′)={10 if (x′,y′)∈S′ if (x′,y′)∈/S′

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

Ψ(r,t)=R⋅F[f(x′,y′)eik2z(x′2+y′2)]\Psi(\mathbf{r}, t) = R \cdot \mathcal{F}[f(x', y') e^{\frac{ik}{2z}(x'^2 + y'^2)}] Ψ(r,t)=R⋅F[f(x′,y′)e2zik(x′2+y′2)]

F[f(x′,y′)eik2z(x′2+y′2)]=∫−∞∞∫−∞∞f(x′,y′)eik2z(x′2+y′2)e−ikxzx′−ikyzy′dx′dy′\mathcal{F}[f(x', y') e^{\frac{ik}{2z}(x'^2 + y'^2)}] = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x', y') e^{\frac{ik}{2z}(x'^2 + y'^2)} e^{-\frac{ikx}{z}x'-\frac{iky}{z}y'} dx'dy' F[f(x′,y′)e2zik(x′2+y′2)]=∫−∞∞∫−∞∞f(x′,y′)e2zik(x′2+y′2)e−zikxx′−zikyy′dx′dy′

为了便于计算，我们可以将连续域内的 (x′,y′)(x', y')(x′,y′) 进行离散化，比如，可以将 (x′,y′)∈[−Lx,Lx]×[−Ly,Ly](x', y') \in [-L_x, L_x] \times [-L_y, L_y](x′,y′)∈[−Lx,Lx]×[−Ly,Ly] 划分成等间距的 Nx,NyN_x, N_yNx,Ny 份：

xnx′={−Lx+nx2LxNx,0≤nx≤Nx−1}x'_{n_x} = \left \{ -L_x + n_x \frac{2L_x}{N_x}, \quad 0 \leq n_x \leq N_x -1 \right \} xnx′={−Lx+nxNx2Lx,0≤nx≤Nx−1}

yny′={−Ly+ny2LyNy,0≤ny≤Ny−1}y'_{n_y} = \left \{ -L_y + n_y \frac{2L_y}{N_y}, \quad 0 \leq n_y \leq N_y -1 \right \} yny′={−Ly+nyNy2Ly,0≤ny≤Ny−1}

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

G(ux,uy)=∑nx=0Nx−1∑ny=0Ny−1g(xnx′,yny′)e−i2πuxnxNx−i2πuynyNyG(u_x, u_y) = \sum_{n_x = 0}^{N_x -1} \sum_{n_y =0}^{N_y -1} g(x'_{n_x}, y'_{n_y})e^{-i2\pi u_x \frac{n_x}{N_x} - i 2 \pi u_y \frac{n_y}{N_y} } G(ux,uy)=nx=0∑Nx−1ny=0∑Ny−1g(xnx′,yny′)e−i2πuxNxnx−i2πuyNyny

在距离为 L 的平面上，

xux=(ux−Nx/2)Lλ2Lxx_{u_x} = \frac{(u_x - N_x / 2)L \lambda}{2 L_x} xux=2Lx(ux−Nx/2)Lλ

yux=(uy−Ny/2)Lλ2Lyy_{u_x} = \frac{ (u_y - N_y / 2)L \lambda}{ 2 L_y} yux=2Ly(uy−Ny/2)Lλ


### 代码转载自：
##  https://rafael-fuente.github.io/solving-the-diffraction-integral-with-the-fast-fourier-transform-fft-and-python.htmlimport numpy as np
import matplotlib.pyplot as pltclass 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.4Nx= 2500
Ny= 1500sheet = Sheet(extentX = [-Lx, Lx] , extentY = [-Ly, Ly], Nx= Nx, Ny= Ny)#slit separation
mm = 1e-3
D = 500 * mmsheet.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

利用 FFT 模拟菲涅尔衍射积分相关推荐

基于Matlab——夫琅禾夫衍射以及菲涅尔衍射
我以往在学习F分析的时候,编写了一个作业代码,在此附上供大家学习交流. 引言: 在傅里叶光学信息基础中,主要研究的是光在传播过程携带的信息如何去检测得到.如果光在自由空间(均匀 ...
matlab菲涅尔衍射_有问必答——SYNOPSYS安装体验课堂——可以设计菲涅尔透镜吗?...
问:SYNOPSYS可以设计菲涅尔透镜吗? 答:在USS中有多种菲涅尔面型,用户输入参数即可. 问:SYNOPSYS中具有的输入方式? 答:大家总是有个误区,以为SYNOPSYS需要输入命令运行,其实 ...
基于Matlab模拟菲涅尔公式
✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.
Chango的数学Shader世界（七）水波模拟-透明水面，菲涅尔(Fresnel)效应
目的: 解析,改进,批评一个国外免费透明水面Shader,进一步了解Shader背后的物理原理. 参考: 菲涅尔反射分析: 我将原水面Shader一再简化,从中抽取最主要的部分,忽略细枝末节,并改掉 ...
【scipy】菲涅尔积分和羊角螺线
文章目录菲涅尔积分简介复数域的菲涅尔积分羊角螺线菲涅尔积分简介菲涅尔积分一般被写为 S ( x ) S(x) S(x)和 C ( x ) C(x) C(x),定义为 S ( x ) = ∫ ...
Unity Shader 水多种元素的实现（反射、折射、菲涅尔、深浅、浪花/泡沫、水波、可交互）
综合效果经过各元素叠加和程序的审美调参后的综合效果交互的水波与边缘浪花的合并需要优化一下反射两种方案: cubeMap 以水面对称设一个摄像机 cubeMap 实现:反射探针生成Cube ...
Unity Shader：实现菲涅尔+色散效果的环境映射以及相关原理解析
文章目录 1,色散在光学中的原理 2,反射的数学计算方法以及用它实现环境映射 3,折射的原理以及色散的实现 4,菲涅尔效果 5,拥有菲涅尔与色散效果的环境映射 1,色散在光学中的原理复色光 --现实 ...
Unity Shader：实现菲涅尔+色散效果以及相关原理解析
1,色散在光学中的原理 2,反射的原理以及环境映射的实现 3,折射的原理以及色散的实现 4,菲涅尔效果 5,将菲涅尔与色散效果增加到环境映射中 1,色散在光学中的原理复色光 --现实生活中 ...
matlab波带片程序,Matlab编程快速实现振幅型菲涅尔波带片的设计
维普资讯 http://doc.docsou.com 第 8卷第 1 5期 20 0 8年 8月科学技术与工程 V 1 8 No 5 o. .1 Au .2 0 g 08 17-89 ...
matlab 伽马波带片,菲涅尔波带片
满意答案 mbtf4535 2018.11.28 采纳率:52% 等级:8 已帮助:310人菲涅尔透镜是是带通光学滤镜,被广泛使用在警报器和相机对焦屏上.菲涅尔波带片既有会聚透镜的作用,又有发 ...

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

利用 FFT 模拟菲涅尔衍射积分相关推荐

最新文章

热门文章