做光学镀膜的人谁还没跟反射率谱线死磕过几个晚上。我记得第一次跑这套多层膜仿真代码的时候脑子里只有一个念头能不能在写膜系之前先用电脑把光谱行为摸个底别等镀膜机开炉了才发现设计反了。答案是能而且用MATLAB就能做得挺顺手。这套光学多层膜系统模拟仿真代码核心就是用传输矩阵法计算平面波在多层薄膜里的传播行为重点处理TE模的反射率、透射率顺带可以抠出电场分布。它适合三类人刚开始接触薄膜光学的学生、需要做膜系设计预研的工程师、以及想把手头的膜系设计与实测谱线互相对照的工艺人员。看懂它你就掌握了镀膜仿真最基础也最常用的一把尺子。1. 光学多层膜仿真的核心思路为什么偏要用传输矩阵法1.1 什么是光学多层膜系统它到底在模拟什么光学多层膜说白了就是把两种或多种不同折射率的材料按设计厚度一层一层堆在基底上。常见的有相机镜头的增透膜、激光器腔镜用的高反膜、冷光镜、滤光片甚至现在光学传感器里用的共振吸收结构本质上都是同一个套路。膜层厚度通常在几十纳米到几微米之间跟可见光波长一个量级所以光进去之后会发生强烈的干涉效应。宏观上表现出来就是某些波长反射增强、某些波长透射增强通过调整每层的折射率和厚度就能“雕塑”出你想要的光谱曲线。这套MATLAB仿真代码做的就是把上述物理过程数值化。它输入的是每一层材料的复折射率、厚度、入射介质折射率、基底折射率、入射角和波长范围输出的是对应的反射率R、透射率T以及吸收率A。TE模是这套代码首先处理的偏振态物理上就是电场矢量垂直于入射面的那一束光比如很多激光系统里s偏振分量就是TE模。1.2 为什么必须用数值方法而不是手推公式单层膜的反射率公式是教科书里的经典内容很多人一开始以为多层膜也能从头到尾推一个封闭表达式出来。实际上当膜层数超过三层解析表达式的复杂度就会急速膨胀而且一旦膜层有吸收折射率是复数或入射角非零手推简直是一场灾难。传输矩阵法的思路非常聪明不直接去解每层内部复杂的多次反射叠加而是把每层等效成一个2×2的矩阵这个矩阵同时包含了光在该层中传播的相位变化和界面处的边界条件。把每一层的矩阵按顺序乘起来得到一个总的2×2矩阵然后从总矩阵元素直接读出反射系数和透射系数。这个方法的物理图像清楚、代码量少、数值稳定性好是当前光学薄膜仿真领域事实上的标准方法。MATLAB擅长矩阵运算复数、矩阵连乘、循环扫描波长这些都是基础操作所以跑起来非常顺手。1.3 为什么我推荐先把TE模单独实现TE模和TM模的差别本质上是边界条件里场分量方向不同导致的。TE模的电场只有垂直于入射面的一个分量边界条件里不涉及电场法向分量的跳变所以推导和编程都相对简洁。先把TE模跑通你就能建立对传输矩阵法的直观理解相位厚度怎么算、导纳怎么代、矩阵连乘怎么排这些核心逻辑在TM模里完全复用只是导纳表达式从 n·cosθ 换成 n/cosθ 而已。另外很多实际应用场景比如某些偏振分光膜、斜入射高反膜、光栅结构里的类导模TE模本身就是关键角色。所以先实现TE模不是偷懒而是用最小成本掌握方法再往TM、椭圆偏振、色散材料、粗糙界面等方向扩展。2. 传输矩阵法与TE模的数学原理把公式变成代码前的必要铺垫2.1 多层膜问题的“降维打击”从麦克斯韦方程组到2×2矩阵光在多层膜中的传播本质上是由麦克斯韦方程组和界面连续条件控制的。如果入射光是单色平面波每个膜层内的电场可以写成前向传播波和反向传播波的叠加。在一个界面上入射、反射、透射波的切向分量必须连续于是可以把这个约束关系表达成矩阵形式把界面一侧的场与另一侧的场联系起来。而光在一层均匀介质内部从界面a走到界面b电场只会经历一个相位变化即乘上 e 的 正负iδ 次方。这个传播过程同样可以用一个对角矩阵表示。把界面矩阵和传播矩阵合在一起就得到了一个层的特征矩阵也就是传输矩阵。[ M_j \begin{bmatrix} \cos\delta_j \dfrac{i\sin\delta_j}{\eta_j} \ i\eta_j\sin\delta_j \cos\delta_j \end{bmatrix} ]其中δ_j是第j层的相位厚度η_j是该层对应偏振态的光学导纳。整个膜系的响应就是所有层的特征矩阵按光传播方向连乘[ M M_1 M_2 \cdots M_m ]这里的矩阵乘法顺序不能乱第一层先乘最后一层最后乘这直接对应光从入射介质进入膜层再到达基底的物理路径。最后反射系数和透射系数可以从总矩阵元素与入射介质、基底的导纳组合得到。这个过程可以类比电路里的级联网络每一层是一个二端口网络整体级联后等效成一个总的二端口网络反射率和透射率就是网络的散射参数。2.2 TE模与TM模的区别为什么导纳表达式不一样导纳这个概念第一次接触会觉得抽象但你可以把它理解为“介质对电磁波的等效阻抗”由折射率和入射角共同决定。对于TE模入射面内电场只有垂直分量磁场的切向分量与电场比值关系会得到[ \eta_{TE} n\cos\theta ]这里的θ是光在该层内的折射角不是入射角。对于TM模因为涉及磁场垂直入射面、电场在入射面内边界条件里的阻抗关系不同结果是[ \eta_{TM} \frac{n}{\cos\theta} ]如果正入射θ0cosθ1TE和TM退化为一模一样。这也是为什么很多人验证代码时习惯从正入射开始这时候偏振无关公式最简单。一旦斜入射就必须分开处理。这套代码里先做了TE模后面要扩展到TM只需要把计算导纳的语句改掉其他逻辑几乎不动。2.3 相位厚度、斯涅尔定律和入射角的连锁计算相位厚度是薄膜干涉的核心量它决定了光在膜层里走一个来回积累了多少相位差[ \delta_j \frac{2\pi n_j d_j \cos\theta_j}{\lambda} ]这里n_j是第j层复折射率d_j是几何厚度λ是入射光在真空中的波长。注意每一层里的折射角θ_j并不等于入射角θ_0需要通过斯涅尔定律逐层追踪[ n_0\sin\theta_0 n_1\sin\theta_1 n_2\sin\theta_2 \cdots ]也就是说虽然光从外部以某个角度入射但进入高折射率材料后折射角会变小进入低折射率材料后会变大。每一层的光路都不一样相位厚度也不一样代码里必须循环逐层算出θ_j再代进相位厚度公式。这个逐层追角的操作是很多人写代码时容易漏掉的地方漏掉之后正入射没问题但斜入射结果就全错了。3. 核心代码模块实现从零搭一个TE模传输矩阵仿真器3.1 函数骨架与输入参数设计我用MATLAB写了一个函数输入是膜层参数和入射条件输出是反射率、透射率和吸收率。函数的好处是方便反复调用比如扫描波长、扫描角度或者扫膜厚只需要在外层套一个for循环就能批量跑不同条件。这里给出核心骨架function [R, T, A] multilayer_TE(n_list, d_list, n0, ns, lambda, theta0_deg) % multilayer_TE 传输矩阵法计算TE模多层膜反射、透射、吸收 % 输入 % n_list - 各膜层复折射率行向量长度m % d_list - 各膜层几何厚度单位nm行向量长度m % n0 - 入射介质折射率可为复数 % ns - 基底折射率可为复数 % lambda - 波长单位nm可以是标量或行向量 % theta0_deg - 入射角单位度 % 输出 % R, T, A - 反射率、透射率、吸收率与lambda同维度输入里我把厚度和波长都统一用nm这样计算相位厚度时两个量纲可以消掉不会出现单位混用导致的量级错误。折射率用复数的场景稍后会说但函数从第一版就支持复数后面处理吸收介质就不用改接口。3.2 入射角度追踪与导纳初始化主函数第一步是把入射角从度转成弧度然后初始化输出数组theta0 theta0_deg * pi / 180; N length(lambda); R zeros(1, N); T zeros(1, N); A zeros(1, N);对每个波长扫描点先算出每一层里面的折射角。注意当膜层有吸收、折射率为复数时直接套实数斯涅尔定律会有点问题这里先按常规做法取实部来算角度满足绝大多数弱吸收介质场景for idx 1:N lam lambda(idx); m length(n_list); theta zeros(1, m); % 逐层用斯涅尔定律追角 theta(1) asin(real(n0) * sin(theta0) / real(n_list(1))); for j 2:m theta(j) asin(real(n_list(j-1)) * sin(theta(j-1)) / real(n_list(j))); end % 入射介质导纳 eta0 n0 * cos(theta0);有个边界情况要小心当入射角太大或者从高折射率介质入射到低折射率介质asin 参数可能略大于1这时会发生全反射。MATLAB 会返回复数数值上不至于直接崩但结果需要你多想一想物理上是否合理。后面我会在常见问题里专门展开。3.3 特征矩阵连乘与反射率计算接下来就是核心循环首先把总矩阵设为单位阵然后逐层构建特征矩阵并连乘M_total eye(2); % 2x2单位矩阵 for j 1:m delta 2 * pi * n_list(j) * d_list(j) * cos(theta(j)) / lam; eta_j n_list(j) * cos(theta(j)); M_layer [cos(delta), 1i * sin(delta) / eta_j; 1i * eta_j * sin(delta), cos(delta)]; M_total M_total * M_layer; end这段代码是整套仿真的心脏。cos(delta)和sin(delta)里的delta直接用复数n计算所以会同时计入吸收导致的振幅衰减。矩阵元素里的1i是MATLAB的虚数单位千万别把1i当成变量i用否则会跟循环索引冲突。基底导纳要单独计算把基底看作半无限介质没有反射回来的波。然后通过总矩阵元素和入射/基底导纳组合出反射系数和透射系数% 基底内的折射角 theta_s asin(real(n0) * sin(theta0) / real(ns)); eta_s ns * cos(theta_s); % 反射系数和透射系数公式 ref_num eta0 * (M_total(1,1) M_total(1,2) * eta_s) - (M_total(2,1) M_total(2,2) * eta_s); ref_den eta0 * (M_total(1,1) M_total(1,2) * eta_s) (M_total(2,1) M_total(2,2) * eta_s); r ref_num / ref_den; trans_num 2 * eta0; t trans_num / ref_den; R(idx) abs(r)^2; T(idx) real(eta_s) / real(eta0) * abs(t)^2; A(idx) 1 - R(idx) - T(idx); end透射率公式里那个 real(eta_s)/real(eta0) 是介质阻抗修正项很多初学的人查了半天发现T和R加起来不是1就是因为漏了这层介质导纳比。有了修正项无吸收膜系里RT才严格等于1。3.4 验证代码正确性的第一步单层膜解析解对照代码写完别急着跑多层先拿单层膜正入射验证。单层膜在中心波长处的反射率存在解析解比如基底为玻璃1.52膜层为氟化镁1.38设计中心波长550nm膜厚取100nm左右正入射。用上面的函数算出来的反射率曲线应该与教科书上单层介质膜的反射率公式一致。我在实际验证时经常这么操作n_list [1.38]; d_list [99.64]; n0 1.0; ns 1.52; lambda linspace(400, 800, 401); theta0_deg 0; [R, T, A] multilayer_TE(n_list, d_list, n0, ns, lambda, theta0_deg); plot(lambda, R);这个膜厚对应的是550nm的四分之一光学厚度550/(4×1.38) ≈ 99.64nm。在550nm附近R出现一个极小值这就是增透谷。拿这个结果与理论解析式比较一致就说明矩阵构建和边界条件处理是对的。这是最快、最有效的“自检”手段胜过任何复杂膜系的直觉判断。3.5 扩展把电场分布也算出来反射率和透射率只是冰山一角很多情况下你需要知道光在膜系内部的电场分布。比如激光损伤研究里损伤经常出现在电场峰值位置设计吸收增强结构时吸收层要放到电场最大处。这个需求在代码里扩展起来不算难核心思路是得到了总矩阵M和反射系数r、透射系数t之后逐层向前递推可以求出每一层入射界面处的复振幅。function E_field calc_Efield(n_list, d_list, n0, ns, lambda, theta0_deg) % 先调用主函数得到反射系数和透射系数再逐层递推电场 % 这里给出递推的核心表达式省略了完整封装 E_plus 1; % 入射电场振幅归一化 E_minus r; % 反射电场振幅 % 在每一层入口处前向波 E_plus * exp(-i*delta_partial) % 后向波 E_minus * exp(i*delta_partial) % 总电场 前向波 后向波这块代码如果展开本质上是对每一层的传播矩阵做逐步“剥洋葱”从空气一侧往基底方向推进把每一层入口和出口的电场算出来。有了这个扩展你的仿真能力就从“只看光谱曲线”跳到“看懂场分布”这对分析干涉效应、吸收层摆放位置、镀膜应力带来的光程变化都有价值。4. 仿真结果怎么看那些藏在曲线里的信息4.1 典型高反膜堆的反射率谱线解读用上面这套代码跑一个最常见的膜系中心波长1064nm(HL)^15 高反膜堆H材料用二氧化钛n2.35L材料用二氧化硅n1.46基底是玻璃正入射。每个周期中H层厚度为1064/(4×2.35)L层厚度为1064/(4×1.46)。你会发现中心波长附近出现一个高反射带带宽大约在几百纳米量级带内反射率非常接近1。这个高反射带的物理来源是布拉格反射每个界面的反射波在中心波长附近同相叠加类似于X射线在晶体上的衍射。通过调整折射率比和周期数你能控制带宽和反射率。折射率比越大反射带越宽周期数越多带内反射率越平、越接近100%。用这个代码去扫描不同周期数的膜系你会直观看到“增加膜层数”不是无限有效膜层数超过一定数目后中心反射率变化非常小但膜层厚度和应力都在涨工程上必须做取舍。4.2 光子禁带与布拉格条件一维光子晶体在结构上跟高反膜堆完全相同区别在于你是否把它放到周期结构语境里理解。当一个周期内两种材料的相位厚度之和满足[ 2(n_H d_H \cos\theta_H n_L d_L \cos\theta_L) m\lambda ]就出现布拉格反射对应频率范围内的光无法在周期结构中传播这就是所谓的光子禁带。斜入射时θ_H和θ_L会随入射角变化禁带位置也会移动这就是很多角度选择滤光片的设计基础。你的代码里直接把入射角从0度扫到60度就能看到反射带逐渐蓝移的整个过程非常直观。4.3 电场分布揭示了什么很多人只盯反射率忽略了电场分布但在一些实际开发里电场分布比反射率更致命。曾经我做介质高反膜激光损伤实验发现膜层总是从某个特定界面开始坏一开始怀疑材料纯度后来用仿真把1064nm波长下的电场分布算出来才发现那个界面临近电场波腹电场强度是入射场的两倍以上。从那以后设计高功率激光膜系时我都会顺手把电场分布打印出来看峰值有没有落在薄膜界面附近。在λ/4高反膜堆里驻波场的波节通常落在高折射率层中心波腹落在低折射率层中心或界面附近。如果你把吸收较强的材料放到波腹位置吸收损耗会被显著放大反过来把吸收小的材料放在波腹区吸收型材料放在波节区就能实现“低吸收、高损伤阈值”的设计。这就是电场分布仿真在镀膜工艺里的价值——它让你看到光在膜系里到底“待在哪里”。5. 实操总结避坑指南与调试建议5.1 我踩过的五个坑从NaN到静默错误第一个坑是矩阵左乘右乘搞反。总矩阵必须按照光传播方向依次乘第一层在最左边。这个顺序错了单层膜结果也可能“碰巧”是对的但膜层一多就彻底乱套。第二次写代码时我就因为循环里 M_total M_total * M_layer 写反耗时大半天查不出来。第二个坑是单位不统一。厚度用微米、波长用纳米算出来的相位厚度差了1000倍反射谱看起来完全不是你以为的那个膜系。现在我的代码里给d_list和lambda都加了注释一律要求nm函数入口处再做一次检查如果发现数量级不对就提示。第三个坑是正入射时cosθ都是1很多分支逻辑不会被触发结果一切正常一旦切到斜入射各种各样的索引错误、NaN就冒出来了。这提示我正入射验证只是一个必要不充分条件斜入射至少选一个30度、60度做持续测试。第四个坑是虚数单位i和循环索引打架。MATLAB内置的虚数单位是1i如果我在某个地方顺手写了i 1这种赋值后面所有1i都会出问题。现在我一直用1i绝不把i当普通变量用。第五个坑是基底导纳的修正因子。不乘real(eta_s)/real(eta0)无吸收膜系的RT会差一点斜入射时差更多。这个因子是能量守恒的直接体现丢了它你会花费很多时间怀疑物理模型实际上只是少了一个因子。5.2 常见问题速查表现象可能原因解决办法结果全是NaN入射角过大或折射率组合导致asin参数1出现全反射检查入射角与折射率对asin结果判断实部合法性RT明显不等于1透射率缺介质导纳修正因子或膜层有吸收检查T公式中是否用了real(eta_s)/real(eta0)正入射正常、斜入射乱套逐层追角漏写了θ_j计算确认每层都用斯涅尔定律算出各自折射角反射谱中心波长偏了膜厚与波长单位不一致统一使用nm检查d_list数值是否对应四分之一光学厚度膜层越多结果越乱矩阵连乘顺序错误检查循环内 M_total M_total * M_layer光从吸收介质入射时结果不可靠复数斯涅尔定律处理不当弱吸收可近似取实部强吸收建议改用复数形式求解cosθ5.3 几个让代码更耐用的建议第一把函数封装好后统一放在一个仿真脚本里用脚本循环扫描波长、角度、膜厚这样参数扫描和参数拟合都方便。第二把结果画成双y轴图一个y轴放R一个放T不要混在一起。第三跑比较宽的波长范围时把波长间隔设小一点尤其扫描带边位置很多特殊结构对波长扫描分辨率很敏感。另外如果你需要跟商用软件比如TFCalc、Essential Macleod对比结果可以把同一膜系的折射率、厚度、入射角输入进去然后对照反射率曲线。如果发现某个波段出现整体偏移先怀疑折射率色散数据不一致如果只是个别尖峰位置有偏差再考虑相位厚度和追角这些细节。我现在跑新膜系设计时一般流程是先在这套MATLAB代码里快速扫几种候选膜系结构确认光谱趋势符合预期再上专业镀膜模拟软件做更完整的反演和公差分析。MATLAB在这里的角色是“快速验证器”虽然它没有商用软件那么多材料库和优化算法但胜在逻辑透明、自由度大想加什么参数随时都能加这对理解物理过程特别有帮助。最后再分享一个小技巧如果你只是想验证某个周期性膜系的光子禁带位置可以把周期数先设小一点比如6到8个周期先看禁带大概位置确定没有把膜的厚度算错再加大周期数跑最终结果。这样每次跑循环都很快调试效率高很多。光学仿真这事耐心和细心永远比硬件配置重要。
MATLAB光学多层膜传输矩阵法仿真：TE模反射率计算实战
做光学镀膜的人谁还没跟反射率谱线死磕过几个晚上。我记得第一次跑这套多层膜仿真代码的时候脑子里只有一个念头能不能在写膜系之前先用电脑把光谱行为摸个底别等镀膜机开炉了才发现设计反了。答案是能而且用MATLAB就能做得挺顺手。这套光学多层膜系统模拟仿真代码核心就是用传输矩阵法计算平面波在多层薄膜里的传播行为重点处理TE模的反射率、透射率顺带可以抠出电场分布。它适合三类人刚开始接触薄膜光学的学生、需要做膜系设计预研的工程师、以及想把手头的膜系设计与实测谱线互相对照的工艺人员。看懂它你就掌握了镀膜仿真最基础也最常用的一把尺子。1. 光学多层膜仿真的核心思路为什么偏要用传输矩阵法1.1 什么是光学多层膜系统它到底在模拟什么光学多层膜说白了就是把两种或多种不同折射率的材料按设计厚度一层一层堆在基底上。常见的有相机镜头的增透膜、激光器腔镜用的高反膜、冷光镜、滤光片甚至现在光学传感器里用的共振吸收结构本质上都是同一个套路。膜层厚度通常在几十纳米到几微米之间跟可见光波长一个量级所以光进去之后会发生强烈的干涉效应。宏观上表现出来就是某些波长反射增强、某些波长透射增强通过调整每层的折射率和厚度就能“雕塑”出你想要的光谱曲线。这套MATLAB仿真代码做的就是把上述物理过程数值化。它输入的是每一层材料的复折射率、厚度、入射介质折射率、基底折射率、入射角和波长范围输出的是对应的反射率R、透射率T以及吸收率A。TE模是这套代码首先处理的偏振态物理上就是电场矢量垂直于入射面的那一束光比如很多激光系统里s偏振分量就是TE模。1.2 为什么必须用数值方法而不是手推公式单层膜的反射率公式是教科书里的经典内容很多人一开始以为多层膜也能从头到尾推一个封闭表达式出来。实际上当膜层数超过三层解析表达式的复杂度就会急速膨胀而且一旦膜层有吸收折射率是复数或入射角非零手推简直是一场灾难。传输矩阵法的思路非常聪明不直接去解每层内部复杂的多次反射叠加而是把每层等效成一个2×2的矩阵这个矩阵同时包含了光在该层中传播的相位变化和界面处的边界条件。把每一层的矩阵按顺序乘起来得到一个总的2×2矩阵然后从总矩阵元素直接读出反射系数和透射系数。这个方法的物理图像清楚、代码量少、数值稳定性好是当前光学薄膜仿真领域事实上的标准方法。MATLAB擅长矩阵运算复数、矩阵连乘、循环扫描波长这些都是基础操作所以跑起来非常顺手。1.3 为什么我推荐先把TE模单独实现TE模和TM模的差别本质上是边界条件里场分量方向不同导致的。TE模的电场只有垂直于入射面的一个分量边界条件里不涉及电场法向分量的跳变所以推导和编程都相对简洁。先把TE模跑通你就能建立对传输矩阵法的直观理解相位厚度怎么算、导纳怎么代、矩阵连乘怎么排这些核心逻辑在TM模里完全复用只是导纳表达式从 n·cosθ 换成 n/cosθ 而已。另外很多实际应用场景比如某些偏振分光膜、斜入射高反膜、光栅结构里的类导模TE模本身就是关键角色。所以先实现TE模不是偷懒而是用最小成本掌握方法再往TM、椭圆偏振、色散材料、粗糙界面等方向扩展。2. 传输矩阵法与TE模的数学原理把公式变成代码前的必要铺垫2.1 多层膜问题的“降维打击”从麦克斯韦方程组到2×2矩阵光在多层膜中的传播本质上是由麦克斯韦方程组和界面连续条件控制的。如果入射光是单色平面波每个膜层内的电场可以写成前向传播波和反向传播波的叠加。在一个界面上入射、反射、透射波的切向分量必须连续于是可以把这个约束关系表达成矩阵形式把界面一侧的场与另一侧的场联系起来。而光在一层均匀介质内部从界面a走到界面b电场只会经历一个相位变化即乘上 e 的 正负iδ 次方。这个传播过程同样可以用一个对角矩阵表示。把界面矩阵和传播矩阵合在一起就得到了一个层的特征矩阵也就是传输矩阵。[ M_j \begin{bmatrix} \cos\delta_j \dfrac{i\sin\delta_j}{\eta_j} \ i\eta_j\sin\delta_j \cos\delta_j \end{bmatrix} ]其中δ_j是第j层的相位厚度η_j是该层对应偏振态的光学导纳。整个膜系的响应就是所有层的特征矩阵按光传播方向连乘[ M M_1 M_2 \cdots M_m ]这里的矩阵乘法顺序不能乱第一层先乘最后一层最后乘这直接对应光从入射介质进入膜层再到达基底的物理路径。最后反射系数和透射系数可以从总矩阵元素与入射介质、基底的导纳组合得到。这个过程可以类比电路里的级联网络每一层是一个二端口网络整体级联后等效成一个总的二端口网络反射率和透射率就是网络的散射参数。2.2 TE模与TM模的区别为什么导纳表达式不一样导纳这个概念第一次接触会觉得抽象但你可以把它理解为“介质对电磁波的等效阻抗”由折射率和入射角共同决定。对于TE模入射面内电场只有垂直分量磁场的切向分量与电场比值关系会得到[ \eta_{TE} n\cos\theta ]这里的θ是光在该层内的折射角不是入射角。对于TM模因为涉及磁场垂直入射面、电场在入射面内边界条件里的阻抗关系不同结果是[ \eta_{TM} \frac{n}{\cos\theta} ]如果正入射θ0cosθ1TE和TM退化为一模一样。这也是为什么很多人验证代码时习惯从正入射开始这时候偏振无关公式最简单。一旦斜入射就必须分开处理。这套代码里先做了TE模后面要扩展到TM只需要把计算导纳的语句改掉其他逻辑几乎不动。2.3 相位厚度、斯涅尔定律和入射角的连锁计算相位厚度是薄膜干涉的核心量它决定了光在膜层里走一个来回积累了多少相位差[ \delta_j \frac{2\pi n_j d_j \cos\theta_j}{\lambda} ]这里n_j是第j层复折射率d_j是几何厚度λ是入射光在真空中的波长。注意每一层里的折射角θ_j并不等于入射角θ_0需要通过斯涅尔定律逐层追踪[ n_0\sin\theta_0 n_1\sin\theta_1 n_2\sin\theta_2 \cdots ]也就是说虽然光从外部以某个角度入射但进入高折射率材料后折射角会变小进入低折射率材料后会变大。每一层的光路都不一样相位厚度也不一样代码里必须循环逐层算出θ_j再代进相位厚度公式。这个逐层追角的操作是很多人写代码时容易漏掉的地方漏掉之后正入射没问题但斜入射结果就全错了。3. 核心代码模块实现从零搭一个TE模传输矩阵仿真器3.1 函数骨架与输入参数设计我用MATLAB写了一个函数输入是膜层参数和入射条件输出是反射率、透射率和吸收率。函数的好处是方便反复调用比如扫描波长、扫描角度或者扫膜厚只需要在外层套一个for循环就能批量跑不同条件。这里给出核心骨架function [R, T, A] multilayer_TE(n_list, d_list, n0, ns, lambda, theta0_deg) % multilayer_TE 传输矩阵法计算TE模多层膜反射、透射、吸收 % 输入 % n_list - 各膜层复折射率行向量长度m % d_list - 各膜层几何厚度单位nm行向量长度m % n0 - 入射介质折射率可为复数 % ns - 基底折射率可为复数 % lambda - 波长单位nm可以是标量或行向量 % theta0_deg - 入射角单位度 % 输出 % R, T, A - 反射率、透射率、吸收率与lambda同维度输入里我把厚度和波长都统一用nm这样计算相位厚度时两个量纲可以消掉不会出现单位混用导致的量级错误。折射率用复数的场景稍后会说但函数从第一版就支持复数后面处理吸收介质就不用改接口。3.2 入射角度追踪与导纳初始化主函数第一步是把入射角从度转成弧度然后初始化输出数组theta0 theta0_deg * pi / 180; N length(lambda); R zeros(1, N); T zeros(1, N); A zeros(1, N);对每个波长扫描点先算出每一层里面的折射角。注意当膜层有吸收、折射率为复数时直接套实数斯涅尔定律会有点问题这里先按常规做法取实部来算角度满足绝大多数弱吸收介质场景for idx 1:N lam lambda(idx); m length(n_list); theta zeros(1, m); % 逐层用斯涅尔定律追角 theta(1) asin(real(n0) * sin(theta0) / real(n_list(1))); for j 2:m theta(j) asin(real(n_list(j-1)) * sin(theta(j-1)) / real(n_list(j))); end % 入射介质导纳 eta0 n0 * cos(theta0);有个边界情况要小心当入射角太大或者从高折射率介质入射到低折射率介质asin 参数可能略大于1这时会发生全反射。MATLAB 会返回复数数值上不至于直接崩但结果需要你多想一想物理上是否合理。后面我会在常见问题里专门展开。3.3 特征矩阵连乘与反射率计算接下来就是核心循环首先把总矩阵设为单位阵然后逐层构建特征矩阵并连乘M_total eye(2); % 2x2单位矩阵 for j 1:m delta 2 * pi * n_list(j) * d_list(j) * cos(theta(j)) / lam; eta_j n_list(j) * cos(theta(j)); M_layer [cos(delta), 1i * sin(delta) / eta_j; 1i * eta_j * sin(delta), cos(delta)]; M_total M_total * M_layer; end这段代码是整套仿真的心脏。cos(delta)和sin(delta)里的delta直接用复数n计算所以会同时计入吸收导致的振幅衰减。矩阵元素里的1i是MATLAB的虚数单位千万别把1i当成变量i用否则会跟循环索引冲突。基底导纳要单独计算把基底看作半无限介质没有反射回来的波。然后通过总矩阵元素和入射/基底导纳组合出反射系数和透射系数% 基底内的折射角 theta_s asin(real(n0) * sin(theta0) / real(ns)); eta_s ns * cos(theta_s); % 反射系数和透射系数公式 ref_num eta0 * (M_total(1,1) M_total(1,2) * eta_s) - (M_total(2,1) M_total(2,2) * eta_s); ref_den eta0 * (M_total(1,1) M_total(1,2) * eta_s) (M_total(2,1) M_total(2,2) * eta_s); r ref_num / ref_den; trans_num 2 * eta0; t trans_num / ref_den; R(idx) abs(r)^2; T(idx) real(eta_s) / real(eta0) * abs(t)^2; A(idx) 1 - R(idx) - T(idx); end透射率公式里那个 real(eta_s)/real(eta0) 是介质阻抗修正项很多初学的人查了半天发现T和R加起来不是1就是因为漏了这层介质导纳比。有了修正项无吸收膜系里RT才严格等于1。3.4 验证代码正确性的第一步单层膜解析解对照代码写完别急着跑多层先拿单层膜正入射验证。单层膜在中心波长处的反射率存在解析解比如基底为玻璃1.52膜层为氟化镁1.38设计中心波长550nm膜厚取100nm左右正入射。用上面的函数算出来的反射率曲线应该与教科书上单层介质膜的反射率公式一致。我在实际验证时经常这么操作n_list [1.38]; d_list [99.64]; n0 1.0; ns 1.52; lambda linspace(400, 800, 401); theta0_deg 0; [R, T, A] multilayer_TE(n_list, d_list, n0, ns, lambda, theta0_deg); plot(lambda, R);这个膜厚对应的是550nm的四分之一光学厚度550/(4×1.38) ≈ 99.64nm。在550nm附近R出现一个极小值这就是增透谷。拿这个结果与理论解析式比较一致就说明矩阵构建和边界条件处理是对的。这是最快、最有效的“自检”手段胜过任何复杂膜系的直觉判断。3.5 扩展把电场分布也算出来反射率和透射率只是冰山一角很多情况下你需要知道光在膜系内部的电场分布。比如激光损伤研究里损伤经常出现在电场峰值位置设计吸收增强结构时吸收层要放到电场最大处。这个需求在代码里扩展起来不算难核心思路是得到了总矩阵M和反射系数r、透射系数t之后逐层向前递推可以求出每一层入射界面处的复振幅。function E_field calc_Efield(n_list, d_list, n0, ns, lambda, theta0_deg) % 先调用主函数得到反射系数和透射系数再逐层递推电场 % 这里给出递推的核心表达式省略了完整封装 E_plus 1; % 入射电场振幅归一化 E_minus r; % 反射电场振幅 % 在每一层入口处前向波 E_plus * exp(-i*delta_partial) % 后向波 E_minus * exp(i*delta_partial) % 总电场 前向波 后向波这块代码如果展开本质上是对每一层的传播矩阵做逐步“剥洋葱”从空气一侧往基底方向推进把每一层入口和出口的电场算出来。有了这个扩展你的仿真能力就从“只看光谱曲线”跳到“看懂场分布”这对分析干涉效应、吸收层摆放位置、镀膜应力带来的光程变化都有价值。4. 仿真结果怎么看那些藏在曲线里的信息4.1 典型高反膜堆的反射率谱线解读用上面这套代码跑一个最常见的膜系中心波长1064nm(HL)^15 高反膜堆H材料用二氧化钛n2.35L材料用二氧化硅n1.46基底是玻璃正入射。每个周期中H层厚度为1064/(4×2.35)L层厚度为1064/(4×1.46)。你会发现中心波长附近出现一个高反射带带宽大约在几百纳米量级带内反射率非常接近1。这个高反射带的物理来源是布拉格反射每个界面的反射波在中心波长附近同相叠加类似于X射线在晶体上的衍射。通过调整折射率比和周期数你能控制带宽和反射率。折射率比越大反射带越宽周期数越多带内反射率越平、越接近100%。用这个代码去扫描不同周期数的膜系你会直观看到“增加膜层数”不是无限有效膜层数超过一定数目后中心反射率变化非常小但膜层厚度和应力都在涨工程上必须做取舍。4.2 光子禁带与布拉格条件一维光子晶体在结构上跟高反膜堆完全相同区别在于你是否把它放到周期结构语境里理解。当一个周期内两种材料的相位厚度之和满足[ 2(n_H d_H \cos\theta_H n_L d_L \cos\theta_L) m\lambda ]就出现布拉格反射对应频率范围内的光无法在周期结构中传播这就是所谓的光子禁带。斜入射时θ_H和θ_L会随入射角变化禁带位置也会移动这就是很多角度选择滤光片的设计基础。你的代码里直接把入射角从0度扫到60度就能看到反射带逐渐蓝移的整个过程非常直观。4.3 电场分布揭示了什么很多人只盯反射率忽略了电场分布但在一些实际开发里电场分布比反射率更致命。曾经我做介质高反膜激光损伤实验发现膜层总是从某个特定界面开始坏一开始怀疑材料纯度后来用仿真把1064nm波长下的电场分布算出来才发现那个界面临近电场波腹电场强度是入射场的两倍以上。从那以后设计高功率激光膜系时我都会顺手把电场分布打印出来看峰值有没有落在薄膜界面附近。在λ/4高反膜堆里驻波场的波节通常落在高折射率层中心波腹落在低折射率层中心或界面附近。如果你把吸收较强的材料放到波腹位置吸收损耗会被显著放大反过来把吸收小的材料放在波腹区吸收型材料放在波节区就能实现“低吸收、高损伤阈值”的设计。这就是电场分布仿真在镀膜工艺里的价值——它让你看到光在膜系里到底“待在哪里”。5. 实操总结避坑指南与调试建议5.1 我踩过的五个坑从NaN到静默错误第一个坑是矩阵左乘右乘搞反。总矩阵必须按照光传播方向依次乘第一层在最左边。这个顺序错了单层膜结果也可能“碰巧”是对的但膜层一多就彻底乱套。第二次写代码时我就因为循环里 M_total M_total * M_layer 写反耗时大半天查不出来。第二个坑是单位不统一。厚度用微米、波长用纳米算出来的相位厚度差了1000倍反射谱看起来完全不是你以为的那个膜系。现在我的代码里给d_list和lambda都加了注释一律要求nm函数入口处再做一次检查如果发现数量级不对就提示。第三个坑是正入射时cosθ都是1很多分支逻辑不会被触发结果一切正常一旦切到斜入射各种各样的索引错误、NaN就冒出来了。这提示我正入射验证只是一个必要不充分条件斜入射至少选一个30度、60度做持续测试。第四个坑是虚数单位i和循环索引打架。MATLAB内置的虚数单位是1i如果我在某个地方顺手写了i 1这种赋值后面所有1i都会出问题。现在我一直用1i绝不把i当普通变量用。第五个坑是基底导纳的修正因子。不乘real(eta_s)/real(eta0)无吸收膜系的RT会差一点斜入射时差更多。这个因子是能量守恒的直接体现丢了它你会花费很多时间怀疑物理模型实际上只是少了一个因子。5.2 常见问题速查表现象可能原因解决办法结果全是NaN入射角过大或折射率组合导致asin参数1出现全反射检查入射角与折射率对asin结果判断实部合法性RT明显不等于1透射率缺介质导纳修正因子或膜层有吸收检查T公式中是否用了real(eta_s)/real(eta0)正入射正常、斜入射乱套逐层追角漏写了θ_j计算确认每层都用斯涅尔定律算出各自折射角反射谱中心波长偏了膜厚与波长单位不一致统一使用nm检查d_list数值是否对应四分之一光学厚度膜层越多结果越乱矩阵连乘顺序错误检查循环内 M_total M_total * M_layer光从吸收介质入射时结果不可靠复数斯涅尔定律处理不当弱吸收可近似取实部强吸收建议改用复数形式求解cosθ5.3 几个让代码更耐用的建议第一把函数封装好后统一放在一个仿真脚本里用脚本循环扫描波长、角度、膜厚这样参数扫描和参数拟合都方便。第二把结果画成双y轴图一个y轴放R一个放T不要混在一起。第三跑比较宽的波长范围时把波长间隔设小一点尤其扫描带边位置很多特殊结构对波长扫描分辨率很敏感。另外如果你需要跟商用软件比如TFCalc、Essential Macleod对比结果可以把同一膜系的折射率、厚度、入射角输入进去然后对照反射率曲线。如果发现某个波段出现整体偏移先怀疑折射率色散数据不一致如果只是个别尖峰位置有偏差再考虑相位厚度和追角这些细节。我现在跑新膜系设计时一般流程是先在这套MATLAB代码里快速扫几种候选膜系结构确认光谱趋势符合预期再上专业镀膜模拟软件做更完整的反演和公差分析。MATLAB在这里的角色是“快速验证器”虽然它没有商用软件那么多材料库和优化算法但胜在逻辑透明、自由度大想加什么参数随时都能加这对理解物理过程特别有帮助。最后再分享一个小技巧如果你只是想验证某个周期性膜系的光子禁带位置可以把周期数先设小一点比如6到8个周期先看禁带大概位置确定没有把膜的厚度算错再加大周期数跑最终结果。这样每次跑循环都很快调试效率高很多。光学仿真这事耐心和细心永远比硬件配置重要。
