对衍射最经典的解释是Huygens-Fresnel原理，Huygens认为波阵面上每一点都会成为新的波源，这些子波源的相互干涉就形成了衍射。这显然是一种离散的观点，仿佛是专门为程序员准备的一样。
假设一束光打在一个方形孔上，这个方形孔被细分成
n
×
n
n\times n
n×n个网格，那么每个网格都相当于是一个小孔，而这些小孔的互相干涉，即为衍射。随着网格不断被细分，最终可以逼近真实的衍射情形。那么，假设矩孔处为等相位面，其网格坐标为
(
i
,
j
)
(i,j)
(i,j)，到衍射屏距离为
d
d
d,那么对于衍射屏上任意一点
P
(
x
,
y
)
P(x,y)
P(x,y)，其光强为
I
(
x
,
y
)
=
I
0
n
2
∑
i
=
0
,
j
=
0
n
cos
(
2
π
(
x
−
i
)
2
+
(
y
−
j
)
2
+
d
2
λ
)
I(x,y)=\frac{I_0}{n^2}\displaystyle\sum_{i=0,j=0}^{n}\cos(\frac{2\pi\sqrt{(x-i)^2+(y-j)^2+d^2}}{\lambda})
I(x,y)=n2I0i=0,j=0∑ncos(λ2π(x−i)2+(y−j)2+d2)
如果这束光来自点光源，则球面波打在方孔之后，考虑到不同格点到光源的距离不同，则其相位亦有不同。假设光源坐标
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0)到矩孔平面的距离为
L
L
L，则
(
i
,
j
)
(i,j)
(i,j)点的相位为
φ
(
i
,
j
)
=
2
π
(
x
0
−
i
)
2
+
(
y
0
−
j
)
2
+
L
2
λ
\varphi(i,j)=\frac{2\pi\sqrt{(x_0-i)^2+(y_0-j)^2+L^2}}{\lambda}
φ(i,j)=λ2π(x0−i)2+(y0−j)2+L2
另外，光强在矩孔上也不再是均匀分布，而是变为
I
(
i
,
j
)
=
I
0
(
i
−
x
0
)
2
+
(
j
−
y
0
)
2
+
L
2
I(i,j)=\frac{I_0}{\sqrt{(i-x_0)^2+(j-y_0)^2+L^2}}
I(i,j)=(i−x0)2+(j−y0)2+L2I0
同理，我们刚刚写下的平面矩孔光强叠加也出现了问题，如果我们默认矩孔上每个格点都是一个点光源，那么打在衍射屏上之后，应该遵从球面波的衰减原则。又因为这种假设其实忽略了从光源射到矩孔过程中的光线的传播方向，所以应该有一个倾斜因子，即
(
i
,
j
)
(i,j)
(i,j)点打在
(
x
,
y
)
(x,y)
(x,y)点的光强为
I
i
j
(
x
,
y
)
=
K
I
(
i
,
j
)
(
i
−
x
)
2
+
(
j
−
y
)
2
+
L
2
I_{ij}(x,y)=K\frac{I(i,j)}{\sqrt{(i-x)^2+(j-y)^2+L^2}}
Iij(x,y)=K(i−x)2+(j−y)2+L2I(i,j)
其中，
K
K
K为倾斜因子。假设一束光在前进的过程中经过一个小孔，然后分裂成了许多束光，那么这些光应该在传播方向上最强，越远离传播方向越弱，但必须为正值，并且在整个空间中的积分为1。也就是说这个
K
K
K需要满足
{
K
m
a
x
=
k
(
0
)
∫
−
π
2
π
2
K
(
θ
)
d
θ
=
1
。
\left\{Kmax=k(0)∫π2−π2K(θ)dθ=1。
\right.
⎩⎪⎨⎪⎧Kmax=∫−2π2πK(θ)dθ=k(0)1。
则
1
2
cos
θ
\frac{1}{2}\cos\theta
21cosθ满足要求，这就是当年Fresnel所猜测的倾斜因子的值。当矩孔的尺度远小于光源、矩孔、衍射屏之间的距离时，这个值几乎等于1。
记
L
i
j
=
(
i
−
x
0
)
2
+
(
j
−
y
0
)
2
+
L
2
d
i
j
(
x
,
y
)
=
(
x
−
i
)
2
+
(
y
−
j
)
2
+
d
2
Lij=√(i−x0)2+(j−y0)2+L2dij(x,y)=√(x−i)2+(y−j)2+d2
Lij=dij(x,y)=(i−x0)2+(j−y0)2+L2(x−i)2+(y−j)2+d2
则可将点光源在衍射屏
(
x
,
y
)
(x,y)
(x,y)处的光强写为
I
(
x
,
y
)
=
∑
i
=
0
,
j
=
0
n
I
0
L
i
j
⋅
d
i
j
(
x
,
y
)
cos
(
2
π
(
d
i
j
(
x
,
y
)
+
L
i
j
)
λ
)
I(x,y)=n∑i=0,j=0I0Lij⋅dij(x,y)cos(2π(dij(x,y)+Lij)λ)
I(x,y)=i=0,j=0∑nLij⋅dij(x,y)I0cos(λ2π(dij(x,y)+Lij))
#基尔霍夫衍射,衍射屏坐标范围-dGrid:dGrid,光源坐标(0,0)
#简单的矩孔衍射,dSource为光源到小孔的距离；dScreen为衍射屏到小孔距离
#dHole为矩孔网格尺寸；dGrid为衍射屏网格尺寸；nGrid为网格数目
def squareDiff(dSource=1,dScreen=1.4,dWave=1e-6,
dHole=3e-5,nGrid=100,dGrid=1e-5):
nX,nY = nGrid*np.array([1,1])
axisX = np.arange(-nX,nX+1)*dGrid
axisY = np.arange(-nY,nY+1)*dGrid
xAxis,yAxis = np.meshgrid(axisX,axisY) #此为衍射屏的x坐标
axisX = np.arange(-nX,nX+1)*dHole
axisY = np.arange(-nY,nY+1)*dHole
xHole,yHole = np.meshgrid(axisX,axisY) #此为矩孔的x坐标
dArrS = np.sqrt(xHole**2+yHole**2+dSource**2) #孔平面到光源的距离
nSide = int(nGrid*2+1) #格点个数
pane = np.zeros([nSide,nSide]) #衍射屏强度
for m in range(nSide):
for n in range(nSide):
dArr = np.sqrt((xHole-xAxis[m,n])**2+(yHole-yAxis[m,n])**2+dScreen**2)
pane[m,n] = np.sum(
np.cos(np.pi*2*(dArr+dArrS)/dWave)/dArr/dArrS)
pane = np.abs(pane)
pane = pane/np.max(np.array(pane))
fig = plt.figure()
ax = axd(fig)
ax.plot_surface(xAxis,yAxis,pane)
plt.show()
return pane, xAxis, yAxis
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
其结果为
分析上面的代码可知，当小孔、衍射屏均为
n
×
n
n\times n
n×n的网格时，需要运行
n
2
n^2
n2次循环，每次循环需要运行
n
2
n^2
n2次三角函数，
3
n
2
3n^2
3n2次除法，
3
n
2
3n^2
3n2次乘法，
5
n
2
5n^2
5n2次加减法，
n
2
n^2
n2次开方。其时间复杂度为
O
(
n
4
)
O(n^4)
O(n4)，效率相当低下。
我们发现当矩孔和衍射屏的网格划分方式相同时，
d
i
j
(
x
,
y
)
d_{ij}(x,y)
dij(x,y)中包含大量的重复数据，如果能够运算一次，反复调用，应该能够为衍射计算提供良好的优化。例如，假设衍射屏和矩孔都是
n
×
n
n\times n
n×n网格，那么衍射屏上第一个数据为
d
i
j
(
1
,
1
)
d_{ij}(1,1)
dij(1,1)矩阵，这个矩阵中包含此后衍射屏上所有的数据。即
d
i
j
(
1
,
1
)
=
[
d
11
,
11
d
12
,
11
…
d
1
n
,
11
d
21
,
11
d
22
,
11
…
d
2
n
,
11
⋮
d
n
1
,
11
d
n
2
,
11
…
d
n
n
,
11
]
d_{ij}(1,1)=[d11,11d12,11…d1n,11d21,11d22,11…d2n,11⋮dn1,11dn2,11…dnn,11]
dij(1,1)=⎣⎢⎢⎢⎡d11,11d21,11dn1,11d12,11d22,11dn2,11……⋮…d1n,11d2n,11dnn,11⎦⎥⎥⎥⎤
则
d
i
j
(
x
,
y
)
=
d
(
∣
x
−
i
+
1
∣
,
∣
y
−
j
+
1
∣
)
,
11
d_{ij}(x,y)=d_{(|x-i+1|,|y-j+1|),11}
dij(x,y)=d(∣x−i+1∣,∣y−j+1∣),11
故可定义矩阵索引
#输入对于M*M矩阵的第一个值到N*N矩阵的距离，返回M(m,n)的距离矩阵
def getDisMat(dMat,N,m,n):
dMat = np.mat(dMat)
A = dMat[1:m,1:n]
B = dMat[1:m,0:N-n+1]
C = dMat[0:N-m+1,1:n]
D = dMat[0:N-m+1,0:N-n+1]
return np.vstack((np.hstack((np.flip(A),np.flipud(B))),np.hstack((np.fliplr(C),D)))) #stack矩阵拼接，flip翻转矩阵
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
相应地算法改为(其他位置不变)
dArrS = np.sqrt(xHole**2+yHole**2+dSource**2) #孔平面到光源的距离
dScreen = np.sqrt( #衍射平上第(0,0)个点的距离矩阵
(xHole-xAxis[0,0])**2+(yHole-yAxis[0,0])**2+dScreen**2)
nSide = int(nGrid*2+1) #格点个数
pane = np.zeros([nSide,nSide]) #衍射屏强度
for m in range(nSide):
for n in range(nSide):
dArr = getDisMat(dScreen,nSide,m,n)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
但这种优化是极其有限的，这是直观无脑的思维方式所带来的麻烦。在接下来的傅里叶光学中，我们将继续处理衍射的计算问题。