MATLAB牛顿环动态仿真：从等厚干涉到曲率半径反演

发布时间：2026/11/5 2:45:59

简介：这份资源围绕物理光学中的等厚干涉现象，用 MATLAB 对牛顿环做动态仿真，适合光电信息科学与工程等专业的学生、实验课教师以及需要课堂演示素材的科研人员参考。压缩包共 4 个文件，约 361.96MB，包含可调参数的 .m 仿真脚本、记录干涉条纹连续变化的 .gif 动图、一份 .docx 对比分析文档以及 .txt 使用说明；脚本负责数值计算与图形绘制，动图便于直接嵌入课件，文档则补充原理对照。使用者可通过脚本自定义初始厚度与厚度变化速率，观察空气膜变厚、变薄过程中光程差与干涉环级数的演化，结合明暗条纹宽度变化理解其与膜厚、曲率半径的对应关系，并借助文档厘清等倾干涉与牛顿环在原理和观察结果上的异同。目前已有 1798 人学习。整体结构精简，既可作为干涉章节的演示素材，也能作为光学仿真入门的练手案例。

1. 牛顿环不只是课本插图：用 MATLAB 把这圈环做成会动的仿真图

牛顿环是光学仿真里性价比最高的题材之一：一块平凸透镜扣在平板玻璃上，钠灯一照，接触点周围就浮出一圈圈同心环。它成本低、现象直观，但课本上那张静态照片恰恰丢掉了最关键的信息——膜厚只要变动半个波长，整组环纹就会向内移动一整级。用 MATLAB 做动态仿真，把中心气隙从 0 逐步调到几个波长，环纹会从边缘向中心收缩并逐个消失，物理量和画面严格对应，这比任何一张静态图都有说服力。

这份资源解决的就是这一步：从光程差公式出发，直接算出可逐帧变化的干涉图并导出动图。它适合三类人：要交光学仿真作业、又不想只贴一张图的学生；需要给学生做动态演示的实验课教师；以及想拿干涉条纹练手相位测量、面形检测的工程师。核心计算只有十几行矩阵运算，不依赖任何工具箱。

2. 物理模型先成立：等厚干涉公式怎么变成一张相位矩阵

写仿真代码之前，先把公式推导一遍，因为后面所有参数错误几乎都能追溯回这一步的某个符号或某个 π。

2.1 光程差、半波损失与暗环半径

平凸透镜的球面在近轴条件下可以写成抛物线：距顶点 r 处的空气膜厚度 e ≈ r²/(2R)，R 是透镜曲率半径。反射光的干涉发生在两个界面上——玻璃下表面与空气的界面（光疏到光密，反射有 π 相变），以及空气与平板玻璃的上界面（光密到光疏，无相变）。所以两束反射光的光程差是：

Δ = 2ne + λ/2

那个 λ/2 就是半波损失，它决定了中心点是暗是亮，也是这类仿真最常翻车的地方。代入 e 之后：

亮纹条件 Δ = kλ → r_k² = (k − 1/2)·λR/n
暗纹条件 Δ = (k + 1/2)λ → r_k² = k·λR/n，k = 0, 1, 2…

k=0 对应 r=0，也就是中心的暗斑。注意这是反射光的结论；如果看透射光，两个界面的相变刚好抵消，中心是亮的，条纹与反射光互补。教材里两种图都出现过，抄错一半就全反了。

抛物面近似什么时候不成立？严格式是 e = R − √(R² − r²)，展开后比近似式多出 r⁴/(8R³) 一项。举个数：R = 50 mm、视场半径 3 mm 时，这一项大约是 81 nm，接近 λ/7，相位上相当于小半级条纹，肉眼能看出环半径偏移。所以小曲率半径 + 大视场的组合要改用严格式，而 R = 1.5 m 这种常规配置下修正项只有纳米量级，近似式完全够用。

2.2 从 r² 到像素网格：三层矩阵的搭建顺序

整个仿真的数据流是一条直线，四步走完：

坐标矩阵：由视场边长 L 和网格边长 N 生成一维坐标向量，再用meshgrid扩成二维的 X、Y。
厚度矩阵：r2 = X.^2 + Y.^2，e = d0 + r2./(2*R)。这里直接算 r² 而不是先开方再平方，省一次sqrt，在大网格上差别不小。
相位矩阵：phi = 2*pi*n*e/lambda + pi，末尾那个+pi就是半波损失。
光强矩阵：双光束干涉且两束光振幅近似相等时，I = I₁ + I₂ + 2√(I₁I₂)·cos φ，归一化后就是cos(phi/2).^2。

最后一步值得多说一句：cos(φ/2)²和(1+cos φ)/2在数学上恒等，但前者少一次加法，而且在 φ 很大时不容易丢精度——这里的 φ 在视场边缘可以到几百弧度（条纹级数高），直接用(1+cos φ)/2会看到明显的数值抖动。

坐标约定要盯紧。常见的写法是x = linspace(-L/2, L/2, N)，这在 N 为偶数时会出现「没有像素落在中心」的问题，中心暗斑变成两个半亮的像素，看起来像一条黑线。更稳的写法是显式按像素中心生成：

dx = L / (N - 1); % 像素间距 x = (0:N-1) * dx - L/2; % 像素中心坐标，N 取奇数时中心严格落在第 (N+1)/2 点

N 取奇数（比如 801）能让中心点落在网格上，环纹对称性看起来才干净。

2.3 采样够不够：dx、最小条纹间距与混叠边界

这是决定仿真「像不像」的硬指标。条纹在径向上的间距不是常数，由 r_k = √(kλR/n) 可以推出相邻暗环的间距：

Δr ≈ λR / (2·r)

也就是说，越靠外条纹越密，最密的地方在视场边缘 r_max 处。要避免混叠，像素间距 dx 至少要比最小条纹间距小一半（奈奎斯特），工程上一般取 1/4 更保险。

按 λ = 589.3 nm、视场边长 6 mm（r_max = 3 mm）、N = 801（dx = 7.5 µm）算一组：

曲率半径 R	边缘条纹间距 Δr	2·dx	采样是否足够
0.10 m	9.8 µm	15 µm	不够，边缘出现锯齿
0.50 m	49 µm	15 µm	够，余量充足
1.50 m	147 µm	15 µm	够，可再放大视场
5.00 m	491 µm	15 µm	够，但视场内环数很少

R = 0.1 m 那一行是个真实陷阱：环数多到三百多级，看着「信息丰富」，实际上边缘已经在混叠，放大后是一片摩尔纹。解决办法要么把 N 提到 1601 以上，要么把视场缩到 2 mm 以内。我在做参数对比实验时习惯先用这张表估一遍，再决定网格数——先算后跑，比跑完发现糊了再回头改要省事得多。

3. 动态仿真实现：帧循环、渲染与双路导出

物理层跑通之后，动态仿真的工作量其实集中在两件事：怎么让帧之间的过渡平滑，以及怎么把屏幕上的东西原样录下来。

3.1 参数集中管理：一个结构体管住所有可调项

散落在脚本里的数字是后期调试的最大障碍。把所有可调量塞进一个结构体，改参数时只动一处，也方便批量扫描：

P.R = 1.5; % 平凸透镜曲率半径 (m) P.lambda = 589.3e-9; % 波长 (m)，钠黄光 P.n = 1.0; % 膜层折射率，空气取 1 P.N = 801; % 网格边长（取奇数，保证中心落点） P.L = 6e-3; % 视场边长 (m) P.nFrame = 60; % 总帧数 P.d0Span = 4; % 中心气隙扫描范围，单位：波长数

d0Span = 4的意思是中心气隙从 0 扫到 4λ。这个量决定了动图里「消失几圈环」——每变化一个 λ，中心处的干涉级就变化两级（因为光程差是 2ne），所以 4λ 大约对应中心区进出 8 级条纹，视觉上足够明显又不至于太快。

核心计算函数单独放一个文件，输入输出都是纯数据，方便后面被脚本、App 或批处理复用：

function I = nr_core(R, lambda, d0, N, L, n) %NR_CORE 计算牛顿环反射光归一化光强 % R 平凸透镜曲率半径 (m) % lambda 真空波长 (m) % d0 透镜顶点与平板之间的中心气隙 (m) % N 网格边长（像素数，建议奇数） % L 视场边长 (m) % n 膜层折射率（空气取 1.0） % 返回 I 归一化光强矩阵，0~1，0 为暗纹 dx = L / (N - 1); x = (0:N-1) * dx - L/2; [X, Y] = meshgrid(x, x); r2 = X.^2 + Y.^2; % 到顶点的距离平方 e = d0 + r2 ./ (2*R); % 空气膜厚度：抛物线近似 phi = 2*pi*n*e/lambda + pi; % 光程差相位，+pi 为下表面的半波损失 I = cos(phi/2).^2; % 等价于 sin^2(pi*n*e/lambda) end

写完先做一次自检：d0 = 0时I的中心值应该是 0（暗），往外第一个极大值出现在 r₁ = √(0.5·λR/n) ≈ 0.66 mm 处。这个数用手算就能验证，比盯着图看靠谱。

3.2 单帧渲染：为什么这里用 imagesc 而不是 imshow

imshow会自动按InitialMagnification缩放显示，实际渲染尺寸依赖窗口大小，getframe抓到的帧尺寸也就跟着窗口变——导出 GIF 时会出现每帧大小不一致、播放器抖动的问题。imagesc直接把矩阵映射到 axes 的数据坐标，配合固定的Position，帧尺寸就锁死了。

fig = figure('Color', 'w', 'Position', [300 200 560 560]); ax = axes('Parent', fig, 'Position', [0.06 0.06 0.88 0.88]); I = nr_core(1.5, 589.3e-9, 0, 801, 6e-3, 1.0); hIm = imagesc(ax, I); colormap(ax, gray(256)); axis(ax, 'image'); % 保证像素为正方形，环才是正圆 axis(ax, 'off'); caxis(ax, [0 1]); % 锁定色标范围

caxis那一行是动态场景的关键。如果不锁定，MATLAB 每帧都会按当前数据的最大最小值重新拉伸色标，结果就是环纹在暗处「呼吸」——明明物理上亮度没变，画面上却一亮一暗，看起来像数值不稳定。旧版本用caxis，较新版本推荐clim，两者等价。

axis image同样不能省。默认的 axes 会按窗口比例拉伸，矩阵是 801×801 但显示区域是 560×480，环就被压成椭圆了。

3.3 帧循环：让气隙 d0 随帧号变化

动态效果的本质就是每帧改一个参数。这里改的是中心气隙 d0，它对应实验里「缓慢抬起透镜」这个动作：

d0List = linspace(0, P.d0Span * P.lambda, P.nFrame); for k = 1:P.nFrame I = nr_core(P.R, P.lambda, d0List(k), P.N, P.L, P.n); set(hIm, 'CData', I); % 只换数据，不重建图像对象 set(get(ax, 'Title'), 'String', ... sprintf('d_0 = %.0f nm', d0List(k) * 1e9)); drawnow; end

用set(hIm,'CData',I)而不是每帧重画imagesc，性能差别很大：重画要销毁并重建整个图像对象、重设色标和坐标轴，60 帧的脚本能差出几秒；在高帧数（比如 300 帧）或者循环里还要叠加工况标注时，这个差距会直接决定脚本能不能顺畅运行。

方向问题值得单独说一句。d0 从 0 增大时，为了维持同一干涉级，所需的 r 变小，环纹向内收缩并逐个消失；反过来从大往小扫，环纹向外扩张，像水波一样从中心推出去。两种都能做，但做实验录像比对时方向必须对得上，否则对比图一眼就假。

时间采样够不够，用相位算：中心处相位是 φ = 2πn·d0/λ，相邻两帧的相位差要控制在 π/4 以内，否则人眼能看出「跳帧」。换算成步长就是 Δd0 < λ/8，也就是总扫描范围除以帧数不能超过这个值：

dStep = P.d0Span * P.lambda / P.nFrame; % 每帧气隙变化量 assert(dStep < P.lambda/8, '帧数不足，动图会出现跳帧');

3.4 导出：GIF 与 MP4 的取舍和各自参数

两条导出路线差别不小，选错了后面全是麻烦。

GIF 路线胜在通用——插进 PPT、贴进网页、发消息都能直接播，代价是 256 色调色板和文件体积。关键参数是DelayTime（单位秒）和LoopCount：

gifName = 'newton_rings.gif'; delay = 0.04; % 25 fps for k = 1:P.nFrame I = nr_core(P.R, P.lambda, d0List(k), P.N, P.L, P.n); set(hIm, 'CData', I); drawnow; f = getframe(ax); % 只抓坐标区，不录窗口边框 [A, map] = rgb2ind(frame2im(f), 256); if k == 1 imwrite(A, map, gifName, 'gif', ... 'LoopCount', Inf, 'DelayTime', delay); else imwrite(A, map, gifName, 'gif', ... 'WriteMode', 'append', 'DelayTime', delay); end end

LoopCount, Inf是无限循环，不写的话默认只播一遍。WriteMode, append从第二帧开始必须写，否则每次都是覆盖，最后只剩一帧。getframe(ax)只截取坐标区，不会把 figure 标题栏和边框录进去；但注意 axes 的 title 属于 axes 的一部分，会被录进来——如果标题文字长度随帧变化（比如数字位数不同），帧的实际裁剪尺寸就会变，GIF 播放时画布会轻微跳动。

MP4 路线质量高得多，尤其是灰阶渐变多的牛顿环，256 色量化很容易出现色带。更干净的做法是绕开屏幕抓取，直接把矩阵转成 RGB 再写进视频：

vw = VideoWriter('newton_rings.mp4', 'MPEG-4'); vw.FrameRate = 25; vw.Quality = 95; open(vw); cmap = gray(256); for k = 1:P.nFrame I = nr_core(P.R, P.lambda, d0List(k), P.N, P.L, P.n); idx = uint8(round(I * 255)) + 1; % 0~1 映射到 1~256 的索引 rgb = ind2rgb(idx, cmap); % N×N×3，不经过屏幕 writeVideo(vw, rgb); end close(vw);

ind2rgb那一步把归一化光强查表成真彩色，整个链路不涉及 figure，帧率稳定、无边框、无分辨率漂移。Quality = 95是我常用的折中值，再往上文件涨得快而肉眼看不出差别。

4. 反推与标定：从仿真环纹里把曲率半径算回来

仿真图好看只是及格线。真正能证明模型对不对的，是把图里的环纹重新量一遍，看能不能把输入参数反解出来——这是把「画图」变成「测量」的分水岭。

4.1 径向平均剖面：先把二维问题压成一维

二维图上直接找环心很麻烦，但牛顿环是旋转对称的，沿任何一条半径切一刀都能拿到完整的一维干涉剖面。更稳的做法是沿整圈做方位平均，把噪声压掉一个量级：

I = nr_core(1.5, 589.3e-9, 0, 801, 6e-3, 1.0); c = (size(I,1) + 1) / 2; % 中心像素索引 dx = 6e-3 / (size(I,1) - 1); rr = (0:floor(size(I,1)/2)-1) * dx; % 采样半径序列 (m) prof = zeros(size(rr)); theta = linspace(0, 2*pi, 720); for j = 1:numel(rr) px = c + rr(j) * cos(theta) / dx; py = c + rr(j) * sin(theta) / dx; prof(j) = mean(interp2(I, px, py, 'linear'), 'omitnan'); end

interp2默认把矩阵当作索引网格，px、py 直接写像素坐标即可，不需要额外的坐标向量。用 720 个角度点做平均，单帧耗时可以忽略，但对真实相机拍到的图（有散斑噪声）效果非常明显。

4.2 峰值提取与 r²-k 线性拟合

有了剖面，暗环就是剖面里的极小值。把剖面取负，极小变极大，直接用findpeaks：

[pkVal, pkIdx] = findpeaks(-prof, ... 'MinPeakProminence', 0.5, ... % 忽略浅坑 'MinPeakDistance', 5); % 相邻环至少隔 5 个采样点 rDark = rr(pkIdx); rDark = rDark(rDark > 0); % 剔除 r=0 的中心暗点 k = (1:numel(rDark))'; p = polyfit(k, rDark.^2, 1); % r_k^2 = k * lambda * R Rfit = p(1) / 589.3e-9; fprintf('拟合曲率半径 R = %.4f m，相对误差 %.3f%%\n', ... Rfit, abs(Rfit - 1.5) / 1.5 * 100);

MinPeakProminence取 0.5 是相对归一化光强 0~1 而言的：暗环接近 0、亮环接近 1，峰谷差接近 1，取 0.5 能把真正的暗环和数值噪声区分开。MinPeakDistance要给一个下限，因为越靠外的环越密，但再密也不会少于几个像素。

这里有个容易忽略的点：中心 r=0 那个暗点在剖面上也是一个极小值，但它对应 k=0，如果不剔除，polyfit的自变量就整体偏移一位，拟合出来的 R 会偏大十几个百分点。这个错误在图上完全看不出来——拟合直线照样很直，残差照样很小，只有 R 值对不上。

正常情况下，801 网格、视场 6 mm 的配置下，Rfit与设定值的偏差应该在 0.1% 以内。如果跑出来偏差 5% 以上，基本可以判定是采样不足（边缘环纹被混叠偏移了位置）或者中心点没落准。

4.3 和实测照片对齐的三步校准

仿真和真图对比时，直接叠图往往对不上，问题通常出在成像链而不是物理模型。按顺序做三步：

步骤	目的	操作	收敛判据
平场校正	去照明不均与暗角	采集无样品背景图，逐像素相除	校正后背景灰度标准差 < 2%
中心定位	确定 r=0 的像素坐标	对若干暗环做圆拟合或灰度重心	拟合残差 < 0.5 px
尺度标定	像素坐标转物理尺寸	用已知 R 的标准件反解 dx	反解 R 与标称值偏差 < 2%