MATLAB牛顿环仿真：从等厚干涉物理模型到图像实现

简介一份针对MATLAB光学仿真教学的doc文档聚焦牛顿环实验的计算机模拟实现。内容面向需要理解光学干涉仿真原理、掌握MATLAB图像可视化与动态仿真技巧的高校学生或科研人员从实验目的、牛顿环干涉原理推导到光强二维分布的可视化与影片动画技术均有完整阐述。文档共1个doc文件压缩包仅154KB便于快速下载查阅。目前已有202人学习使用适合作为光学实验课程的辅助参考资料。具体而言文档详细讲解了利用imshow实现灰度图像显示干涉条纹的方法并给出基于影片动画技术模拟空气膜厚度变化时条纹动态移动的代码思路能帮助读者快速搭建牛顿环仿真程序加深对干涉现象的理解。 做光学实验的同学对牛顿环应该都不陌生一块平凸透镜压在平面玻璃板上单色光一照显微镜下就能看见一圈圈明暗相间的同心圆环。原理很好懂可一旦真走进实验室调光路、寻接触点、数环号、测直径每一步都能折腾半天。相比之下用MATLAB做牛顿环仿真不仅能把原理用代码完整复现一遍还能随时改波长、改曲率半径直观看到参数和条纹之间的对应关系。这篇文章就把我从物理模型到仿真代码的实现过程完整梳理一遍适合正在写物理实验报告、准备光学课程设计或者刚接触MATLAB图像仿真的同学参考。牛顿环仿真的本质并不复杂核心就是算一个二维的光强分布函数再把它显示成图像。但要把这个分布算对、画好、还能做进一步数据提取中间有不少细节值得展开说。以下内容按我从建模到调试的实际顺序来写每个环节都会给出可运行的代码和参数选择依据。1. 牛顿环的物理模型先从干涉说起1.1 空气膜厚度与光程差的关系仿真不是凭空画几圈圆环而是要先把物理过程翻译成数学表达式。牛顿环的核心装置是曲率半径R很大的平凸透镜和一块平面玻璃板透镜凸面与玻璃板相切中间夹着一层厚度不均匀的空气薄膜。设某个位置离接触点的径向距离为r由于透镜曲率半径通常在一米以上远大于观察区域尺寸膜厚d(r)可以用抛物面近似d(r) r² / (2R)这个公式其实是从圆的方程推出来的。把透镜截面看作一个半径为R的圆圆心在接触点下方R处则凸面高度y与横向距离r满足(r)² (y - R)² R²展开后忽略高阶小量就得到y ≈ r²/(2R)。这里的y就是该点处空气膜的厚度。接下来是干涉条件。单色光垂直入射后会在空气膜的上表面透镜凸面和下表面玻璃板表面各反射一束光。关键点在于一束光在玻璃→空气界面反射另一束在空气→玻璃界面反射后者属于光疏到光密的反射会产生半波损失所以两束反射光之间额外多出了λ/2的光程差。于是总光程差写为Δ(r) 2d(r) λ/2接触点处d0Δλ/2对应干涉极小所以牛顿环中心是暗斑。这个结论和实验室观察完全一致也是你验证自己仿真是否正确的第一道关口如果仿真图中心是亮的那多半是把半波损失处理反了。1.2 从光程差到位图的桥梁有了光程差下一步就是算干涉光强。双光束干涉的强度公式是I I₀cos²(φ/2)这里的φ是两束反射光之间的相位差φ 2πΔ/λ。把Δ(r)代入并化简相位差可以写成φ 2πr²/(λR) π再做一步三角化简光强分布最终变成I(r) I₀sin²(πr²/(λR))这个式子非常简洁整个仿真的“发动机”就是它。每当πr²/(λR) kπ也就是r √(kλR)时光强为0出现暗环当sin²等于1时则是亮环。这也解释了为什么牛顿环从中心向外越排越密r正比于√k环与环之间的径向间隔并不是常数。仿真牛顿环也因此变成了“算一个二维的sin²图案”而不是真的去模拟光的传播过程。这是MATLAB做干涉仿真最常见的思路——把光强分布公式离散到像素网格上再用图像函数显示出来。搞懂这一点下面的代码就顺理成章了。2. 仿真参数设计与网格离散化2.1 五个必须想清楚的参数动手写代码前先把参数定下来。我做这版仿真用的典型参数如下参数典型取值说明波长 λ589.3 nm钠黄光实验室常用光源曲率半径 R1 m平凸透镜常见量级网格数 N1024生成N×N的像素矩阵物理尺寸 L10 mm观察区域边长中心光强 I₀1归一化处理这里最容易被忽视的就是“物理尺寸L”和“网格数N”的配合。它决定了仿真图上能看到多少个环也直接决定外圈条纹会不会混叠。比如R1m、λ589.3nm、L10mm时边缘r5mm处的最大环级数约为r²/(λR)≈42环用1024个像素采样每个条纹周期还能分到约24个像素显示效果不错。如果把L拉到20mm外圈局部条纹周期会小于单个像素尺寸这时候就会出现摩尔纹和假条纹具体解决办法放到第5章说。2.2 网格坐标为什么要从负值到正值仿真图像的中心必须对应透镜接触点因此网格坐标要把原点设在正中间。用linspace(-L/2, L/2, N)生成一维坐标向量再用meshgrid扩展成二维坐标矩阵X和Y最后通过r sqrt(X.^2 Y.^2)算出每个像素到中心的径向距离。这段代码是后续所有计算的基础写成矩阵形式会比双重for循环快上几个数量级。初学者容易踩的坑是直接把x和y的范围设置成0到L这样图像中心就会落在某个角上出来的条纹也变成四分之一个环。坐标从负到正是保证“圆心在画面正中”的最直接方式千万别图省事。2.3 半波损失到底该落在哪里很多人写仿真会对半波损失怎么处理犯迷糊。关键在于半波损失引入的额外相位不是可选项而是必须加否则中心会变成亮斑与实际实验完全相反。在具体写法上有两种等价路径要么在光程差里加λ/2得到cos²表达式要么在相位差里直接加π。两种方式最终化简结果相同但注意不能两种都加否则相位多出一倍图像就会完全错乱。我的习惯是始终回到“光程差2dλ/2”这个物理式子统一用相位差公式计算光强。这样即使之后想加反射率差异、表面缺陷项代码结构也最不容易出错。3. 核心代码实现从公式到图像3.1 网格生成与参数定义直接给一段可运行的完整代码每段都加上注释。%% 参数定义 lambda 589.3e-9; % 波长单位米 R 1.0; % 凸透镜曲率半径单位米 N 1024; % 网格行数和列数 L 10e-3; % 观察区域边长单位米 %% 网格生成 x linspace(-L/2, L/2, N); % 一维坐标从-5mm到5mm [X, Y] meshgrid(x, x); % 二维网格坐标 r sqrt(X.^2 Y.^2); % 每个像素到中心的径向距离网格生成看着简单但有一点值得注意如果你模拟的是狭长观察窗口可以把meshgrid的两个向量设置成不同长度就能得到椭圆或矩形区域内的条纹片段。这在仿真某些非标准实验装置时很实用。3.2 干涉光强计算与显示%% 光强计算 d r.^2 / (2 * R); % 空气膜厚度 delta 2 * d lambda / 2; % 光程差包含半波损失 phase 2 * pi * delta / lambda; % 反射光相位差 I0 1; % 中心入射光强 I I0 * cos(phase / 2).^2; % 牛顿环光强分布这段代码算出I之后直接用imshow显示。注意imshow对double型矩阵默认把[0,1]映射到黑到白如果I的取值范围超出这个区间画面会严重失真。所以显示参数要写成空矩阵[]让MATLAB自动拉伸显示范围%% 可视化 figure; imshow(I, []); axis on; xlabel(像素); ylabel(像素); title(牛顿环MATLAB仿真 (\lambda589.3nm, R1m)); colormap(gray);如果想在报告里更直观展示明暗环边界可以用imagesc搭配axis equal再套一个parula或jet伪彩色图。实验室里看到的是黑白条纹但伪彩色在分清各级暗环时其实更好用。3.3 画剖面图验证条纹密度变化仿真图像好看只是第一步我要做的第二件事是沿中心引一条线画出光强随半径的变化也就是径向剖面图。这个图能直观看出“越往边缘条纹越密”的规律。%% 径向剖面图 profile_row I(N/2 1, :); % 取正中间一行 radius_mm x / 1e-3; % 坐标转为毫米 figure; plot(radius_mm, profile_row, LineWidth, 1.2); xlim([0, L/2/1e-3]); % 只看中心到右边缘 xlabel(半径 (mm)); ylabel(归一化光强); title(牛顿环径向光强分布剖面);从剖面图上你能清楚看到第一个暗环大约出现在0.767mm处第二个暗环约在1.086mm第三个约在1.329mm完全符合r_k√(kλR)的推算结果。每往外一级环间距越来越窄这就是牛顿环区别于均匀条纹干涉图的关键特征。4. 从仿真到实验半径提取、噪声与参数影响4.1 自动提取暗环半径并与理论值对比仿真做到位之后下一步可以做个“虚拟实验”让程序自动识别暗环位置再和理论公式比对。MATLAB里用islocalmin函数可以直接找到局部极小值点但需要注意直接对整行数据找极小值时外圈由于振荡密集可能出现假极点最好只取中心到右边缘这一段并只保留明显的最小值位置。%% 自动提取暗环半径 profile_row I(N/2 1, :); % 中心行 idx_min find(islocalmin(profile_row)); % 找所有局部极小值 idx_min idx_min(idx_min N/2); % 只取中心右侧 r_fit abs(x(idx_min)); % 暗环对应的径向距离 %% 理论值对比 k (1:numel(r_fit)); r_theory sqrt(k * lambda * R); % 暗环半径公式 T table(r_theory / 1e-3, r_fit / 1e-3, ... VariableNames, {理论半径mm, 仿真半径mm}); disp(T);跑完这段代码表格里的数值基本一致误差就是像素离散化带来的。把这种自动提取思路迁移到真实实验照片上只需要先做标定把像素尺寸换算成实际毫米数再对灰度图做平滑和二值化处理就可以实现光斑半径的半自动测量。这也是我在课程设计里把仿真结果和实验照片联动的做法。4.2 加入噪声模拟实验环境真实实验里不可能拍出完美干净的条纹。透镜表面有灰尘、光源强度不均匀、相机有底噪这些都会叠加在图像上。仿真时可以在光强矩阵上叠一个高斯白噪声更接近实际拍摄效果noise_level 0.08; I_noise I noise_level * randn(size(I)); imshow(I_noise, []);加完噪声后再去提取暗环半径你会发现靠外圈的暗环因为对比度下降识别误差明显增大。这个现象帮我理解了一个实验技巧为什么测量牛顿环时通常取中间几级环而不是最外圈——因为外圈条纹越密受噪声影响越严重读数可靠性越低。4.3 参数扫描R和λ对条纹的影响仿真最大的优势是参数可以任意扫描。把R从0.5m逐步增加到2m保持其他参数不变你会看到环间距明显变大把λ从400nm扫到700nm同样能看到条纹疏密变化。用subplot把几种情况拼在一起对比比读公式直观太多。更进一步可以对不同R下提取的暗环半径做拟合验证r_k²与k的线性关系斜率就是λR。这是一套非常标准的数据处理流程先在仿真里跑通到实验室用真实图像做同样处理心里会踏实很多。%% 参数扫描示例 R_scan [0.5, 1.0, 1.5, 2.0]; for idx 1:4 subplot(2, 2, idx); I_tmp sin(pi * r.^2 / (lambda * R_scan(idx))).^2; imshow(I_tmp, []); title(sprintf(R%.1f m, R_scan(idx))); end5. 常见问题与排查技巧实录5.1 仿真图全黑或全白这个问题九成是归一化或数据类型问题。imshow对double型矩阵默认认为取值范围是[0,1]如果I的最大值超过1会被截断成一片白如果I包含负值又会变成一片黑。最简单的解决方法是imshow(I, [])让显示范围自动拉伸到全范围。另一个常见原因是单位不统一把λ直接写成589.3而不是589.3e-9相位差会大出好几个数量级图像变成随机噪声。每次写参数之前先确认单位是米还是毫米能省掉一大半调试时间。5.2 外圈出现摩尔纹牛顿环条纹不是均匀间隔的从内到外越来越密。当设置的物理尺寸L过大而网格数N不够时最外圈条纹的局域周期可能小于单个像素于是产生混叠画面上出现不存在的粗细交替条纹。解决思路有三个减小L、增大N、或者减小R/λ让区域内条纹总数变少。也可以对图像做高斯平滑来掩盖部分混叠但那是权宜之计根子上的问题要靠合理选参解决。5.3 环纹看起来是椭圆先检查坐标网格是否等比例。X和Y都用linspace(-L/2, L/2, N)时两个方向坐标范围相同网格本身是正方形。但显示阶段如果坐标轴比例不对比如像素纵横比不是1:1圆环就会被拉成椭圆。这时加上axis equal和pbaspect([1 1 1])问题基本能解决。另一类可能原因是用了imagesc时坐标范围设置不对或者后续做插值处理时不小心改变了数组尺寸。总之先排查显示设置再回头检查物理模型。5.4 代码运行慢如果现在还写着双重for循环去遍历每个像素赶紧改成矩阵运算。1000×1000的网格用for循环算光强MATLAB要跑好几秒而向量化写法几乎瞬间完成。我在调参阶段一般用N512确认参数没问题再升到2048做最终出图。还有一种做法是先计算半平面再用flipud或fliplr镜像拼出整张图省一半计算量。当然现代MATLAB的向量化效率已经很高这类优化更多是为了交互流畅不是必需品。5.5 中心位置不是纯黑严格来说用了sin²公式后中心值理论上是0显示出来应该是纯黑。但由于离散化的坐标并不一定能精确落在r0那个像素上中心像素可能显示为一个很暗的灰点。这属于正常现象不表示物理模型错了。如果报告需要完美图像可以在显示时对中心几个像素做特殊处理或者干脆接受这一个小瑕疵。真实实验中中心往往也不是理想暗斑因为接触点本身可能有灰尘或微小变形。最后再分享一点个人体会。做牛顿环仿真最忌讳的就是照着网上代码跑一遍、看到图就完事。我第一次做的时候也是先跑出漂亮条纹觉得大功告成可当我把仿真图放进实验报告、被问了一句“第一个暗环半径你算过吗”时才真正愣住。后来老老实实从光程差公式开始推一步步验证中心暗斑、暗环间距和理论值的关系才真正理解这个实验。建议你先用本文代码跑一遍然后去改R、改λ、加噪声、做半径提取。参数每改一次你对等厚干涉的理解就会深一层。本文还有配套的精品资源点击获取