简介这份MATLAB光学仿真资源面向光电信息、物理及相关专业学习者生动还原了牛顿环这一等厚干涉现象的动态形成过程。通过精确计算光程差并控制厚度变化项目以动图形式呈现干涉环级数与明暗条纹随膜厚增加或减薄而连续演变帮助理解曲率半径、折射率测量等光学测量原理。资源包内共4个文件包含可直接运行的MATLAB脚本、动态演示GIF、等倾干涉与牛顿环对比分析的Word文档以及含使用说明的TXT文件压缩包整体约362MB。目前已有1794人学习下载。脚本可直接运行支持调整初始厚度、变化速率等参数以对比不同条件下的干涉图样适合课堂演示、实验教学与自学复盘对比文档则清晰梳理了两类干涉在原理和应用上的差异是光学数值模拟与应用研究的实用案例。1. 光学仿真的第一道坎牛顿环不是几个同心圆那么简单光学仿真的第一个作业多半是从牛顿环开始的。用 MATLAB 做动态仿真把牛顿环的干涉图像做成一张动图表面上是为了交作业时比静态图好看实际上它逼你把干涉条件从公式翻译成像素级的相位计算——这是后面做波前仿真、干涉测量、薄膜设计都要过的坎。我见过不少人第一版就把中心那个黑斑做成白斑不是代码写错是半波损失压根没进公式。这篇文章从光程差推导开始给出可复现的 MATLAB 代码把曲率半径、波长这些参数逐个变成动画帧最后导出 GIF 动图。适合正在做物理光学课程设计的人、想验证干涉理论的工程师以及所有被「环为什么会动」困扰的初学者。2. 先把干涉条件翻译成像素相位牛顿环的离散化建模2.1 光程差与半波损失为什么中心一定是暗斑牛顿环属于等厚干涉。平凸透镜放在平板玻璃上中间夹着一层空气膜。以接触点为中心半径 r 处的空气膜厚度 t 可以近似为 r² / (2R)其中 R 是透镜曲率半径。这个近似只在 r 远小于 R 时成立——实验中 R 通常是 0.5 m 到几米r 只有几毫米误差在 1e-8 量级完全够用。反射光的干涉必须算上半波损失。光从空气射向玻璃在玻璃表面反射时相位突变 π相当于多走了半个波长。总光程差写成 Δ 2t λ/2空气折射率近似为 1。暗环条件 Δ mλ代入得 2t (m - 1/2)λ亮环条件 Δ (m 1/2)λ代入得 2t mλ。把 t r²/(2R) 代进去就得到暗环半径和亮环半径的解析式暗环半径r_k sqrt((k - 1/2)λR)k 1, 2, 3...亮环半径r_k sqrt(kλR)当 r 0 时光程差只剩 λ/2反射光干涉相消所以牛顿环中心永远是暗的。这是判断仿真结果对不对的第一个锚点。后面写代码时会用 sin²(π r²/(Rλ)) 表示反射光强度这个函数在 r0 处取 0天然满足中心暗斑。你要是用了 cos²中心就变成白斑典型翻车现场。2.2 meshgrid 构建二维场并计算相位连续公式翻译成 MATLAB 代码核心就两步用 meshgrid 生成二维坐标网格然后对整个矩阵批量计算相位。不要用双重 for 循环逐点算MATLAB 矩阵运算一次搞定。lambda 632.8e-9; % 波长单位米氦氖激光常见波长 R 1.0; % 透镜曲率半径单位米 L 8e-3; % 模拟区域边长单位米 N 2000; % 每边采样点数 x linspace(-L/2, L/2, N); [X, Y] meshgrid(x); % 生成 N×N 二维坐标网格 r2 X.^2 Y.^2; % 半径平方矩阵 phase pi * r2 / (R * lambda); % 反射光相位差的一半 I sin(phase).^2; % 反射光干涉强度中心为 0meshgrid(x) 相当于 meshgrid(x, x)生成两个 N×N 矩阵 X 和 Y。r2 是每个像素到中心接触点的距离平方。这里最容易被忽略的是坐标单位linspace 生成的 x 是米lambda 也是米如果混用毫米和纳米条纹密度会差几千倍。习惯上我在变量名后面直接注释单位避免算完发现图像是一张纯色板。2.3 采样率与像素标定让理论环半径可以直接比对生成图像之前先算一组采样约束。相邻条纹间距在半径 r 处大约等于 Rλ/(2r)越靠外沿条纹越密。模拟区域边长为 L边缘处半径 r ≈ L/2最小条纹间距 d_min ≈ Rλ/L。MATLAB 图像每个像素对应物理尺寸 dp L/N要让条纹看得清至少满足 dp ≤ d_min / 2也就是 N ≥ 2L²/(Rλ)。拿上面参数算L 8e-3R 1λ 632.8e-9得到 N ≥ 202。我取了 N2000采样余量很足。如果 R 小到 0.2 mL 又开到 20 mm同样公式要求 N 接近 100002000 就会看到摩尔纹。这个检查必须写在参数前面而不是等图像出问题再回头调。像素标定也很简单每像素对应的物理尺寸 dp L/N。暗环第 k 个的理论像素半径 sqrt((k - 1/2)λR) / dp。用图像工具读出环的像素位置跟这个值比对误差在 1~2 像素内基本正常超过 5 像素就该查坐标单位。参数取值作用单位陷阱lambda632.8e-9决定条纹疏密必须转成米R0.5~2.0控制环半径整体尺度和 r 保持一致的单位L8e-3模拟视野边长过大会导致条纹过密N2000网格密度受 N ≥ 2L²/(Rλ) 约束3. 静态牛顿环最小代码20 行生成第一幅干涉图3.1 单帧生成代码与运行顺序前面已经把相位矩阵算好了静态图只是加一个显示层。完整脚本如下直接复制到脚本文件或实时编辑器都能跑close all; lambda 632.8e-9; % 氦氖激光波长单位米 R 1.0; % 曲率半径单位米 L 8e-3; % 模拟区域边长单位米 N 2000; % 每边采样点数 x linspace(-L/2, L/2, N); [X, Y] meshgrid(x); r2 X.^2 Y.^2; I sin(pi * r2 / (R * lambda)).^2; imagesc(x*1e3, x*1e3, I); % 坐标轴转成毫米显示 axis image; colormap(gray); title(牛顿环反射光干涉); xlabel(x (mm)); ylabel(y (mm));运行顺序上先算 r2再算 I最后 imagesc。这里有个细节imagesc 的 x 轴和 y 轴标签我用毫米显示但 I 矩阵本身对应的物理坐标仍然是米只是显示层做了缩放。axis image 保证像素是正方形不会把圆环拉成椭圆。3.2 imagesc 与 imshow 的方向、坐标与色阶差异静态图出来了接下来要留意显示函数的方向问题。imagesc(x, y, C) 显示时 C 的每一列对应 x 坐标每一行对应 y 坐标。对牛顿环这种中心对称图像转不转置肉眼看不出来但如果你想叠加一条半径方向的剖面线坐标对不上就会很痛苦。我一般会显式转置一次imagesc(x1e3, x1e3, I.)确保后续用 xline 或 plot 叠加时方向一致。imshow 也可以显示但它不做坐标映射像素位置直接对应矩阵下标。如果代码后面要读取暗环像素位置imshow 反而方便因为它省掉了一次坐标换算。我的习惯是出图给别人看用 imagesc自己做验证用 imshow然后直接在矩阵下标层面操作。色阶方面gray colormap 是 64 级灰度干涉条纹包络很光滑。不要用 jet 或 parula会把干涉图染成伪彩色掩盖细节。如果想让暗环更锐利可以把 colormap 改为 gray(256)并配合 caxis 把强度范围压缩到 [0, 1]。3.3 反射光与透射光互补一份代码出两幅图反射光和透射光强度是互补的。反射光用 sin²透射光用 cos²两者之和恒等于 1。这可以做一个很直观的对照实验同一组 R、lambda把 cos² 的结果画在旁边你会看到反射光的暗环位置正是透射光的亮环位置。I_ref sin(pi * r2 / (R * lambda)).^2; I_tra cos(pi * r2 / (R * lambda)).^2; subplot(1, 2, 1); imagesc(x*1e3, x*1e3, I_ref); axis image; colormap(gray); title(反射光); subplot(1, 2, 2); imagesc(x*1e3, x*1e3, I_tra); axis image; colormap(gray); title(透射光);注意中心位置反射光中心是暗的透射光中心是亮的。这个互补关系本身就是检查代码有没有写反的试金石。如果两幅图中心同为暗斑那多半是半波损失处理出了问题。4. 动态仿真牛顿环扫描曲率半径并导出 GIF 动图4.1 扫描参数怎么选R、λ、膜厚三个方向的视觉差异动态仿真的本质是让干涉条件里的某个物理量随时间变化。常见做法是扫描曲率半径 R、扫描波长 λ、或者给相位加一个随时间变化的倾斜项。三种扫描的视觉效果完全不同扫描参数视觉表现背后物理曲率半径 R 增大环向外扩散条纹变疏r_k ∝ sqrt(R)波长 λ 增大环向外扩散条纹变疏r_k ∝ sqrt(λ)相位加线性倾斜项环整体平移中心点移动等效于楔角变化我最常用的是扫描曲率半径 R。原因有两个一是 R 的物理意义直观平凸透镜换一枚曲率半径更大的环就变稀二是 R 出现在分母上数值变化对条纹密度的影响比 λ 更明显动画效果更易观察。实际演示时我一般让 R 从 0.5 m 匀速变化到 2.0 m整个动画时长控制在 2 秒左右。4.2 帧循环与 GIF 写入的完整脚本动态仿真的脚本结构是预先算好不变量x、r2在循环里更新变量参数重新计算强度矩阵抓取当前画面追加写入 GIF 文件。下面这段是我常用的模板gifFile newton_rings_R_scan.gif; RStart 0.5; % 起始曲率半径 REnd 2.0; % 结束曲率半径 nFrames 100; % 总帧数 fig figure(Color, w); for i 1:nFrames alpha (i - 1) / (nFrames - 1); % 对数扫描避免前几帧相位变化太剧烈 R RStart * (REnd / RStart)^alpha; I sin(pi * r2 / (R * lambda)).^2; imagesc(x*1e3, x*1e3, I); axis image; colormap(gray); title(sprintf(R %.2f m, R)); drawnow; frame getframe(fig); [ind, cmap] rgb2ind(frame2im(frame), 256); if i 1 imwrite(ind, cmap, gifFile, gif, LoopCount, Inf, DelayTime, 0.02); else imwrite(ind, cmap, gifFile, gif, WriteMode, append, DelayTime, 0.02); end end这段代码有几个参数需要说明。nFrames 取 100配合 DelayTime0.02播放时长约 2 秒看起来流畅又不卡顿。LoopCountInf 让 GIF 无限循环。R 用对数扫描而不是线性扫描这点很关键——线性扫描时 R 从 0.5 到 2.0前几帧 R 变化比例很大条纹移动速度不均匀看起来像「先快后慢」对数扫描保证每帧相位变化比例恒定动画节奏均匀。4.3 让环「平滑扩散」而不闪烁相位步进的三个约束动态仿真的坑不在代码能不能跑而在画面会不会闪烁。闪烁的根源是相邻帧之间条纹移动距离超过半个条纹周期视觉上会看到明暗环在抖动而不是平滑扩散。要避免这个必须控制相位步进。以扫描 R 为例边缘处像素的相位 φ ≈ π (L/2)² / (Rλ)R 变化时相位变化量近似为 Δφ φ · |ΔR| / R。要让 Δφ 小于 π帧数得满足 nFrames ≥ φ_max · ln(REnd/RStart) / π。拿前面的参数算一下φ_max 约 159 radln(4) / π 约 0.44算出来至少需要 70 帧。这解释了为什么 30 帧的 GIF 会闪100 帧就基本平滑。除了帧数还有两个小习惯能减少闪烁。第一GIF 写入前用 rgb2ind 统一减色到 256 色并保持 colormap 不变不然每一帧颜色表都变视频会像霓虹灯一样跳。第二不要每帧重新调用 colormap(gray)该函数会重映射颜色表在 getframe 前加一次 drawnow 即可颜色表只在循环开始前设置一次。5. 牛顿环动态仿真避坑指南现象、原因、三行解决5.1 中心黑斑变白斑半波损失多写或少写一个 π现象仿真出的牛顿环中心是亮斑跟教材里的暗斑完全相反。即使动态演示中环的扩散方向正确这个错误依然很显眼。 原因反射光强度公式用了 cos²(π r2/(Rλ))而正确的应该是 sin²。cos² 在 r0 处取 1中心自然变亮。半波损失的本质是反射相位突变 π公式里少处理这半个周期。 解决把强度公式改回 sin(pi * r2 / (R * lambda)).^2。如果想用 cos² 写就手动给相位加 πcos(pi * r2 / (R * lambda) pi/2).^2但这属于自找麻烦直接用 sin² 最省心。5.2 满屏摩尔纹几条环都数不出来采样约束没算现象画面上全是细密的水波纹或者画出来是几条歪歪扭扭的弧形跟牛顿环完全没关系。图像分辨率明明很高却看不出环形条纹。 原因网格密度 N 低于 2L²/(Rλ) 的约束条纹最密处的空间频率超过采样率发生混叠。常见于把 L 从毫米改成厘米量级、或把 R 调小到 0.2 以下却忘了同步调整 N。 解决在参数区加一行守卫代码采样不足时直接报错dp L / N; % 像素对应的物理尺寸 dMin R * lambda / L; % 模拟区域边缘的条纹间距 if dp dMin / 2 error(采样不足请将 N 至少提高到 %d, ceil(2 * L * L / (R * lambda))); end5.3 GIF 卡顿或文件巨大颜色表和帧尺寸是元凶现象动图文件动辄几十 MB播放器打开就卡或者 GIF 文件很小但画面全是色块条纹边缘像锯齿。 原因每帧都用 rgb2ind 生成新的颜色表颜色表一变GIF 需要为每一帧保存完整映射压缩率断崖下跌。另一个原因是 getframe 抓的是整个 figure如果窗口中有标题栏、坐标轴标签等多余元素帧的尺寸会远超图像本身。 解决循环外先设定固定的 256 色灰度表循环内只做 append 写入写 GIF 前用 figure(Position) 固定窗口尺寸删掉不必要的装饰元素。实测同样参数下固定颜色表能让文件体积缩到原来的三分之一左右。5.4 环没向内收反而向外爆扫描方向与相位符号现象R 从 0.5 扫到 2.0按理环应该向外扩展、条纹变稀但画面上环却在向中心收缩而且越来越密。动态方向跟理论完全相反。 原因相位表达式里 R 放错了位置。r_k sqrt(kλR)R 增大环半径必然增大。如果代码里误写成 phase pi * r2 * R / lambda效果就是 R 增大时相位增大同一位置的条纹序数变大看起来环向中心挤。 解决检查 phase 中 R 是否在分母。同时用解析环位置做静态验证先固定 R0.5 和 R2.0 各生成一张静态图用 2.3 节的像素标定法数三条暗环的像素半径确认它们分别对应 sqrt((k-1/2)λR) 的预测值再去做动态扫描。6. 进阶用滑动条实时调参并用像素行核对暗环半径6.1 基于普通 figure 的交互式动态仿真静态 GIF 适合展示但自己调试时每改一次 R 都要重跑脚本太浪费时间。我一般用一个加滑动条的交互窗口滑一下 R 图像立刻更新。普通 figure 的 uicontrol 写法最稳兼容性也最好fig figure(MenuBar, none); ax axes(Units, normalized, Position, [0.1 0.15 0.8 0.8]); ph imagesc(ax, x*1e3, x*1e3, I); axis(ax, image); colormap(ax, gray); sld uicontrol(Style, slider, ... Min, 0.2, Max, 3.0, Value, 1.0, ... Units, normalized, Position, [0.1 0.05 0.8 0.04], ... Callback, (src, ~) updateRing(src.Value, ph)); function updateRing(R, ph) I sin(pi * r2 / (R * lambda)).^2; set(ph, CData, I); end关键是提前拿到 imagesc 的句柄 ph回调里只更新 CData不重建坐标轴这样滑动时响应无延迟。R 的滑动范围设在 0.2 到 3.0足以覆盖大多数平凸透镜的典型曲率。这个交互窗口比 GIF 更适合课堂演示学生拖一下滑块就能直观看到环疏密随 R 的变化。6.2 自检方法取一行灰度三点拟合暗环半径最后一个自检技巧从生成的强度矩阵中取中心行用 findpeaks 定位暗环的像素位置再跟理论公式比对。这个方法每次跑动态仿真前我都会先跑一遍确保公式没写歪。row I(round(N/2), :); % 取中心行的强度分布 [pks, locs] findpeaks(-row); % 找极小值即暗环位置 rPix abs(x(locs)); % 暗环的物理半径 rTheo sqrt((0.5:numel(rPix)) * lambda * R); % 理论暗环半径 table(rPix(:), rTheo, VariableNames, {measured, theoretical})findpeaks 对边界假峰会误报所以我在 rPix 里加一个条件只保留 rPix 大于 3*dp 且小于 L/2 的位置。实测中心行能数出十几个环跟理论值误差基本在 1 个像素以内。这个习惯帮我避开了很多次玄学问题——图像看着正常、但数字对不上最后发现是坐标某个地方少乘了 1e3。养成任何动态仿真都带一个解析公式校验线的习惯后这类问题基本能在五分钟内定位。希望帮到你。本文还有配套的精品资源点击获取