简介这份压缩包提供基于传输矩阵方法计算多层介质堆叠透射率与反射率光谱的Matlab程序主要面向计算机、电子信息工程、数学等专业的大学生课程设计、期末大作业与毕业设计也适合光学方向初学者快速入门。包内共两个文件包括一个M脚本代码文件与一张结果示意图压缩包大小仅25KB轻量易用。代码采用参数化编程注释详细参数可灵活修改兼容多个常用Matlab版本用户可直接运行脚本获得光谱曲线。实现中光波在每一层介质内的相位变化与振幅变化被整合为矩阵运算通过矩阵乘积即可得到整个堆叠的总体透射率与反射率计算过程清晰、便于扩展。目前已有162人学习下载适合需要快速验证理论、开展光学仿真实验或完成课程设计的学生与研究人员。借助这份程序使用者能更直观地理解介质参数对光波传播的影响掌握多层膜系光学特性的计算思路为后续光学元件设计及更高级别光电子工程问题解决提供实用工具。1. 多层介质堆叠的透射率和反射率光谱为什么绕不开传输矩阵做光学薄膜、光子晶体、滤光片或激光器反射镜的朋友大概率都下载过类似《使用传输矩阵方法的多层介质堆叠的透射率和反射率光谱Matlab代码.rar》这样的资源。解压之后里面的 m 文件往往互相调用、变量命名随性跑出来又不知道对不对。问题不在代码本身而在于传输矩阵方法TMM这套思想没理清它用一组 2×2 矩阵描述电磁波在每一层介质中的传播与界面边界条件最终把透射率、反射率光谱浓缩成几个矩阵元素。相比 FDTD 或 COMSOLTMM 不需要网格、没有时域迭代几百个波长点几毫秒就能算完特别适合膜系设计的初期筛选。这篇博文从特征矩阵推导讲到 Matlab 实现再给出调参和排错的方法目标是让读者能自己写出可复用的多层堆叠光谱计算函数而不是依赖来路不明的 rar 压缩包。2. 传输矩阵方法的特征矩阵推导从单层到堆叠2.1 光学导纳与单层特征矩阵传输矩阵的核心是把每一层介质抽象成一个 2×2 矩阵连接层两侧的切向电场和切向磁场。设介质折射率为 (n)入射角在层内为 (\theta)定义光学导纳[ \eta \begin{cases} n\cos\theta \text{TE 偏振} \ n / \cos\theta \text{TM 偏振} \end{cases} ]这个量实际上就是该层中切向磁场与切向电场的比值。对无磁损耗、无铁磁性的介质这个公式足够精确。层内传播相位为[ \delta \frac{2\pi n d\cos\theta}{\lambda} ]其中 (d) 是物理厚度(\lambda) 是真空波长。单层介质的特征矩阵写成[ M \begin{bmatrix} \cos\delta i\sin\delta / \eta \ i\eta\sin\delta \cos\delta \end{bmatrix} ]为什么是这种形式把电磁场在层内上下界面的切向分量写成行波叠加经过指数函数的和差化简自然得到 (\cos\delta) 和 (\sin\delta) 的组合。记住这个矩阵后面所有代码都建立在这两行公式上。2.2 多层堆叠与总反射透射系数多层介质的意思就是入射介质与基底之间夹着若干层薄膜。把每一层的 (M_i) 按光线传播顺序从左到右连乘得到膜系总特征矩阵 (M_{\text{total}} M_1 M_2 \cdots M_N)。接着定义[ \begin{bmatrix} B \ C \end{bmatrix} M_1 M_2 \cdots M_N \begin{bmatrix} 1 \ \eta_s \end{bmatrix} ]其中 (\eta_s) 是基底的导纳。这样处理之后膜系连同基底被等效成一个导纳 (Y C/B)。反射系数和透射系数分别为[ r \frac{\eta_0 B - C}{\eta_0 B C}, \qquad t \frac{2\eta_0}{\eta_0 B C} ]这里的 (\eta_0) 是入射介质的导纳。反射率和透射率取模方[ R |r|^2, \qquad T \frac{\operatorname{Re}(\eta_s)}{\eta_0} |t|^2 ]如果基底无吸收(\operatorname{Re}(\eta_s) \eta_s)。这个递推式的顺序不能反入射侧的先乘否则边界条件会错。式中的 (B) 和 (C) 虽然在很多教材里被称作“等效参数”但实际求解时只需按顺序做矩阵乘不需要显式求逆。2.3 斜入射与偏振的导纳修正上面所有公式都只依赖于导纳 (\eta)所以处理斜入射时不需要改动矩阵结构只需要把每层内的折射角算出来。根据 Snell 定律[ n_i \sin\theta_i n_0 \sin\theta_0 ]于是[ \cos\theta_i \sqrt{1 - \left(\frac{n_0 \sin\theta_0}{n_i}\right)^2} ]当介质有损耗时折射率是复数开方后的 (\cos\theta_i) 也是复数它决定了波在层内的衰减方向。TE 和 TM 的差异全在导纳公式里TE 用 (\eta_i n_i\cos\theta_i)TM 用 (\eta_i n_i/\cos\theta_i)。入射介质的导纳同样按这个规则取。理解了这一层后面写 Matlab 代码时就不会把偏振条件写错。3. Matlab 代码实现建立可复用的透射率和反射率函数3.1 函数接口设计与单位约定写传输矩阵代码的第一原则接口要能直接替代“查表”。输入是波长数组、入射角、偏振、各层折射率和厚度输出是相同长度的透射率、反射率、吸收率数组。厚度和波长必须同单位建议统一用纳米因为可见光波长在 380~780 nm膜厚也常用几十到几百纳米数值不会出现 10 的负九次方这种容易出错的形式。下面的参数表对应一个完整函数调用写代码时保持参数顺序稳定尤其是nLayer和dLayer必须一一对应且从入射介质一侧开始排列。参数含义类型说明lambda波长数组double与厚度同单位建议 nmtheta0入射角double单位度标量pol偏振charTE 或 TM支持 s/pn0入射介质折射率double一般空气取 1nLayer各层复折射率向量double 向量按入射侧到出射侧排列dLayer各层厚度向量double 向量与波长同单位nSub基底折射率double如玻璃 1.523.2 核心函数 tmm_multilayer这里给出一个简洁版本直接对波长数组做向量化运算不需要循环扫描每个波长点function [T, R, A] tmm_multilayer(lambda, theta0, pol, n0, nLayer, dLayer, nSub) % TMM_MULTILAYER 传输矩阵法计算多层介质堆叠的透射率和反射率光谱 % 输入 lambda 为行向量厚度 dLayer 与 lambda 同单位。 % 输出 R、T 为与 lambda 等长的数组A 1 - R - T。 th deg2rad(theta0); sin_th n0 * sin(th); % Snell 不变量 % 基底导纳 cosSub sqrt(1 - (sin_th ./ nSub).^2); if strcmpi(pol, TE) || strcmpi(pol, s) eta0 n0 * cos(th); etaSub nSub .* cosSub; else eta0 n0 / cos(th); etaSub nSub ./ cosSub; end B ones(size(lambda)); C etaSub; for k 1:numel(nLayer) nk nLayer(k); dk dLayer(k); cosk sqrt(1 - (sin_th ./ nk).^2); if strcmpi(pol, TE) || strcmpi(pol, s) etak nk .* cosk; else etak nk ./ cosk; end delta 2 * pi * nk .* cosk .* dk ./ lambda; cp cos(delta); sp sin(delta); newB cp .* B (1i * sp ./ etak) .* C; newC 1i * etak .* sp .* B cp .* C; B newB; C newC; end r (eta0 * B - C) ./ (eta0 * B C); t 2 * eta0 ./ (eta0 * B C); R abs(r).^2; T real(etaSub ./ eta0) .* abs(t).^2; A 1 - R - T; end这段代码的巧妙之处在于没有构造二维矩阵而是把 2×2 矩阵乘法展开成B和C两个数组的递推。B初始化为全 1C初始化为基底导纳在每一层里计算新的newB和newC对应特征矩阵左乘当前列向量。因为所有变量都是数组一次调用就能算出整个波长范围的光谱。3.3 代码中的向量化逻辑说明很多人习惯写成双层for循环先循环波长再循环层数。对于 1000 个波长点和 20 层膜那就要执行两万次矩阵运算Matlab 的解释型循环会让速度慢一个数量级。这里的向量化思路是把波长当成数组维度矩阵乘法中用到的cos(delta)、sin(delta)都变成数组乘法和加法按元素自动广播。层数只有几十层外层循环次数不多内层全部是数组运算实际测试中 1000 个波长点加 20 层膜耗时通常在 10 毫秒以内。需要注意eta0是标量而B、C是数组Matlab 会自动把标量广播到数组所以eta0 * B没有任何歧义。r和t都是复数数组取模方即可得到能量比值。这套写法在 Matlab R2016b 之后都支持隐式扩展老版本可能需要手动补.*或repmat但思路不变。4. 光谱扫描与多层实例把代码跑出可用的 R(T) 曲线4.1 实例一单层增透膜验证程序正确性写代码第一件事不是直接算复杂膜系而是先拿单层膜验证公式。取入射介质为空气 (n_01)基底玻璃 (n_s1.52)单层 MgF2折射率 (n1.38)。在正入射和 532 nm 波长下理论最优厚度满足 (nd \lambda/4)lambda linspace(400, 700, 301); % nm nLayer 1.38; dLayer 532 / (4 * 1.38); % 约 96.4 nm [T, R, A] tmm_multilayer(lambda, 0, TE, 1, nLayer, dLayer, 1.52); plot(lambda, R, lambda, T);单层增透膜在中心波长处的反射率接近零因为 (n \sqrt{n_0 n_s} \approx 1.23) 时完全匹配MgF2 的 1.38 已经能显著压制反射。如果输出在 532 nm 附近的反射率低于 0.02说明特征矩阵公式、导纳定义和厚度换算都没有问题。这里R T A 1自动满足A几乎为 0因为所有介质无吸收。4.2 实例二分布式布拉格反射镜DBR的高反带DBR 是最典型的多层介质堆叠高低折射率交替每层光学厚度都是四分之一波长。设计中心波长 550 nm高折射率层用 TiO2 取 (n_H2.35)低折射率层用 SiO2 取 (n_L1.46)基底玻璃12 层lambda linspace(450, 700, 501); nH 2.35; nL 1.46; lambda0 550; nLayer repmat([nH, nL], 1, 6); % 12 层 dLayer repmat([lambda0/(4*nH), lambda0/(4*nL)], 1, 6); [T, R] tmm_multilayer(lambda, 0, TE, 1, nLayer, dLayer, 1.52); plot(lambda, R, LineWidth, 1.5);高反带中心应该落在 550 nm 附近峰值反射率随层数增加迅速逼近 1。12 层时带内反射率通常在 0.99 以上带宽约为中心波长的 10% 到 20%。如果发现反射峰偏移第一检查各层的n*d是否等于 137.5 nm也就是 (550/4)第二检查折射率是否用了设计波长而非色散实部。4.3 入射角、偏振与厚度数组的对应关系斜入射时特征矩阵里的相位项不再等于 (2\pi nd/\lambda)而是多了 (\cos\theta_i) 因子。相同物理厚度下斜入射等效于光学厚度变小所以 DBR 反射带会蓝移TM 偏振的蓝移程度比 TE 更明显。这一现象在代码里不需要额外处理只要把theta0从 0 改为 30 度再观察结果。厚度数组dLayer始终是物理厚度不要手算角度修正那是公式内部的事情。多层堆叠的顺序也很容易踩坑。nLayer(1)对应靠近入射介质的那一层nLayer(end)靠近基底。如果下载的代码采用相反顺序结果在对称膜系里可能看不出差别但遇到非对称结构就会完全错误。建议在函数开头加一行断言assert(numel(nLayer) numel(dLayer))至少能挡住长度不一致的低级错误。5. 偏振、金属损耗与数值稳定传输矩阵代码里的 4 个坑5.1 导纳公式反转斜入射结果错得离谱TE 和 TM 的导纳经常被混淆尤其在代码里用if分支切换时一个字母写反就会导致高角度反射率计算错误。正入射时 (\cos\theta1)两种偏振完全相同所以很多人在 0 度角下方验证不出问题一旦斜入射反射率随角度变化趋势就会反掉。要记住的判断方法是TE 的电场垂直入射面磁场与界面切向分量关系直接导纳是 (n\cos\theta)TM 的磁场垂直入射面导纳要除以 (\cos\theta)。如果看到 45 度角计算结果中 TM 反射率高于 TE基本可以断定反了。5.2 复数折射率下的开方分支选择金属和损耗介质折射率是复数例如金在 600 nm 处约为 (0.25 3.3i)。计算层内角度时cosk sqrt(1 - (sin_th ./ nk).^2);开方结果可能有两个分支物理上要求波在损耗介质中沿传播方向衰减对应 (\cos\theta_i) 的虚部为正。Matlab 的sqrt对复数选择实部为正的分支在大多数情况下已经满足要求但当折射率虚部很大时表达式 (1 - (\dots)^2) 可能处于副分支附近偶尔会出现不连续。稳妥做法是计算后检查虚部符号if imag(cosk) 0 cosk -cosk; end这句话只对复数cosk起作用无损耗介质中cosk是实数imag为零不会误伤。5.3 厚度与波长单位不一致导致相位错误特征矩阵中 (\delta 2\pi n d \cos\theta / \lambda)只要 (d) 和 (\lambda) 同单位结果就正确。最常见的错误是波长用纳米、厚度用微米算出来的相位偏大三到四个数量级delta变成几十甚至上百弧度正弦余弦快速振荡反射率曲线像噪声。排查方法很简单在某个波长点手动打印delta。对于四分之一波长层在中心波长处delta应该等于 (1.5708) 附近也就是 (\pi/2)。如果打印结果是 1570那一定单位不匹配。5.4 厚金属层下的指数溢出与散射矩阵替代传输矩阵中的 (\cos\delta) 和 (\sin\delta) 在 (\delta) 为纯虚数时等价于 (\cosh) 和 (\sinh)会随金属厚度指数增长。比如 200 nm 厚的金膜在可见光波段cos(delta)可能达到 (10^{20}) 量级导致B、C溢出或输出NaN。这是特征矩阵法的固有缺点不是代码 bug。解决手段有两种一是把厚金属层拆成多个薄层每个子层厚度小于趋肤深度这种方法治标不治本但代码改动最小二是改用散射矩阵法 SMM数值稳定且同样基于传输思想。如果只是算常规光学膜系金属层厚度通常不超过 100 nm特征矩阵法仍然够用。遇到溢出时优先检查层厚和折射率虚部不要先怀疑矩阵公式。6. 批量扫角度并接上优化多层膜光谱计算的进阶用法6.1 一次扫上百个入射角在做减反膜或滤光片设计时角度响应和波长响应同样重要。最简单的批量扫角度是把theta0放进循环对每个角度调用一次函数然后拼成二维矩阵thetaList 0:2:60; Rmap zeros(numel(lambda), numel(thetaList)); for j 1:numel(thetaList) [T, R] tmm_multilayer(lambda, thetaList(j), TM, ... 1, nLayer, dLayer, 1.52); Rmap(:, j) R; end surf(thetaList, lambda, Rmap, EdgeColor, none);得到的就是角度-波长反射率图谱可以直观看到 TM 偏振的高反带蓝移和 p 偏振 Brewster 角附近的反射谷。注意循环内不要重复分配大数组提前初始化Rmap能让速度提升一倍以上。若追求极致可以把函数内部继续向量化到角度维度但多数设计场景下几十个角度的循环已经足够快。6.2 用 Matlab 优化工具箱反演膜厚自动寻优是传输矩阵代码最有价值的地方。比如要设计两层减反膜让 450~700 nm 平均反射率最低可以写一个目标函数function meanR designObj(x) dLayer [x(1), x(2)]; nLayer [2.05, 1.46]; % 例如 Ta2O5 / SiO2 [T, R] tmm_multilayer(... lambda, 0, TE, 1, nLayer, dLayer, 1.52); meanR mean(R); end x0 [50, 90]; opt optimset(Display, iter, TolFun, 1e-4); xOpt fminsearch(designObj, x0, opt); dH xOpt(1); dL xOpt(2);fminsearch是 Matlab 优化工具箱或基础模块里的无导数算法适合这种目标函数不平滑、解析梯度拿不到的膜厚问题。如果想加约束比如总厚度不能超过 500 nm改用fmincon并设置线性约束矩阵。目标函数每次调用都会执行一次全光谱 TMM 计算几百次迭代在毫秒级单次耗时下通常几十秒内收敛。6.3 把能量守恒当验证标准最后留一个实用技巧调用函数后先画一下A 1 - R - T。无吸收膜系中A应该在 (10^{-14}) 量级是纯粹的浮点误差有吸收层时A大于零代表膜系吸收物理上必定小于 1。如果A出现负值或超过 1说明某个波长点的数值已经溢出或导纳分支选错。这个检查放在所有优化循环之前执行一次能避免带着 bug 去求解参数。传输矩阵方法的输出没有实验噪声结果一致性就是程序正确性的最后一道保险。本文还有配套的精品资源点击获取