MATLAB圆孔菲涅尔衍射仿真：从原理到代码实现与避坑指南

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 参数体系的搭建逻辑

做仿真之前，必须先确定一套合理的参数体系。圆孔菲涅尔衍射涉及的核心参数包括：

这里最关键的是菲涅尔数N_F = a²/(λz)。这个无量纲数决定了衍射所处的区域：N_F >> 1时接近几何光学，N_F ~ 1时是菲涅尔衍射的典型区域，N_F << 1时过渡到夫琅禾费衍射。我在设计参数时，会先算一下菲涅尔数，确保它落在感兴趣的范围内。比如a=1mm，λ=632.8nm，z=50cm，算出来N_F ≈ 3.16，属于典型的菲涅尔衍射区域，中心会出现明显的泊松亮斑和同心环。

3. 核心原理拆解：菲涅尔衍射到底在算什么

3.1 从惠更斯-菲涅尔原理到菲涅尔-基尔霍夫积分

惠更斯-菲涅尔原理的表述很直观：波前上每一点都可以看作一个新的子波源，这些子波源发出的球面波在空间某点相干叠加，就形成了该点的光场。菲涅尔在这个基础上引入了干涉的概念，但缺乏严格的数学基础。基尔霍夫后来用格林定理严格推导了衍射积分公式，给出了倾斜因子K(θ) = (1+cosθ)/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.8nm，z=50cm，a=1mm为例，Δ₁ ≤ 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的大小。如果嫌慢，可以把观察平面的循环也向量化，但内存消耗会急剧增加。我的建议是先用小尺寸（N2=101）调试，确认结果正确后再加大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²)，当N2=301、N1=501时，运算量大约是2.3×10¹⁰次，MATLAB里跑起来确实慢。优化方案有几个：

第一，利用圆孔的轴对称性，把二维积分简化为一维积分。对于圆孔，观察平面上的光场只依赖于径向坐标ρ = sqrt(x²+y²)，所以可以用一维的贝塞尔函数积分来代替二维积分。这样计算量降到O(N2 × N1)，速度提升几个数量级。

第二，用GPU加速。如果你的MATLAB支持Parallel Computing Toolbox，把gpuArray加到矩阵上，速度能提升10倍以上。

第三，降低观察平面的采样数。对于初步验证，N2=101就够了，确认结果正确后再加大。

5.4 常见问题速查表

5.5 几个容易被忽略的物理细节

第一个是倾斜因子。很多简化教程直接令K=1，这在傍轴近似下问题不大，但如果观察平面边缘的观察角超过10度，误差就不可忽略了。我的做法是始终保留K=(1+cosθ)/2，计算量增加很少，但精度提升明显。

第二个是孔径边缘的相位突变。实际圆孔的边缘不是理想的阶跃函数，存在一定的过渡区域。在仿真中，如果网格间距太大，边缘的相位突变会导致数值振荡。解决办法是适当加密网格，或者在孔径函数上加一个平滑窗。

第三个是偏振效应。严格的矢量衍射理论需要考虑偏振，但标量近似在傍轴条件下已经足够。如果你的应用场景涉及大角度衍射或者高数值孔径，就需要升级到矢量模型。

6. 参数扫描与结果分析：从仿真数据中读出物理

6.1 传播距离对衍射图样的影响

固定孔径半径a=1mm，波长λ=632.8nm，让传播距离z从10cm变化到200cm，观察中心光强和第一暗环半径的变化。这个扫描用直接积分法做会比较慢，我建议用轴对称简化后的一维积分，速度快很多。

从物理上预期：随着z增大，菲涅尔数N_F减小，衍射图样逐渐从菲涅尔区过渡到夫琅禾费区。中心光强的振荡幅度逐渐减小，最终趋于稳定。第一暗环半径近似正比于z，因为夫琅禾费衍射的暗环半径公式是1.22λz/(2a)。

我在实际扫描中发现一个有趣的现象：当z恰好使得圆孔包含奇数个半波带时，中心光强达到极大值；包含偶数个半波带时，中心光强达到极小值。这个振荡在z较小时非常明显，z增大后逐渐衰减。你可以把这个振荡曲线和理论上的半波带数公式对比，验证仿真的正确性。

6.2 孔径半径对衍射图样的影响

固定z=50cm，让a从0.2mm变化到2mm。菲涅尔数N_F = a²/(λz)随a²增长，所以孔径越大，菲涅尔数越大，衍射越接近几何光学。具体表现是：中心亮斑的尺寸减小，旁瓣的强度减弱，整体图样越来越像圆孔的几何投影。

这个扫描对于理解“什么时候可以用几何光学近似”很有帮助。我的经验是，当N_F > 10时，几何光学近似已经相当好了；当N_F < 0.1时，必须用夫琅禾费衍射；中间区域就是菲涅尔衍射的天下。

6.3 波长对衍射图样的影响

固定a=1mm，z=50cm，让λ从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或者离散求和都能快速计算。我实测过，对于N2=301的观察平面，一维积分比二维积分快大约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，对于z=50cm来说，相对误差0.4%，影响不大。但如果z只有10cm，这个误差就不可忽略了。

9. 我踩过的几个坑和对应的解决方案

第一个坑是网格不够密导致的相位混叠。我一开始用N1=101，孔径平面尺寸取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法做交叉验证，再用轴对称简化做快速参数扫描。三步走下来，你对菲涅尔衍射的理解会深入很多。最后再分享一个小技巧：把仿真结果和理论上的半波带数对比，如果中心光强的振荡周期和半波带数的变化一致，说明你的仿真基本靠谱了。