基于python的角谱方法
我正在尝试用角谱方法数值传播一个给定的(电场)场。为此，我阅读了“傅立叶光学的原理和应用”(罗伯特·K·泰森)，第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：
回答 2
Stack Overflow用户
发布于 2021-01-07 03:51:26
我使用Python完成了角谱方法的完整实现：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", pad=(3* mm, 3* mm), Nx=1500, Ny=1500
)
rgb = F.compute_colors_at(20*cm)
F.plot(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", pad=(3* mm, 3* mm), Nx=1400, Ny=1400
)
rgb = F.compute_colors_at(20*cm, spectrum_divisions=40)
F.plot(rgb, xlim=[-4.5, 4.5], ylim=[-4.5, 4.5])这会导致：
代码的问题在于没有正确使用快速傅立叶变换。
我的实现的核心(角谱的传播)在monochromatic_simulator.py中的MonochromaticField propagate方法中
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)) 您可以看到，在执行快速傅立叶变换以匹配傅立叶变换的定义之后，我使用了方法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方法中的最后几行。希望这能对你有所帮助！
Stack Overflow用户
发布于 2020-01-12 03:49:34
要做到这一点可能会很棘手。这就是它：
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，将形状加倍，然后对结果进行切片。
https://stackoverflow.com/questions/58528110
复制相似问题
腾讯云开发者
