从惠更斯原理到傅里叶变换用Python模拟光的衍射全过程附代码当一束光穿过狭缝时为什么会在屏幕上形成明暗相间的条纹这个看似简单的现象背后隐藏着从惠更斯原理到傅里叶变换的深刻物理原理。本文将带你用Python代码一步步构建光的衍射模型从基础理论到完整实现让你不仅能理解原理还能亲手模拟出各种衍射图样。1. 光学衍射的物理基础衍射是光波遇到障碍物时偏离直线传播的现象。要理解它我们需要从几个关键概念开始惠更斯原理每个波阵面上的点都可以看作新的球面波源惠更斯-菲涅尔原理这些次级波源的相干叠加决定了光场分布角谱理论将光场分解为不同方向传播的平面波用数学表达惠更斯-菲涅尔原理U(P) (1/iλ) ∬ U(Q) e^(ikr)/r K(θ) dS其中K(θ)是倾斜因子λ是波长r是传播距离。2. 从物理到计算衍射的数值实现2.1 离散化光场在计算机中我们需要将连续的光场离散化。假设孔径平面和观察平面都采样为N×N的网格import numpy as np N 1024 # 采样点数 L 1e-3 # 模拟区域大小(m) dx L/N # 采样间隔 x np.linspace(-L/2, L/2, N) # 坐标轴 X, Y np.meshgrid(x, x) # 二维网格2.2 构建孔径函数定义一个圆形孔径作为示例def circular_aperture(radius): return np.where(X**2 Y**2 radius**2, 1, 0) aperture_radius 100e-6 # 100微米 U circular_aperture(aperture_radius)3. 衍射计算的三种方法3.1 直接积分法按照惠更斯-菲涅尔原理直接计算def direct_integration(U, wavelength, z): k 2*np.pi/wavelength result np.zeros_like(U, dtypecomplex) for i in range(N): for j in range(N): r np.sqrt((X-X[i,j])**2 (Y-Y[i,j])**2 z**2) result U[i,j] * np.exp(1j*k*r)/r return result/(1j*wavelength)这种方法计算量巨大实际中很少使用。3.2 角谱传播法利用傅里叶变换实现高效计算def angular_spectrum(U, wavelength, z): k 2*np.pi/wavelength fx np.fft.fftfreq(N, dx) fy np.fft.fftfreq(N, dx) FX, FY np.meshgrid(fx, fy) # 传递函数 H np.exp(1j*k*z*np.sqrt(1 - (wavelength*FX)**2 - (wavelength*FY)**2)) H[np.isnan(H)] 0 # 处理超出传播范围的频率 # 计算角谱 U_fft np.fft.fft2(U) U_prop np.fft.ifft2(U_fft * H) return U_prop3.3 菲涅尔近似当传播距离足够大时可以使用菲涅尔近似def fresnel_diffraction(U, wavelength, z): k 2*np.pi/wavelength fx np.fft.fftfreq(N, dx) fy np.fft.fftfreq(N, dx) FX, FY np.meshgrid(fx, fy) # 菲涅尔传递函数 H np.exp(1j*k*z) * np.exp(-1j*np.pi*wavelength*z*(FX**2 FY**2)) # 计算衍射场 U_fft np.fft.fft2(U * np.exp(1j*k/(2*z)*(X**2 Y**2))) U_prop np.fft.ifft2(U_fft * H) return U_prop * np.exp(1j*k/(2*z)*(X**2 Y**2))4. 夫琅禾费衍射的特殊情况当满足z ≫ (a²)/λa为孔径尺寸时可以使用夫琅禾费近似def fraunhofer_diffraction(U, wavelength, z): k 2*np.pi/wavelength fx np.fft.fftfreq(N, dx) fy np.fft.fftfreq(N, dx) FX, FY np.meshgrid(fx, fy) # 直接计算傅里叶变换 U_fft np.fft.fftshift(np.fft.fft2(U)) # 缩放因子 scale np.exp(1j*k*z) * np.exp(1j*k/(2*z)*(X**2 Y**2)) / (1j*wavelength*z) return U_fft * scale5. 完整模拟流程与可视化将上述方法整合成完整模拟流程import matplotlib.pyplot as plt # 参数设置 wavelength 500e-9 # 500nm z 0.1 # 传播距离0.1m # 计算衍射场 U_prop angular_spectrum(U, wavelength, z) # 可视化 plt.figure(figsize(10,5)) plt.subplot(121) plt.imshow(np.abs(U), cmapgray) plt.title(Aperture) plt.subplot(122) plt.imshow(np.abs(U_prop)**2, cmapgray) # 光强分布 plt.title(Diffraction Pattern) plt.show()典型衍射图样对比孔径类型菲涅尔衍射夫琅禾费衍射单缝条纹逐渐展宽清晰正弦条纹圆孔同心圆环艾里斑图案方孔复杂过渡态二维sinc函数6. 高级应用与优化技巧6.1 计算效率优化使用FFT加速计算时要注意采样定理dx ≤ λz/(2L)零填充防止混叠GPU加速使用CuPy替代NumPy# 使用零填充的FFT def padded_fft(U, pad_factor2): old_shape U.shape new_shape (old_shape[0]*pad_factor, old_shape[1]*pad_factor) padded np.zeros(new_shape, dtypecomplex) padded[:old_shape[0], :old_shape[1]] U return np.fft.fft2(padded)6.2 相位恢复技术从强度图像重建相位def gerchberg_saxton(I1, I2, z, wavelength, iterations100): # I1: 输入面强度 # I2: 输出面强度 U1 np.sqrt(I1) # 初始猜测 for _ in range(iterations): U2 angular_spectrum(U1, wavelength, z) U2 np.sqrt(I2) * np.exp(1j*np.angle(U2)) # 保持相位 U1_new angular_spectrum(U2, wavelength, -z) U1 np.sqrt(I1) * np.exp(1j*np.angle(U1_new)) return U17. 实际案例CD/DVD衍射分析以CD刻线为例模拟其衍射图样# CD刻线参数 track_pitch 1.6e-6 # 道间距 duty_cycle 0.5 # 占空比 # 创建光栅 grating np.zeros((N,N)) for i in range(N): if i % int(track_pitch/dx) int(duty_cycle*track_pitch/dx): grating[i,:] 1 # 计算衍射 U_CD angular_spectrum(grating, wavelength, 0.3) # 可视化 plt.imshow(np.log(1np.abs(U_CD)**2), cmapviridis) plt.colorbar() plt.title(CD Diffraction Pattern)8. 常见问题与调试技巧出现高频振荡检查采样是否满足Nyquist条件增加零填充结果不对称确保使用fftshift正确显示频谱检查孔径是否居中能量不守恒验证传递函数的单位arity检查边界处理提示调试时先用简单孔径如单缝验证再尝试复杂结构9. 扩展应用方向计算全息设计特定衍射图样的相位板光学测量通过衍射分析表面形貌成像系统模拟透镜的PSF和MTF新型材料超表面衍射特性研究# 计算透镜的PSF def lens_psf(focal_length, wavelength, diameter): # 创建圆形孔径 aperture circular_aperture(diameter/2) # 添加二次相位透镜相位 phase np.exp(-1j*np.pi/(wavelength*focal_length)*(X**2 Y**2)) U_lens aperture * phase # 计算远场衍射 return fraunhofer_diffraction(U_lens, wavelength, focal_length)10. 性能对比与选择指南不同方法的适用场景方法适用距离计算复杂度精度直接积分任意距离O(N⁴)最高角谱传播任意距离O(N²logN)高菲涅尔近似z ≫ a²/λO(N²logN)中等夫琅禾费近似z ≫ a²/λO(N²logN)近场较低选择建议近场10倍瑞利距离角谱法中等距离菲涅尔近似远场夫琅禾费近似11. 现代光学仿真工具链完整的光学仿真通常需要几何光学Zemax, Code V物理光学本文方法或专用软件VirtualLab渲染引擎POV-Ray, Blender Cycles数据分析Python科学计算栈# 与Zemax数据交互示例 def read_zemax_wavefront(filename): # 解析Zemax波前数据 pass def compare_with_zemax(): zemax_data read_zemax_wavefront(lens.zbf) our_result lens_psf(focal_length100e-3, wavelength550e-9, diameter10e-3) plt.figure() plt.subplot(121) plt.imshow(np.abs(zemax_data)**2) plt.subplot(122) plt.imshow(np.abs(our_result)**2) plt.show()12. 从模拟到实验的注意事项实验室观测衍射图样时激光安全始终使用适当功率的激光器对准精度使用精密光学调整架环境振动在光学平台上进行实验探测器校准考虑CCD的响应非线性注意模拟结果与实验的差异可能来自光源的非单色性光学元件像差环境杂散光13. 教育资源与进阶学习推荐学习路径基础教材《光学原理》Born Wolf《傅里叶光学导论》Goodman在线课程MIT OpenCourseWare 光学工程Coursera计算摄影学开源项目PyOptica光学仿真库LightPipes光束传播工具箱学术前沿超表面衍射光学计算成像技术# 使用LightPipes的示例 import LightPipes as lp def lightpipes_example(): size 10*mm wavelength 633*nm N 500 # 创建光束 F lp.Begin(size, wavelength, N) # 添加孔径 F lp.CircAperture(1*mm, 0, 0, F) # 传播 F lp.Fresnel(10*cm, F) # 可视化 lp.Intensity(0, F)