简介光学仿真中多光束干涉是物理光学教学与研究的重要模型。针对不同反射率下透射光条纹强度分布随参数变化的典型问题这份 Matlab 资源提供了可直接运行的仿真脚本与配套结果面向光学、物理光学及光电信息科学与工程方向的学生、教师和科研人员可用来复现教材图 4.3 的曲线并进一步观察反射率对条纹对比度、间距及相位的影响。资源共 4 个文件压缩包约 77KB主要包含一个 m 格式仿真脚本、一张 jpg 干涉条纹结果图、一份 docx 强度分布分析说明和一份 txt 说明文件脚本允许调整反射率等参数基于多光束叠加计算透射光场图像直观展示条纹形态文档则解释条纹间距、对比度变化及其物理解读。目前已有 1596 人学习下载小体积高聚焦的构成适合课程演示、课后验证或科研预研既能快速还原教科书曲线也能为多层薄膜、非均匀介质等更复杂干涉仿真提供扩展起点。1. 光学仿真入门实战用MATLAB复现多光束干涉透射光强曲线做光学仿真的人十有八九都会在某个深夜对着MATLAB窗口发愁多光束干涉的透射光强曲线到底该怎么从公式变成一张能放进报告里的图我第一次跑通这个项目时最大的感受是理论和代码之间的距离往往就隔着一行向量的写法。这个MATLAB仿真脚本做的事情很聚焦在F-P腔模型下计算不同反射率R对应的透射光强分布把物理光学里那张经典的图4.3复现出来让条纹从宽平的包络一步步变成尖锐的梳齿。适合正在学物理光学、做干涉实验预处理或者想练手MATLAB数值仿真的开发者。2. 物理模型推导F-P腔里的多光束叠加与Airy透射公式2.1 F-P腔多光束干涉的物理图像Fabry-Perot腔是产生多光束干涉最经典的装置两块反射率相同的平行平面镜中间夹着一层介质入射光进入腔体后会在两个内表面之间来回反射。每反射一次就有一小部分光从出射面透射出去。于是出射方向上不再只有一束光而是无数束振幅按几何级数衰减、相位依次递增的透射子波。这些子波在透镜焦平面上相干叠加形成多光束干涉条纹。多光束干涉和双光束干涉最大的区别在条纹锐度上。双光束干涉的强度分布是余弦形式的宽条纹而多光束干涉在反射率足够高时会出现非常尖锐的透射峰。原因在于参与叠加的光束数量发生了质变。反射率越低腔内的光往返次数越少有效参与叠加的光束越少结果越接近双光束反射率越高光束在腔内往返次数越多透射子波的有效数量越大相干叠加时对相位偏离δ2kπ的条件就越敏感。这个“敏感度”正是F-P腔能用来做高分辨率光谱分析的基础。仿真里我们最关心的就是透射光强随相位差的变化关系。相位差δ来自相邻两束透射光之间的光程差写成δ 4π n d cosθ / λ其中n是腔内介质折射率d是腔长θ是腔内折射角λ是真空波长。正入射时θ0cosθ1。当δ恰好等于2π的整数倍时所有透射子波同相叠加透射光强达到最大。2.2 透射光强公式与Airy函数推导透射光强有一个标准的复振幅叠加过程。设反射系数为r、透射系数为t反射率R|r|²透射率T|t|²。入射光振幅为E0经过第一次透射后振幅为E0·t·t后续每一束都比前一束多一次往返即多乘一个反射率R的因子同时多积累一个相位差δ。把所有透射复振幅按几何级数求和取模平方后得到I_t/I_0 T² / [(1-R)² 4R·sin²(δ/2)]这是Airy函数的标准形式。如果腔体没有吸收和散射损耗能量守恒要求RT1T(1-R)代入后T²(1-R)²公式可以进一步化简为只含R的版本I_t/I_0 1 / [1 F·sin²(δ/2)]其中F 4R/(1-R)²F叫锐度参数也常写作F-P腔的“精细度系数”。注意这个F和后面要说的Finesse精细度不是同一个量后者等于π√F/2差了一个π/2的因子读文献时很容易混淆。看两个极端就能理解曲线的性格。R0时F0分母恒为1I_t/I_0恒等于1说明没有反射就谈不上干涉入射光全部透射频谱是一条平直线。R越接近1F增长得极快比如R0.9时F360R0.98时F9800。此时只要sin²(δ/2)偏离0一点点分母就会被放大到很大透射率急剧下降只有在δ2mπ处分母回到1透射率保持峰值1。这就是“暗背景上顶着极亮尖峰”的物理来源。2.3 仿真中相位差δ的两种定义方式写代码之前必须先决定横轴怎么取值。第一种做法是把d、λ直接代入δ 4πndcosθ/λ得到的是一个极大数值。比如λ632.8nm、n1、d1mm时δ约等于2万弧度画出来的横轴范围巨大所有条纹挤在一起几乎看不见。第二种做法是把δ当作抽象自变量直接取δ ∈ [0, 4π]画出两个自由光谱范围内的透射率分布。教材插图几乎都采用第二种方式。等价的说法是把腔长工作点锁定在某个整数级次附近然后只看这个级次左右不超过一个周期的局部频谱。这样画出来的曲线保留了所有物理特征同时横轴刻度直观便于标注峰值位置。代码里我也采用这种做法物理参数lambda、n、d保留在参数区用来做量级估算和换算而不会让它们的乘积直接撑爆横轴。3. MATLAB仿真实现参数化脚本与向量化绘制3.1 参数区与相位差数组的构建整个仿真可以拆成三个部分参数区、计算区、绘图区。下面三块代码拼起来就是完整的脚本。参数区负责定义相位差范围和要扫描的反射率集合这部分决定了曲线覆盖哪些物理场景。% par_scan.m 参数区 lambda 632.8e-9; % 波长 632.8 nm n 1.0; % 腔内介质折射率空气取1 d 1e-3; % 腔长 1 mm theta 0; % 正入射腔内折射角为0 % 单次往返相位差量级估算仅用于理解物理尺度 % delta_single 4*pi*n*d*cos(theta)/lambda; % 约 2e4 rad % 以抽象相位差为横轴覆盖两个自由光谱范围 delta linspace(0, 4*pi, 4000); % 要扫描的反射率覆盖低反射到高反射 Rs [0.04, 0.5, 0.9, 0.98];参数说明delta从0取到4π是因为相邻峰之间的相位差为2π画两个周期可以让周期性特征一眼可见。4000个采样点是为了在高反射率下仍然能分辨峰的形态后面避坑章节会详细说明为什么点太少会翻车。Rs数组里4个取值很有代表性0.04模拟低反射率情形0.5是中间状态0.9进入高反射率区间0.98接近真实镀膜镜面。如果你有实验数据这里的Rs可以换成实测反射率。3.2 多反射率曲线计算与向量化写法计算区用循环遍历每个反射率内部对delta数组做向量化运算。向量化是MATLAB仿真性能和代码可读性的分水岭下面这个写法比逐点for循环快一个数量级而且逻辑和公式直接对应。% 预分配结果矩阵行对应反射率列对应相位差 I zeros(length(Rs), length(delta)); for k 1:length(Rs) R Rs(k); F 4 * R / (1 - R)^2; % Airy锐度系数 I(k, :) 1 ./ (1 F * sin(delta / 2).^2); % 归一化透射光强 end逻辑说明F 4R/(1-R)²是从Airy公式整理出来的锐度系数它把反射率映射成条纹尖锐程度。sin(delta/2).^2里delta是1×4000的行向量sin运算、平方运算都按元素进行点运算符必须写全。1 ./ (1 ...)同样是逐元素除法如果写成1 /会变成矩阵除法MATLAB会报维度错误或者算出完全错误的结果。预分配I矩阵避免循环中数组动态增长这是一个很小的习惯但数据量大时能明显感受到卡顿差异。参数说明公式里没有任何多余因子因为透射光强已经被归一化到入射光强I_0。也就是说无论反射率多少曲线的最大峰值都应该是1。这是后面检查代码是否写错的一个重要判据如果max(I(k,:))不等于1说明公式一定有问题。3.3 绘图与图4.3样式还原绘图区把四条曲线画在同一张图上横轴除以π显示为0到4正好对应教材插图的坐标习惯。figure(Color, w, Position, [100 100 800 600]); hold on; colors lines(4); % 生成4个区分度较高的颜色 for k 1:length(Rs) plot(delta / pi, I(k, :), LineWidth, 1.6, ... Color, colors(k, :), ... DisplayName, sprintf(R %.2f, Rs(k))); end xlabel(相位差 \delta / \pi); ylabel(透射光强 I_t / I_0); title(多光束干涉透射光强分布曲线); xlim([0 4]); ylim([0 1.05]); legend(Location, northeast); grid on;逻辑说明plot里传入delta/pi把横轴从弧度映射到π的倍数刻度显示0、1、2、3、4与教材图4.3完全一致。ylim上界设置成1.05而不取1是怕峰顶贴到图像边界影响观察留出5%的余量是绘图小习惯。legend的DisplayName由sprintf生成反射率数值直接出现在图例中这样四条曲线分别对应哪个参数一目了然。运行后你会看到一条接近水平的线、一条圆润的周期曲线、两条越来越尖的梳齿状曲线它们叠在一起就是图4.3的标准样式。4. 仿真结果分析精细度、对比度与条纹锐化规律4.1 不同反射率曲线的演变规律运行上面代码先从视觉上总结四条曲线的变化。R0.04时曲线几乎是一条水平线只在谷底有微弱的凹陷R0.5时波峰波谷清晰可见但峰形圆润R0.9时峰开始变得尖锐背景明显压低R0.98时曲线变成一排极细的梳齿峰与峰之间几乎贴着横轴。为了量化这些视觉印象我整理了下面的参数表全部由公式直接计算得出反射率 R谷值 I_min条纹对比度Finesse精细度半高宽 FWHM (rad)0.040.8521.170.65无谷值高于0.50.500.1119.04.40.720.900.002836129.80.1050.980.00019901155.50.020对比度定义为I_max/I_minI_min 1/(1F)。注意R0.04那一行谷值是0.852也就是说透射光强最低也有85%条纹对比度只有1.17肉眼几乎无法分辨出干涉条纹。这说明低反射率下多光束干涉退化成双光束干涉的包络甚至比双光束还要平。R0.5时谷值降到0.111对比度达到9倍已经能明显看到峰谷结构。4.2 精细度Finesse与自由光谱范围FSR表里的Finesse精细度是我们最关心的性能指标定义是Finesse π√F / 2 π√R / (1-R)。从0.04到0.98Finesse从0.65暴涨到155.5增长是指数级别的。它的物理含义可以理解为FSR与条纹半高宽的比值也就是在一个自由光谱范围内能分辨多少个峰。Finesse越高F-P腔的分辨能力越强。在曲线上可以直观读出两个量。自由光谱范围FSR对应相邻两个峰之间的相位差固定为2π与反射率无关。半高宽FWHM则是峰值一半高度的宽度随R增大急剧变窄。R0.5时FWHM约0.72radR0.9时缩到0.105radR0.98时只有0.020rad——相当于一个周期内能放下约314个半高宽。如果换成频率轴峰越窄就意味着能分辨的两条谱线频率间隔越小这是F-P干涉仪作为光谱仪的核心指标。4.3 峰值恒为1的物理直觉有一个现象特别值得停下来想一想无论反射率是0.04还是0.98透射峰的峰值始终是1。初看很反直觉——反射率这么高光都被反射回去了透射峰怎么可能还是全部透射关键在于δ2mπ时腔内所有透射子波严格同相叠加透射方向发生完全相干相长反射方向完全相消能量守恒要求所有入射能量从透射端出来。这个论断有严格的能量守恒做背书同时也是代码自检的利器。在仿真脚本里加上assert(max(I(k,:)) 0.9999 max(I(k,:)) 1.0001)如果断言失败大概率是公式里的反射率或透射率关系写错了。我每次改参数或扩展功能后都会跑一遍这个断言它能挡住很多低级错误。5. 避坑指南MATLAB光学仿真里的四个常见翻车现场5.1 曲线更像余弦波横轴范围和采样密度的问题现象第一次运行脚本画出来的曲线怎么看都像余弦波完全没有“梳齿”的感觉R0.9和R0.98的曲线几乎重合。原因横轴范围设得太大或者采样点太少。比如delta取linspace(0, 200*pi, 1000)每个峰只分配到约10个点峰值和谷值之间只有几个采样点尖锐的峰形被欠采样抹平了。反射率越高这个问题越严重0.98的峰半高宽只有0.02rad如果相位差点间距大于0.1rad连峰值都可能刚好落在采样点之间。解决先把横轴限制在[0, 4π]再用linspace(0, 4*pi, 4000)保证每个峰上有足够多的点。高反射率R0.95建议把点数提高到8000。确认公式正确后再根据实际需求扩到多周期或多自由光谱范围。5.2 峰值大于1反射率与透射率的关系写反现象R0.9时曲线峰值跑到1.1甚至1.5明显超出物理上最大的透射率1。原因这是公式里反射率和透射率关系处理出了问题。常见的错误是把透射率T单独放在分子上写成I_t/I_0 T / ((1-R)² 4R·sin²(δ/2))有些版本还会把T近似成1-R再放错位置。空气在没有吸收时T1-R代入后最终公式里R是唯一的参数不会再有任何额外的透射因子。解决用原始形式T²/[(1-R)²4R·sin²(δ/2)]一步步代入最后化简成1/[1F·sin²(δ/2)]。写代码时直接采用化简后的形式从源头上杜绝T因子混入。并在代码后面加一行max(I(k,:))的检查峰值必然落在0.9999到1.0001之间。5.3 峰位整体平移δ定义没有统一现象曲线形状完全正确但峰不在δ0、2π、4π的位置而是整体平移了一段距离比如跑到0.74π处。原因不同教材对δ的起点定义有差异。有的把δ定义为往返一次的光程差对应的相位差有的定义为两束相邻透射光的相位差还有的从第一束透射光和第三束透射光来定义。差一个π的整数倍甚至差一个系数2都会导致曲线平移或者周期翻倍。解决写代码前先固定自己的定义并把它写进注释。本文采用的δ是相邻两束透射光之间的相位差所以峰出现在δ2mπ。如果要和某本书上的插图对比先确认那本书的δ定义再决定是否需要把横轴平移或者缩放。5.4 R接近1时出现NaN除零边界没处理现象把R设成1或者非常接近1比如R1.0脚本直接报NaN或者Inf整条曲线画不出来。原因公式分母里有(1-R)²R1时分母为零F变成无穷大MATLAB返回Inf后续计算产生NaN。实际镀膜镜面反射率永远达不到1但代码不会替你判断这个边界。解决在参数区加一个保护R min(R, 0.999)或者用assert(all(R 1))让脚本在非法参数时直接报错。我在之前坑里学到的是参数区要有意识地设边界而不是等MATLAB报错后再回去猜原因。6. 二次开发从静态曲线走向交互式F-P仿真面板6.1 把R变成滑块实时刷新静态脚本只能看预设的4条曲线改成交互式后就能拖动滑块观察反射率连续变化对条纹的影响。用uifigure和uislider搭建一个轻量面板核心回调函数只有十几行fig uifigure(Name, 多光束干涉透射仿真); g uigridlayout(fig, [2 1]); ax uiaxes(g); sld uislider(g, Limits, [0.01 0.99], Value, 0.90); sld.ValueChangingFcn (sld, ~) redraw(ax, sld.Value); function redraw(ax, R) delta linspace(0, 4*pi, 8000); F 4 * R / (1 - R)^2; It 1 ./ (1 F * sin(delta / 2).^2); plot(ax, delta/pi, It, LineWidth, 1.6); ylim(ax, [0 1.05]); xlabel(ax, 相位差 \delta / \pi); ylabel(ax, 透射光强 I_t / I_0); title(ax, sprintf(R %.3f Finesse %.1f, R, pi*sqrt(R)/(1-R))); end这个回调每触发一次就重算一整条曲线8000个点的向量化运算耗时不到10毫秒拖动时完全感觉不到卡顿。标题栏实时显示当前反射率对应的Finesse把“调参”和“看指标”同步起来适合当课堂演示工具。6.2 把透射和反射曲线放在同一坐标系无吸收时反射光强I_r 1 - I_t可以直接复用透射曲线Ir F * sin(delta/2).^2 ./ (1 F * sin(delta/2).^2); hold(ax, on); plot(ax, delta/pi, It, delta/pi, Ir); legend(ax, {透射, 反射});两条曲线会在每个峰处形成“透射峰对应反射谷”的镜像关系这个可视化对理解能量守恒特别直观。反射率低于0.1时反射谷很浅透射峰接近水平线反射率接近1时反射谷几乎触底。如果加上吸收损耗alpha需要把公式改成I_tI_r 1并引入吸收率参数这就离真实镀膜仿真非常近了。这套脚本我其实是从一个更早的作业版本迭代过来的。那一次因为峰值超过1我对着屏幕怀疑了半天人生最后一步步拆开公式才发现是括号写错了位置。从那以后我每次写光学仿真代码都会强制自己先写量纲注释、再写数值计算最后把“峰值应该等于多少”这种物理判据直接写成断言放进代码里作为自检清单。希望帮到你。本文还有配套的精品资源点击获取