简介面向光学仿真与透镜设计的Matlab案例资源针对半径1mm平面波经曲率半径25mm、中心厚度3mm平凸透镜的聚焦问题基于基尔霍夫—菲涅尔衍射积分公式仿真几何焦平面上的聚焦光斑强度分布计算光斑半径并与理论值对比误差适合光学工程、物理光学课程学习者和光学设计入门者参考。包体精简共2个文件1个可直接运行的.m仿真脚本负责衍射积分计算与光斑强度绘制1份docx原理文档解释公式推导、参数设定和误差分析步骤压缩包仅18KB便于快速下载与阅读。已有893人学习使用说明其对理解透镜聚焦与波动光学仿真有切实帮助。借助该案例读者不仅能掌握基尔霍夫衍射积分在Matlab中的数值实现还能学会从仿真结果中提取光斑半径并进行理论校核为后续光学系统仿真积累经验。1. 从平凸透镜聚焦光斑说起为什么几何光学算不出艾里斑半径一块半径 1 mm、凸面曲率半径 25 mm、中心厚度 3 mm 的平凸透镜放在 632.8 nm 平面波前面几何焦平面上的光斑到底多大几何光学只能告诉你焦点在透镜后某个位置光斑半径是零。真实情况是衍射把光摊开形成艾里斑半径量级在十几微米。这份 MATLAB 光学仿真源码包用波动理论写焦平面强度分布再和 1.22λz/D 的理论值对误差适合光学设计课程设计、物理光学大作业以及想看清基尔霍夫衍射积分怎么离散化的人。先别急着按运行采样、相位符号和焦平面距离对不齐结果会差出几倍。2. 基尔霍夫-菲涅尔衍射积分与平凸透镜相位变换2.1 衍射积分从哪来菲涅尔近似什么时候能用基尔霍夫衍射公式把孔径上每一点当成次级子波源观察点场强是所有这些子波复振幅的叠加。写成旁轴形式就是菲涅尔衍射积分$$U(x,y,z)\frac{e^{ikz}}{i\lambda z}\iint U_0(\xi,\eta)\exp\left(\frac{ik}{2z}[(x-\xi)^2(y-\eta)^2]\right)d\xi d\eta$$其中 $U_0$ 是透镜出射面的复振幅$z$ 是传播距离$\lambda$ 是波长$k2\pi/\lambda$。这个近似成立的条件是 $z^3 \gg \frac{\pi}{4\lambda}[(x-\xi)^2(y-\eta)^2]^2_{\max}$。按本文参数孔径半径 1 mm$z$ 取 48 mm$\lambda632.8$ nm菲涅尔数 $N_Fa^2/(\lambda z)\approx 33$远大于 1积分适用。源码里如果直接套夫琅禾费远场公式光斑半径会算错因为 48 mm 并不满足远场条件。2.2 平凸透镜的厚度函数与相位屏平凸透镜中心厚度 $t_03$ mm凸面曲率半径 $R25$ mm径向坐标 $r$。旁轴近似下厚度分布为$$d(r)t_0-\frac{r^2}{2R}$$光穿过透镜后附加相位延迟 $\phi(r)k(n-1)d(r)$$n$ 取 1.5 时常数项不影响强度二次项可以写成 $-\frac{k r^2}{2f}$其中 $fR/(n-1)50$ mm。在 MATLAB 里构造相位屏lambda 632.8e-9; % 波长He-Ne 常用值 k 2*pi/lambda; R 25e-3; % 凸面曲率半径 t0 3e-3; % 中心厚度 n 1.5; % 玻璃折射率 D 2e-3; % 透镜孔径直径 f R/(n-1); % 薄透镜焦距 N 512; % 入射面采样点数 L 3e-3; % 入射面边长要盖住 2 mm 孔径 xi linspace(-L/2, L/2, N); dxi xi(2)-xi(1); [Xi, Yi] meshgrid(xi, xi); r2 Xi.^2 Yi.^2; aperture r2 (D/2)^2; thickness t0 - r2/(2*R); phase k*(n-1)*thickness; U0 aperture .* exp(1i*phase);这段代码先构造方形网格再用圆形孔径裁掉透镜外区域。L必须大于透镜直径否则边缘衍射会被计算窗口截断。dxi后面做积分权重用不能漏。phase前面的符号决定了波前是发散还是汇聚如果跑出来的光斑越传越大先把phase的符号反过来试。2.3 几何焦平面位置薄透镜公式与厚透镜后焦距差多少薄透镜公式给出 $f50$ mm但透镜中心厚度 3 mm 已经不能忽略。平凸透镜的后焦距 $BFL$ 近似为$$BFLf\left(1-\frac{(n-1)t_0}{nR}\right)$$代入得 $BFL50\times(1-0.5\times3/(1.5\times25))\approx48$ mm。如果脚本里直接把观察面放在 50 mm光斑会略微弥散计算半径偏大。常见做法是扫一遍 $z$取中心强度最大或二阶矩半径最小的位置作为实际焦平面。量公式数值薄透镜焦距 $f$$R/(n-1)$50 mm厚透镜后焦距 $BFL$$f(1-(n-1)t_0/(nR))$48 mm理论第一暗环半径$1.22\lambda z/D$18.5 μm$z48$ mm理论 FWHM 半径$1.03\lambda z/D$15.6 μm$z48$ mm3. MATLAB 实现从平面波到焦平面光斑强度分布3.1 网格与采样参数怎么定入射面边长 $L3$ mm采样 $N512$网格间距 $d\xi5.86$ μm每个波长约 9 个采样点满足 $\lambda/2$ 以上的采样要求。观察面不能取整个 3 mm 宽否则光斑只占几个像素。艾里斑半径约 18 μm观察半宽取 100 μm 足够。观察面采样 $M512$间距约 0.39 μm可以分辨半高全宽。下面参数表可以直接抄进脚本参数含义取值lambda波长632.8e-9 mR凸面曲率半径25e-3 mt0中心厚度3e-3 mn折射率1.5D孔径2e-3 mz传播距离48e-3 mN入射面采样512M观察面采样512obs_half观察半宽100e-6 m注意观察面半宽必须大于 3 倍理论艾里斑半径否则半高位置找不到find会返回空值。3.2 直接积分法的矩阵实现直接写双重循环在 MATLAB 里会慢到无法接受利用菲涅尔积分的可分离性把二维积分拆成两个矩阵乘法z 48e-3; % 厚透镜后焦距 M 512; obs_half 100e-6; xobs linspace(-obs_half, obs_half, M); dxobs xobs(2)-xobs(1); % 菲涅尔积分核行对应观察点列对应源点 Kx exp(1i*k/(2*z) * (xobs(:) - xi).^2) * dxi; Ky exp(1i*k/(2*z) * (xobs(:) - xi).^2) * dxi; % 复振幅Kx * U0 * Ky.注意 Ky 做转置而不是共轭转置 Uf exp(1i*k*z)/(1i*lambda*z) * (Kx * U0 * Ky.); I abs(Uf).^2; I I / max(I(:)); % 归一化方便取半高Kx和Ky的每一行是一个观察点的相位因子乘dxi相当于积分权重。Ky.是普通转置如果用会变成共轭转置y 方向相位会反号光斑会被扭曲。Uf前面的 $e^{ikz}/(i\lambda z)$ 是菲涅尔积分的全局因子虽然不影响归一化强度但做能量对比时要保留。512×512 矩阵乘法在普通笔记本上几秒内能跑完若内存吃紧把M降到 256观察半宽保持不变。3.3 FFT 快速法什么时候能用坐标怎么缩放入射面采样升到 2048 时矩阵乘法会变慢可以改用 FFT 做卷积。菲涅尔衍射本质是卷积频域传递函数为fx (-N/2:N/2-1)/(N*dxi); [FX, FY] meshgrid(fx, fx); H exp(1i*k*z)/(1i*lambda*z) .* exp(-1i*pi*lambda*z*(FX.^2FY.^2)); Uf_fft ifft2(fft2(U0) .* fft2(ifftshift(H))) * dxi^2;ifftshift(H)把零频移到数组左上角符合 MATLAB 的fft2排列。dxi^2是二维积分权重。FFT 法的输出坐标是 $x\lambda z f_x$不是原来的xobs这一点很容易搞混。如果只是看强度分布形状FFT 法够用如果要和理论半径做精确对比建议先用直接积分法在 512 网格上校准再切 FFT 法做参数扫描。角谱法也可以做类似传播但这里 $z48$ mm 远大于波长菲涅尔近似已经足够。4. 光斑半径计算与理论值对比4.1 光斑半径的三种定义光学仿真里说“光斑半径”必须带定义否则误差没法比。常见三种第一暗环半径 $r_11.22\lambda z/D$半高全宽半径 $r_{\text{FWHM}}1.03\lambda z/D$二阶矩半径 $r_{\text{rms}}$。第一暗环对旁瓣和噪声敏感二阶矩对观察窗口截断敏感。课程设计里最稳的是半高全宽因为只用到中心主瓣。下面代码从强度矩阵提取半高全宽[~, idx] max(I(:)); [cy, cx] ind2sub(size(I), idx); row I(cy, :); x_centered xobs - xobs(cx); half 0.5; left find(row(1:cx) half, 1, last); right cx - 1 find(row(cx:end) half, 1, first); if isempty(left) || isempty(right) error(半高位置未找到把 obs_half 从 100e-6 调到 200e-6); end FWHM x_centered(right) - x_centered(left); r_fwhm FWHM / 2; r_theory 1.03 * lambda * z / D; err abs(r_fwhm - r_theory) / r_theory * 100; fprintf(仿真 FWHM 半径 %.2f μm理论 %.2f μm误差 %.2f%%\n, ... r_fwhm*1e6, r_theory*1e6, err);find从中心向左、向右分别找强度降到 0.5 的位置。如果观察面太窄left或right为空直接报错比默默返回错误值好。x_centered把坐标平移到峰值位置避免峰值不在正中心时半径算偏。4.2 数值结果与误差分析按前面参数跑一组仿真 FWHM 半径约 16.2 μm理论值 15.6 μm误差约 3.8%。这个量级说明菲涅尔积分和透镜相位屏基本对齐。误差主要来自四处误差来源表现排查方法厚透镜主平面偏移焦平面位置偏 1–2 mm扫z找最小 FWHM球差边缘光线焦点更近比较孔径 1 mm 和 2 mm 的结果采样不足光斑出现高频波纹把N从 512 升到 1024相位符号反了光斑越传越大交换exp(1i*phase)的符号观察面截断半高位置找不到加大obs_half球差在平凸透镜里不可忽略。平面波从平面入射、凸面出射时边缘光线会聚点比旁轴光线更靠近透镜。如果只想要理论艾里斑对比可以把孔径缩小到 1 mm球差影响会明显下降仿真半径会更接近 $1.03\lambda z/D$。4.3 用二阶矩半径做交叉验证半高全宽只用了中心行建议再用二阶矩算一次。二阶矩半径对主瓣能量分布更敏感Xc xobs - xobs(cx); [XX, YY] meshgrid(Xc, Xc); r2_map XX.^2 YY.^2; P sum(I(:)); r_rms sqrt(sum(r2_map(:) .* I(:)) / P);如果 $r_{\text{rms}}$ 和 $r_{\text{FWHM}}$ 的比值在 1.2–1.5 之间说明主瓣形状接近艾里斑。如果比值远大于 1.5通常是有旁瓣能量漏进观察窗或者焦平面没对准。5. 进阶技巧让光学仿真结果可复现5.1 能量守恒自检菲涅尔近似下总能量应该近似守恒。用入射面和观察面的强度分别乘网格面积P_in sum(abs(U0(:)).^2) * dxi^2; P_out sum(abs(Uf(:)).^2) * dxobs^2; fprintf(能量比 %.6f\n, P_out / P_in);如果能量比偏离 1 超过 5%先检查观察面是否覆盖了全部主瓣和第一旁瓣。dxobs^2是观察面二维权重漏掉它会得到量纲错误的结果。直接积分法在观察面不够大时能量比会小于 1这是截断造成的不是积分公式错。5.2 零填充与边界效应用 FFT 法时fft2默认做循环卷积光斑跑到另一侧会叠回来。零填充可以缓解Npad 2*N; U0_pad zeros(Npad, Npad); U0_pad(1:N, 1:N) U0;把U0_pad代替U0做 FFT观察面中心区域的结果更干净。零填充不增加物理信息只是把循环卷积的边界推远。如果直接积分法结果和 FFT 法差很多优先检查ifftshift和零填充而不是改波长。5.3 直接积分与 FFT 交叉验证最后留一个校验动作在观察面上取 64×64 个点用直接积分算一遍再和 FFT 结果在相同坐标处比中心强度。如果两者偏差大于 2%优先检查H的ifftshift和dxi^2因子而不是修改透镜曲率半径。本文还有配套的精品资源点击获取