1. 从一次仿真翻车说起为什么要啃圆孔菲涅尔衍射几年前帮一个做光学检测的朋友复现一个孔径衍射实验他信誓旦旦地说“近场衍射嘛套个夫琅禾费公式就完事了”结果仿真出来的光斑和实验拍到的完全对不上——实验图中心亮斑周围有一圈圈明暗相间的同心环而他的仿真结果只有一团模糊的亮斑。问题出在哪他把菲涅尔衍射当成了夫琅禾费衍射来处理。这两者的区别恰恰是圆孔衍射仿真里最容易踩的坑。圆孔菲涅尔衍射这件事说白了就是一束光打在一个不透明屏上的圆形小孔上在孔后有限距离处观察到的光强分布。它和远场的夫琅禾费衍射最大的不同在于观察屏离孔的距离不能近似为无穷远波前曲率不能忽略必须用菲涅尔-基尔霍夫衍射积分或者菲涅尔半波带法来处理。这个内容适合谁看光学、物理、光电信息专业的本科生做课程设计研究生做仿真验证以及做光学检测、激光光束分析的工程师做快速原型验证。MATLAB在这里的角色就是一个“数值实验台”你不需要真的去搭光路就能把不同孔径、不同波长、不同传播距离下的衍射图样算出来。我自己的体会是圆孔菲涅尔衍射这个题目看起来经典但真正动手写代码的时候坑比想象的多积分区域怎么离散、采样间隔怎么选、FFT用不用、坐标网格怎么建、强度归一化怎么做每一步都有讲究。下面我把整个项目的设计思路、核心原理、代码实现和踩坑经验完整地拆一遍。2. 整体设计思路为什么选MATLAB而不是其他工具2.1 菲涅尔衍射的数学本质与计算策略选型圆孔菲涅尔衍射的数学表达核心是菲涅尔-基尔霍夫衍射积分U(P) (1/(jλ)) ∬ U₀(Q) · exp(jkr) / r · K(θ) dS其中U₀(Q)是孔径平面上的光场分布r是孔径上某点到观察点的距离K(θ)是倾斜因子λ是波长。对于圆孔这种具有轴对称性的孔径理论上可以用贝塞尔函数解析求解但一旦涉及到偏心观察、多孔叠加、或者孔径内有相位调制解析法就力不从心了。所以数值计算是更通用的路线。数值计算有两条路可走直接积分法和FFT传播法。直接积分法就是老老实实把孔径平面划分成网格对每个观察点做二重数值积分。优点是物理意义清晰观察平面的坐标可以任意设定不受FFT网格限制缺点是计算量大如果孔径网格取500×500观察点取200×200那就是200×200×500×500次运算MATLAB里跑起来很慢。FFT传播法的思路是利用菲涅尔衍射的卷积形式把衍射积分转化成傅里叶变换用FFT加速。优点是快缺点是观察平面的采样间隔和范围由FFT的网格决定不够灵活而且容易产生混叠。我最终选的是直接积分法为主、FFT法做交叉验证的策略。原因很简单圆孔菲涅尔衍射的观察区域通常不会太大直接积分法在MATLAB里用矩阵化运算优化之后速度完全可以接受而且直接积分法的物理直观性强适合教学和调试。FFT法用来做快速验证两者结果对上了心里就有底了。2.2 为什么不用现成的光学仿真软件有人可能会问Zemax、VirtualLab这些专业光学软件不香吗香但有两个问题。第一这些软件对菲涅尔衍射的处理往往封装得很深你很难看到中间过程对于理解物理本质帮助有限。第二做课程设计或者快速验证的时候装一个几GB的软件、学一套复杂的操作界面时间成本太高。MATLAB的优势在于代码透明、参数可调、结果可视化灵活而且大部分理工科学生已经有MATLAB基础上手门槛低。你可以在几十行代码里把整个物理过程完整地表达出来这是专业软件做不到的。另外MATLAB的矩阵运算天然适合做网格化数值计算。把孔径平面和观察平面都离散成二维矩阵衍射积分就变成了矩阵乘法或者逐元素运算代码写起来非常简洁。配合imagesc、surf、plot这些可视化函数光强分布、相位分布、截面曲线都能快速画出来。2.3 参数体系的搭建逻辑做仿真之前必须先确定一套合理的参数体系。圆孔菲涅尔衍射涉及的核心参数包括参数名符号典型取值物理意义波长λ632.8 nm决定衍射尺度圆孔半径a0.5~2 mm孔径尺寸传播距离z10~200 cm观察屏到孔的距离孔径平面采样数N₁501积分精度观察平面采样数N₂301输出分辨率孔径平面尺寸L₁4a~6a覆盖孔径及边缘观察平面尺寸L₂根据z调整覆盖衍射图样这里最关键的是菲涅尔数N_F a²/(λz)。这个无量纲数决定了衍射所处的区域N_F 1时接近几何光学N_F ~ 1时是菲涅尔衍射的典型区域N_F 1时过渡到夫琅禾费衍射。我在设计参数时会先算一下菲涅尔数确保它落在感兴趣的范围内。比如a1mmλ632.8nmz50cm算出来N_F ≈ 3.16属于典型的菲涅尔衍射区域中心会出现明显的泊松亮斑和同心环。3. 核心原理拆解菲涅尔衍射到底在算什么3.1 从惠更斯-菲涅尔原理到菲涅尔-基尔霍夫积分惠更斯-菲涅尔原理的表述很直观波前上每一点都可以看作一个新的子波源这些子波源发出的球面波在空间某点相干叠加就形成了该点的光场。菲涅尔在这个基础上引入了干涉的概念但缺乏严格的数学基础。基尔霍夫后来用格林定理严格推导了衍射积分公式给出了倾斜因子K(θ) (1cosθ)/2。对于圆孔菲涅尔衍射孔径函数可以写成P(x₀, y₀) 1, 当 x₀² y₀² ≤ a² P(x₀, y₀) 0, 其他观察平面上的光场就是孔径函数与菲涅尔传播核的卷积。在傍轴近似下r可以展开为r ≈ z [(x-x₀)² (y-y₀)²] / (2z)这个近似成立的条件是z³ [(x-x₀)² (y-y₀)²]² / (8λ)对于毫米级孔径和厘米级传播距离这个条件通常满足。3.2 菲涅尔半波带法的直观解释菲涅尔半波带法是一个很好的辅助理解工具。把圆孔分成若干个同心环带每个环带到观察点的光程差为半个波长。相邻半波带在观察点产生的光场相位相反所以它们相互抵消。如果圆孔恰好包含奇数个半波带中心点就是亮点如果包含偶数个半波带中心点就是暗点。这就是为什么改变传播距离z中心点的光强会周期性地明暗变化。半波带的半径公式是ρ_k sqrt(kλz (kλ/2)²) ≈ sqrt(kλz)第k个半波带的面积近似相等所以每个半波带对中心点的贡献幅度基本相同。这个结论在写代码验证的时候非常有用你可以通过改变z观察中心点光强随菲涅尔数的振荡来验证你的仿真是否正确。3.3 数值离散化中的关键细节直接积分法最核心的问题是如何离散化。孔径平面的网格间距Δ₁和观察平面的网格间距Δ₂需要满足采样定理。对于菲涅尔衍射相位因子exp(jk[(x-x₀)²(y-y₀)²]/(2z))在孔径边缘变化最快所以网格间距需要足够小使得相位变化不超过π。粗略的估计是Δ₁ ≤ λz / (2a)以λ632.8nmz50cma1mm为例Δ₁ ≤ 0.158mm。如果孔径平面尺寸取6mm那么N₁至少需要6/0.158 ≈ 38个点。实际中为了精度我会取N₁501Δ₁0.012mm远小于临界值确保相位采样充分。观察平面的网格间距Δ₂也有类似的要求但通常观察平面的尺寸比孔径大所以Δ₂可以适当放宽。不过如果观察平面太大边缘区域的相位变化也会加快需要相应增加N₂。注意很多人写代码时只关注孔径平面的采样忽略了观察平面的采样结果边缘区域出现明显的数值误差。我的经验是观察平面的采样数不要低于201如果观察范围超过孔径的10倍建议加到501。4. 实操过程从零搭建圆孔菲涅尔衍射仿真4.1 环境准备与基础参数设定MATLAB版本建议R2018b以上主要用到的是基础矩阵运算和绘图函数不需要额外的工具箱。如果你想用FFT法做验证需要Signal Processing Toolbox但直接积分法不需要任何工具箱。先定义基础参数% 基础物理参数 lambda 632.8e-9; % 波长单位米 k 2 * pi / lambda; % 波数 a 1e-3; % 圆孔半径单位米 z 0.5; % 传播距离单位米 % 孔径平面网格 N1 501; % 采样点数 L1 6 * a; % 孔径平面尺寸 x1 linspace(-L1/2, L1/2, N1); y1 x1; [X1, Y1] meshgrid(x1, y1); dx1 x1(2) - x1(1); % 观察平面网格 N2 301; L2 0.02; % 观察平面尺寸根据实际情况调整 x2 linspace(-L2/2, L2/2, N2); y2 x2; [X2, Y2] meshgrid(x2, y2); % 圆孔孔径函数 aperture double(X1.^2 Y1.^2 a^2);这里有几个细节值得说。第一孔径平面尺寸取6a而不是2a是为了让孔径边缘的场分布有足够的过渡区域避免周期性边界效应。第二观察平面尺寸L2需要根据传播距离z和孔径a来估算一个经验公式是L2 ≈ 2λz/a 2a。以当前参数算2×632.8e-9×0.5/1e-3 2e-3 ≈ 2.63mm所以L2取20mm已经足够覆盖主瓣和几个旁瓣了。4.2 直接积分法的矩阵化实现直接积分法的朴素实现是双重循环但在MATLAB里这样写会非常慢。正确的做法是利用矩阵化运算把二重积分转化成矩阵乘法。% 预分配观察平面光场 U2 zeros(N2, N2); % 矩阵化计算 % 将孔径平面坐标拉成列向量 x1_vec X1(:); y1_vec Y1(:); ap_vec aperture(:); % 对每个观察点计算 for m 1:N2 for n 1:N2 x_obs X2(m, n); y_obs Y2(m, n); % 计算距离r r sqrt((x_obs - x1_vec).^2 (y_obs - y1_vec).^2 z^2); % 倾斜因子 cos_theta z ./ r; K (1 cos_theta) / 2; % 积分核 integrand ap_vec .* exp(1j * k * r) ./ r .* K; % 数值积分矩形法 U2(m, n) sum(integrand) * dx1^2 / (1j * lambda); end end % 光强 I2 abs(U2).^2; I2_norm I2 / max(I2(:));这段代码虽然有两层循环但内层是向量化运算实际跑起来在普通笔记本上大约需要十几秒到几十秒取决于N2的大小。如果嫌慢可以把观察平面的循环也向量化但内存消耗会急剧增加。我的建议是先用小尺寸N2101调试确认结果正确后再加大N2。实操心得r的计算中z²项不能省略。很多人为了省事直接用z代替r这在傍轴近似下勉强可以但倾斜因子K的计算会出错导致边缘区域的强度偏差。我实测过忽略z²项在观察平面边缘会造成约5%的强度误差。4.3 光强分布的可视化与截面分析算完光强之后可视化是关键。我通常画三张图二维伪彩色图、三维曲面图、中心截面曲线。% 二维伪彩色图 figure; imagesc(x2*1e3, y2*1e3, I2_norm); colormap(hot); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(圆孔菲涅尔衍射光强分布); % 中心截面曲线 figure; plot(x2*1e3, I2_norm(ceil(N2/2), :), b-, LineWidth, 1.5); xlabel(x (mm)); ylabel(归一化光强); title(中心截面光强分布); grid on; % 三维曲面 figure; surf(x2*1e3, y2*1e3, I2_norm); shading interp; colormap(jet); xlabel(x (mm)); ylabel(y (mm)); zlabel(归一化光强); title(三维光强分布);截面曲线是最有用的分析工具。你可以从曲线上读出中心亮斑的宽度、第一暗环的位置、旁瓣的峰值强度。这些数据和理论值对比就能验证仿真的正确性。比如第一暗环的位置理论上对应贝塞尔函数J₁的第一个零点大约在1.22λz/(2a)处。以当前参数算1.22×632.8e-9×0.5/(2×1e-3) ≈ 0.193mm。你可以在截面曲线上找第一个极小值的位置看看是否接近这个值。4.4 FFT传播法的交叉验证为了确认直接积分法的结果可靠我用FFT法做一次交叉验证。FFT法的核心是把菲涅尔衍射积分写成卷积形式然后用FFT加速。% FFT法参数 M 1024; % FFT点数 L 0.02; % 观察平面尺寸 dx L / M; x_fft (-M/2 : M/2-1) * dx; [X_fft, Y_fft] meshgrid(x_fft, x_fft); % 传递函数 H exp(1j * k * z) / (1j * lambda * z) * ... exp(1j * k / (2*z) * (X_fft.^2 Y_fft.^2)); % 孔径频谱 U1_fft fftshift(fft2(ifftshift(aperture))); % 传播 U2_fft ifftshift(ifft2(fftshift(H .* U1_fft))); % 光强 I2_fft abs(U2_fft).^2; I2_fft_norm I2_fft / max(I2_fft(:));FFT法的结果和直接积分法在中心区域应该高度一致边缘区域可能有轻微差异主要来自FFT的周期性边界效应。如果差异很大说明采样不足或者观察平面尺寸设置不合理。5. 常见问题与排查技巧实录5.1 衍射图样不对称或出现条纹这是最常见的现象。原因通常有三个一是孔径平面的网格不对称比如N1取了偶数导致网格中心不在原点二是观察平面的坐标定义有误比如用了linspace(0, L, N)而不是linspace(-L/2, L/2, N)三是FFT法中的fftshift和ifftshift用反了。排查方法很简单先检查孔径函数aperture的对称性用sum(aperture, 1)看每一列的和是否关于中心对称。然后检查观察平面坐标是否关于零对称。最后检查FFT的移位操作是否正确。避坑技巧N1和N2都取奇数这样网格中心恰好落在原点对称性最好。我一开始用偶数调了半天才发现问题出在这里。5.2 中心光强随传播距离的振荡不符合预期菲涅尔衍射的一个经典特征是中心点光强随传播距离z周期性振荡周期对应菲涅尔数的变化。如果你算出来的曲线是单调的或者振荡周期不对大概率是菲涅尔数算错了或者传播距离的范围设置不合理。排查步骤先手算几个z值对应的菲涅尔数N_F a²/(λz)确认它们覆盖了从大到小的范围。然后检查代码中z的单位是否统一都是米。最后检查积分公式中的相位因子是否正确特别是k*r这一项r的计算是否包含了z²。5.3 计算速度太慢的优化方案直接积分法的计算复杂度是O(N2² × N1²)当N2301、N1501时运算量大约是2.3×10¹⁰次MATLAB里跑起来确实慢。优化方案有几个第一利用圆孔的轴对称性把二维积分简化为一维积分。对于圆孔观察平面上的光场只依赖于径向坐标ρ sqrt(x²y²)所以可以用一维的贝塞尔函数积分来代替二维积分。这样计算量降到O(N2 × N1)速度提升几个数量级。第二用GPU加速。如果你的MATLAB支持Parallel Computing Toolbox把gpuArray加到矩阵上速度能提升10倍以上。第三降低观察平面的采样数。对于初步验证N2101就够了确认结果正确后再加大。5.4 常见问题速查表问题现象可能原因排查方法解决方案图样不对称网格中心不在原点检查坐标定义N取奇数坐标用linspace(-L/2,L/2,N)中心光强单调变化菲涅尔数范围不对手算N_F调整z范围确保N_F跨越1边缘出现异常亮斑采样不足减小dx1增加N1或减小L1FFT法结果与直接法差异大周期性边界效应检查观察平面尺寸增大FFT点数或加窗计算速度极慢循环过多用profiler分析矩阵化、GPU加速、轴对称简化光强归一化后最大值不在中心相位计算错误检查k*r项确保r包含z²相位符号正确5.5 几个容易被忽略的物理细节第一个是倾斜因子。很多简化教程直接令K1这在傍轴近似下问题不大但如果观察平面边缘的观察角超过10度误差就不可忽略了。我的做法是始终保留K(1cosθ)/2计算量增加很少但精度提升明显。第二个是孔径边缘的相位突变。实际圆孔的边缘不是理想的阶跃函数存在一定的过渡区域。在仿真中如果网格间距太大边缘的相位突变会导致数值振荡。解决办法是适当加密网格或者在孔径函数上加一个平滑窗。第三个是偏振效应。严格的矢量衍射理论需要考虑偏振但标量近似在傍轴条件下已经足够。如果你的应用场景涉及大角度衍射或者高数值孔径就需要升级到矢量模型。6. 参数扫描与结果分析从仿真数据中读出物理6.1 传播距离对衍射图样的影响固定孔径半径a1mm波长λ632.8nm让传播距离z从10cm变化到200cm观察中心光强和第一暗环半径的变化。这个扫描用直接积分法做会比较慢我建议用轴对称简化后的一维积分速度快很多。从物理上预期随着z增大菲涅尔数N_F减小衍射图样逐渐从菲涅尔区过渡到夫琅禾费区。中心光强的振荡幅度逐渐减小最终趋于稳定。第一暗环半径近似正比于z因为夫琅禾费衍射的暗环半径公式是1.22λz/(2a)。我在实际扫描中发现一个有趣的现象当z恰好使得圆孔包含奇数个半波带时中心光强达到极大值包含偶数个半波带时中心光强达到极小值。这个振荡在z较小时非常明显z增大后逐渐衰减。你可以把这个振荡曲线和理论上的半波带数公式对比验证仿真的正确性。6.2 孔径半径对衍射图样的影响固定z50cm让a从0.2mm变化到2mm。菲涅尔数N_F a²/(λz)随a²增长所以孔径越大菲涅尔数越大衍射越接近几何光学。具体表现是中心亮斑的尺寸减小旁瓣的强度减弱整体图样越来越像圆孔的几何投影。这个扫描对于理解“什么时候可以用几何光学近似”很有帮助。我的经验是当N_F 10时几何光学近似已经相当好了当N_F 0.1时必须用夫琅禾费衍射中间区域就是菲涅尔衍射的天下。6.3 波长对衍射图样的影响固定a1mmz50cm让λ从400nm变化到700nm。菲涅尔数N_F a²/(λz)随λ增大而减小所以波长越长衍射越明显。具体表现是中心亮斑变大旁瓣变强整体图样更加“发散”。这个扫描在实际应用中很有意义。比如你做激光光束分析不同波长的激光器对应的衍射图样是不同的不能混用。我在帮朋友做光学检测时就遇到过用错波长导致仿真和实验对不上的情况。7. 代码优化与工程化建议7.1 利用轴对称性加速计算圆孔衍射的轴对称性是一个巨大的优化机会。观察平面上的光场只依赖于径向坐标ρ所以可以把二维积分简化为一维积分U(ρ) (2π/(jλz)) ∫₀ᵃ exp(jk(ρ²r₀²)/(2z)) · J₀(kρr₀/z) · r₀ dr₀其中J₀是零阶贝塞尔函数。这个一维积分用MATLAB的integral或者离散求和都能快速计算。我实测过对于N2301的观察平面一维积分比二维积分快大约100倍。% 轴对称一维积分 rho linspace(0, L2/2, N2); r0 linspace(0, a, N1); dr0 r0(2) - r0(1); U_radial zeros(size(rho)); for i 1:length(rho) integrand exp(1j * k * (rho(i)^2 r0.^2) / (2*z)) .* ... besselj(0, k * rho(i) * r0 / z) .* r0; U_radial(i) sum(integrand) * dr0 * 2 * pi / (1j * lambda * z); end I_radial abs(U_radial).^2; I_radial_norm I_radial / max(I_radial);7.2 参数化函数封装把整个仿真封装成函数方便批量扫描参数function [I2, x2, y2] fresnel_circular(a, z, lambda, N1, N2, L2) % 输入参数 % a - 圆孔半径 % z - 传播距离 % lambda - 波长 % N1 - 孔径平面采样数 % N2 - 观察平面采样数 % L2 - 观察平面尺寸 % 输出 % I2 - 归一化光强分布 % x2, y2 - 观察平面坐标 k 2 * pi / lambda; L1 6 * a; x1 linspace(-L1/2, L1/2, N1); [X1, Y1] meshgrid(x1, x1); dx1 x1(2) - x1(1); aperture double(X1.^2 Y1.^2 a^2); x2 linspace(-L2/2, L2/2, N2); [X2, Y2] meshgrid(x2, x2); U2 zeros(N2, N2); x1_vec X1(:); y1_vec Y1(:); ap_vec aperture(:); for m 1:N2 for n 1:N2 r sqrt((X2(m,n) - x1_vec).^2 (Y2(m,n) - y1_vec).^2 z^2); K (1 z ./ r) / 2; integrand ap_vec .* exp(1j * k * r) ./ r .* K; U2(m,n) sum(integrand) * dx1^2 / (1j * lambda); end end I2 abs(U2).^2; I2 I2 / max(I2(:)); end封装成函数之后参数扫描就变成了简单的循环调用代码整洁很多。7.3 结果保存与报告生成做课程设计或者项目报告的时候通常需要把多组参数的结果整理成对比图。我习惯用subplot把不同z值的结果画在一起或者用montage函数拼接多张图。MATLAB的saveas和print函数可以把图保存成高分辨率图片方便插入报告。% 批量扫描并保存 z_list [0.1, 0.2, 0.5, 1.0, 2.0]; figure; for i 1:length(z_list) [I2, x2, ~] fresnel_circular(1e-3, z_list(i), 632.8e-9, 501, 201, 0.02); subplot(1, length(z_list), i); imagesc(x2*1e3, x2*1e3, I2); colormap(hot); axis equal tight; title(sprintf(z %.1f cm, z_list(i)*100)); end saveas(gcf, fresnel_z_scan.png);8. 从仿真到实验几个衔接要点仿真做得再好最终还是要和实验对比。这里分享几个从仿真到实验的衔接经验。第一实验中的圆孔不可能是完美的圆形边缘总有毛刺或者不圆度。如果你的仿真和实验在旁瓣区域对不上先检查圆孔的加工质量。我遇到过用针孔代替圆孔的情况结果衍射图样出现了明显的六边形对称性因为针孔不是圆的。第二实验中的光源不是理想的单色平面波。激光器有发散角光束不是严格的平面波。如果你的仿真用的是平面波入射而实验用的是高斯光束中心区域的图样会有差异。解决办法是在仿真中把入射光改成高斯光束只需要在孔径函数上乘以一个高斯包络。第三实验中的观察屏有颗粒噪声CCD相机的像素响应也不完全线性。对比仿真和实验时不要期望像素级完全一致主要看条纹的位置、间距和相对强度。第四传播距离的测量误差会直接影响菲涅尔数。如果你的仿真和实验在中心光强的振荡周期上对不上先检查z的测量精度。我用卷尺量z的时候误差大概在±2mm对于z50cm来说相对误差0.4%影响不大。但如果z只有10cm这个误差就不可忽略了。9. 我踩过的几个坑和对应的解决方案第一个坑是网格不够密导致的相位混叠。我一开始用N1101孔径平面尺寸取4a算出来的图样在边缘区域出现了明显的波纹。后来把N1加到501波纹就消失了。原因是孔径边缘的相位变化太快采样不足导致混叠。判断采样是否足够的简单方法把N1加倍如果结果不变说明采样够了如果结果变了说明还不够。第二个坑是观察平面尺寸设置不当。我一开始把L2设得很大想看到完整的衍射图样结果边缘区域的强度计算出现异常。原因是观察平面太大时边缘点的观察角很大傍轴近似不再成立倾斜因子的计算也变得敏感。解决办法是把L2控制在合理范围内一般不超过孔径的20倍。第三个坑是FFT法的周期性边界效应。用FFT法时如果孔径平面尺寸不够大孔径的周期性复制会导致衍射图样出现虚假的干涉条纹。解决办法是把孔径平面尺寸加大到孔径的8~10倍或者在孔径函数上加一个平滑窗。第四个坑是单位不统一。MATLAB里所有的长度单位必须统一我习惯全部用米。但有时候从文献里抄参数文献用的是毫米或者微米忘记换算就会导致结果完全错误。比如波长632.8nm写成632.8e-9m是对的写成632.8就会算出荒谬的结果。第五个坑是光强归一化的方式。有些人用I2 / sum(I2(:))做归一化这是错误的因为总能量在传播过程中是守恒的但观察平面的尺寸和孔径平面的尺寸不同总能量不能直接比较。正确的做法是用I2 / max(I2(:))做峰值归一化或者用理论上的总功率做绝对归一化。10. 这个仿真还能怎么扩展圆孔菲涅尔衍射是一个很好的起点掌握了之后可以往几个方向扩展。第一个方向是多孔衍射。把单个圆孔改成多个圆孔观察干涉和衍射的叠加效果。这个在光子晶体、超表面等领域有实际应用。代码上只需要修改孔径函数把多个圆孔的并集作为新的孔径。第二个方向是环形孔径。把实心圆孔改成环形观察中心光强的变化。环形孔径在光学系统中常用于遮拦中心光束比如反射式望远镜的副镜遮拦。第三个方向是相位型孔径。在圆孔内加入相位调制比如螺旋相位板产生涡旋光束。这个在光镊、光通信领域很热门。代码上只需要在孔径函数上乘以exp(j·m·θ)其中m是拓扑荷数。第四个方向是部分相干光衍射。把完全相干的平面波改成部分相干光观察相干度对衍射图样的影响。这个需要引入互相干函数计算量会大一些但物理上更有意思。第五个方向是矢量衍射。当孔径尺寸接近波长时标量近似不再成立需要考虑偏振效应。这个需要用到矢量衍射理论代码复杂度会显著增加但精度也更高。我个人在实际操作中的体会是圆孔菲涅尔衍射这个题目虽然经典但真正把每一个细节都搞清楚需要反复调试和验证。我建议你先用直接积分法把基本流程跑通然后用FFT法做交叉验证再用轴对称简化做快速参数扫描。三步走下来你对菲涅尔衍射的理解会深入很多。最后再分享一个小技巧把仿真结果和理论上的半波带数对比如果中心光强的振荡周期和半波带数的变化一致说明你的仿真基本靠谱了。