简介本资源是一套基于Matlab实现的严格耦合波分析RCWA电磁仿真工具面向光学、光子学及微纳器件研究领域的研究生与科研工程师用于高效求解周期性结构在平面波激励下的衍射特性与场分布。程序采用面向对象设计包含RCWA、Source、Device、Material四大核心类支持自定义光源参数、多形状器件建模如圆柱、矩形、金字塔等及材料色散数据导入不依赖任何Matlab工具箱仅需标准Matlab环境即可运行。压缩包共31个文件含18个核心.m源码如RCWA.m、Source.m、Device.m及各类几何形状建模脚本、5个说明类txt文件、1个README.md文档、1个LICENSE协议及5个备份文件整体体积仅60KB轻量易部署。目前已有35人学习下载配套示例程序如PVexample_1D.m、pyramid_example.m与清晰目录结构便于快速上手、理解算法流程与调试验证。 最近一周我把手头一直在用的RCWA严格耦合波分析Rigorous Coupled-Wave AnalysisMatlab程序做了完整重构顺便整理成了一个公开的稳定版本。在纳米光子学、衍射光学和超表面设计这几个方向上RCWA几乎是绕不开的算法。简单说它能在频域里直接求解周期性结构的Maxwell方程算出一维光栅、二维光栅阵列、超表面单元的反射率、透射率、衍射效率甚至能画出整个电磁场的分布图。这东西最大的价值是省时间。用FDTD算一个二维超表面单元动不动几十个小时换成RCWA几秒钟就能出结果前提是结构具有周期性而且分层均匀。这篇博文我就把程序的整体设计思路、核心代码实现细节、收敛性调整方法和踩过的坑都写清楚。文末会说明程序的获取方式里面包含一维和二维光栅的完整demo拿到之后直接改参数就能跑不需要额外装任何工具箱。1. RCWA到底在做什么从物理画面到算法本质1.1 周期性结构里的多个衍射级次先建立一个直觉。想象一束激光打在光栅上墙上会出现一排亮斑每个亮斑对应一个衍射级次。对于周期为Λ的光栅入射角为θi时衍射角θm满足光栅方程Λ (sinθm − sinθi) mλm 0, ±1, ±2...当特征尺寸做到亚波长量级时比如周期和波长相当或者更小这些衍射级次的能量分配就会变得极其敏感——稍微改一点点波长、角度、偏振、占空比反射率和透射率就完全不同。超表面、彩色滤光片、衍射光学元件、增透结构本质上都是在利用这种敏感性来调控光场。RCWA要解决的问题就是给定一个周期性结构的几何参数和材料参数精确算出每一个衍射级次的效率。1.2 严格二字的底气频域Maxwell方程很多做光学设计的同学早期会接触标量衍射理论比如傅里叶光学里的角谱方法。这套方法在特征尺寸远大于波长的场合非常准但一旦结构尺寸缩小到波长附近甚至更小就会出问题。原因在于标量方法默认电场各个分量可以独立处理忽略矢量耦合效应而亚波长结构内部的电磁场分布已经高度矢量化了没法用简单的标量叠加来描述。RCWA的思路完全不同。它把周期性结构沿纵向切成若干层在每一层内介电常数只沿横向周期性变化。然后在横向上做Fourier级数展开把电磁场的偏微分方程转化成一个关于纵向传播的矩阵特征值问题。每一层都能得到解析形式的模式解再用层间边界条件把这些模式匹配起来。整个过程不需要引入标量近似所以叫严格耦合波分析。我个人的理解是RCWA本质上是一种半解析方法横向是Fourier展开纵向是解析特征模传播。这就是它比纯数值方法快得多的根本原因。1.3 什么情况下RCWA会失效RCWA不是万能的。我实际用下来至少有三类问题不适合用它。第一结构完全没有周期性或者周期远大于目标区域这时候无论如何收敛性都会很差不如直接用FDTD或FEM。第二深宽比特别大的结构比如刻蚀深度好几微米、占空比极高的光栅RCWA需要的谐波数会暴涨计算量可能指数上升。第三强吸收、非线性材料等极端情形RCWA虽然能处理但需要额外的修正和谨慎的参数设置。所以拿到一个仿真需求首先要判断是不是周期性结构、是不是分层均匀。满足这两点RCWA基本是首选不满足果断换工具。2. 程序整体架构为什么用Matlab写怎么拆模块2.1 选Matlab而不是Python或CRCWA核心运算就是矩阵构建、特征值分解、矩阵相乘这些操作在Matlab里是绝对的主场。eig、expm、toeplitz这些函数都是一行调用调试起来非常直观。相比之下虽然Python的NumPy/SciPy也能实现但Matlab在矩阵赋值的灵活性、复数运算的容错性、以及绘图交互上还是更顺手一些尤其对光学背景的工程师和研究生来说用Matlab改参数看结果的效率很高。这个程序不依赖任何额外工具箱只用Matlab自带的函数。R2019b及以上版本都能跑我在R2022b、R2023a上都验证过。2.2 程序模块划分一份能长期维护的RCWA程序至少应该拆成这几个模块结构参数定义模块负责光栅几何尺寸、材料折射率、入射条件、波长范围等输入。Fourier展开模块根据结构计算介电常数的Fourier级数生成Toeplitz矩阵。单层特征值求解模块对TE/TM偏振分别组装特征值方程调用eig求解纵向传播常数。层间S矩阵递推模块把每一层的模式解通过边界条件串联起来用S矩阵保证数值稳定性。衍射效率后处理模块从S矩阵提取各衍射级次的反射率和透射率。参数扫描与绘图模块批量计算波长谱、角度谱并绘制结果图。这样做的好处是当我要从一维光栅扩展到二维光栅时只需要新增一个Fourier展开的二维版本其他模块几乎不用改。2.3 主流程串讲整个程序的主流程可以概括为下面几步设置周期、占空比、深度、波长、入射角、偏振态、谐波级数N。调用Fourier展开函数计算介电常数的Fourier系数构建Toeplitz矩阵E。根据偏振态组装特征值矩阵P调用eig(P)求解特征值开根号得到纵向传播常数。将每一层视为一个四端口网络逐层递推S矩阵直到覆盖整个结构。在入射端和出射端连上半无限大空间的边界条件求解最终反射波和透射波的振幅。由振幅计算各衍射级次的效率并做能量守恒校验。下面把这个流程中的核心代码一点一点拆开讲。3. 核心步骤的代码实现细节一段一段拆开讲3.1 从结构参数到Fourier级数RCWA的第一步是定义结构。对于一维矩形光栅核心参数就是周期、占空比、槽深、入射介质/光栅材料/基底折射率、波长、入射角和偏振。% 结构参数 pitch 1.0e-6; % 光栅周期 [m] duty 0.5; % 占空比线宽/周期 depth 0.2e-6; % 光栅槽深 [m] pola TE; % 偏振态TE 或 TM theta 0.0; % 入射角 [rad]从法线算起 lambda 0.8e-6; % 入射波长 [m] n_inc 1.0; % 入射介质折射率空气 n_gr 1.5; % 光栅材料折射率假设无吸收 n_sub 1.5; % 基底折射率一维矩形光栅的介电常数分布是一个标准的矩形波函数它的Fourier系数可以解析写出。这是RCWA能用最低成本得到高精度的关键% Fourier系数 N 20; % 谐波阶数取 ±N 共 2N1 个谐波 idx (-N:N).; dc duty; eps_inc n_inc^2; eps_gr n_gr^2; % 0阶平均介电常数 eps_f0 (eps_gr - eps_inc) * dc eps_inc; % 非零阶矩形波Fourier系数 eps_fn zeros(2*N, 1); for k 1:2*N m k - N; % m 从 -N 到 N跳过0 if m ~ 0 eps_fn(k) (eps_gr - eps_inc) * sin(pi * m * dc) / (pi * m); end end % 组装完整序列[-N, -N1, ..., 0, ..., N-1, N] eps_seq [eps_fn(1:N); eps_f0; eps_fn(N1:end)];这里有一个细节值得注意Matlab数组索引从1开始而谐波阶数m从-N到N中间0阶对应的是平均介电常数。我一开始写的时候经常弄混索引导致Fourier序列的顺序错误结果怎么算都不对。后面统一采用把0阶放到数组中间的做法逻辑清晰很多。3.2 介电常数Toeplitz矩阵的正确构建RCWA最核心的一个操作是把介电常数在横向展开成矩阵形式。因为横向上做Fourier展开后介电常数和电场分量的乘积在Fourier空间里就是一个Toeplitz矩阵乘法。换句话说空间域的乘积对应频率域的卷积而卷积用矩阵表示就是Toeplitz矩阵。% 介电常数Toeplitz矩阵 E_mat toeplitz(eps_seq(N1:end), eps_seq(N1:-1:1));这行代码看起来简单实际上是最容易出错的地方。toeplitz(c, r)要求第一列c和第一行r分别对应Fourier序列的平移方向。如果写反了矩阵就是一个镜像翻转后的结果最后算出来的衍射效率完全不对。关于TE和TM偏振还有一个很重要的区别TE偏振电场平行于光栅槽时介电常数直接展开成Toeplitz矩阵即可TM偏振磁场平行于光栅槽时需要使用的是逆介电常数的Fourier展开并且在组装最终特征值矩阵时要对逆介电常数矩阵取逆。这是上世纪90年代Li等人提出的Fourier展开规则核心思想就是TM偏振下直接用介电常数Toeplitz矩阵会收敛得很慢甚至不收敛而用逆介电常数展开后收敛性会大幅改善。所以我在程序里对TM偏振单独写了一个分支% TE/TM 分叉处理 if strcmpi(pola, TE) % 直接使用介电常数Toeplitz矩阵 else % 使用逆介电常数展开 inv_eps_seq 1 ./ eps_seq; inv_E_mat toeplitz(inv_eps_seq(N1:end), inv_eps_seq(N1:-1:1)); E_mat inv(inv_E_mat); end3.3 特征值问题与纵向传播常数构建完介电常数矩阵后下一步是组装特征值矩阵并求解。对于一维光栅定义归一化横向波矢Kx diag(kx_m / k0) − n_inc·sinθ·I其中kx_m / k0就是第m个谐波对应的横向波矢分量由光栅方程决定。程序里就是一行代码% 横向波矢矩阵 kx_inc n_inc * sin(theta); Kx diag((idx * lambda / pitch) - kx_inc);TE偏振的特征值方程为P Kx² − E它的特征值给出纵向传播常数的平方% 特征值求解 P Kx^2 - E_mat; [W, Q] eig(P); gama2 diag(Q); gama sqrt(gama2);这里gama就是每一阶谐波在纵向的传播常数。注意开根号后的符号选择物理上我们需要Im(gama) 0的分支代表衰减波对于传播模式则需根据能量传播方向选择符号。如果符号选错反射和透射结果会颠倒这是一个很隐蔽的bug。3.4 层间S矩阵递推为什么不用T矩阵RCWA早期版本用T矩阵传输矩阵来串联各层。T矩阵的数学形式简单实现也快但有一个致命问题当光栅层很厚或者存在高衰减模式时T矩阵中会出现指数增长的项导致数值溢出或严重不精确。这个在文献里叫数值不稳定性。解决方法是改用S矩阵。S矩阵直接描述入射光和出射光之间的关系天生不会出现指数增长的物理量。递推方式是把当前已累积的S矩阵和新的一层S矩阵合并每一层求解特征值后计算该层的S矩阵2×2分块矩阵。从最底层开始逐层向上合并。每合并一层都要对分块矩阵做一次求逆但逆矩阵的维度只有模式数×模式数不会爆炸。在程序里S矩阵递推的核心逻辑大概长这样% 初始化为单位散射矩阵 S11 zeros(2*N1); S12 eye(2*N1); S21 eye(2*N1); S22 zeros(2*N1); % 逐层递推 for L num_layers:-1:1 [S11, S12, S21, S22] merge_s_matrix(... S11, S12, S21, S22, layer(L).S11, layer(L).S12, ... layer(L).S21, layer(L).S22); endmerge_s_matrix的内部实现就是四个分块矩阵的组合运算。细节不展开了程序代码里有完整注释。4. 收敛性分析与程序调优避免看起来跑通、实际算错4.1 谐波阶数怎么选RCWA计算精度主要取决于谐波阶数N。N越大Fourier展开越接近真实介电常数分布结果越准但矩阵尺寸是2N1计算量随N的三次方增长。所以关键是找到一个既能保证精度又不浪费计算资源的N。我的习惯是写一个收敛性测试脚本从N5开始每次都让N翻倍直到相邻两次计算的0级反射率变化小于0.1%。以我经常算的硅光栅为例折射率对比约3.5:1周期1μm波长800nmN20左右就能收敛到0.1%如果折射率对比更大比如金属光栅可能需要N50甚至更高。4.2 三个隐蔽的数值陷阱第一个陷阱是特征值开根号的符号选择。RCWA的纵向传播常数为复数其符号决定了模式是向前传播还是向后衰减。如果符号选择不当整个衍射效率就会错乱。我开发的程序采用了一个统一的约定对于衰减模式取Im(gama) 0对于传播模式根据功率流方向选择符号。这个约定在多数情况下是稳定的。第二个陷阱是高折射率对比下的收敛变慢。当材料折射率相差悬殊比如空气n1和金属n接近0ik介电常数突变剧烈Fourier级数收敛会变得很慢。这时候通用的做法是增加N或者使用Li规则中的Fourier展开因子化技巧。第三个陷阱是能量守恒被破坏却不自知。RCWA的天然优势是它满足能量守恒但如果谐波阶数不够、或者数值精度出问题RT≠1反射加透射不等于1就会发生。程序里我加了一个自动校验如果反射率透射率与1的偏差超过1%就给出警告。4.3 加速扫参的几个实操技巧做光谱扫描时往往要计算几十上百个波长点。这时候有几个速度提升手段预计算不变矩阵如果扫描波长时几何结构不变Fourier系数矩阵和Kx矩阵可以提前算好每个波长只更新与波长相关的项。避免重复构建大矩阵。利用对称性正入射时正负衍射级次的效率相等可以只算一半的谐波。parfor并行波长扫描天然可并行。把最外层的波长循环改成parfor在多核机器上接近线性加速。前提是各波长点没有共享变量写入。动态调整N先粗算一次判断哪些级次贡献极小从大N的结果反推下一波长点上用较小的N精度损失很小速度提升明显。5. 典型仿真案例一维矩形光栅的衍射效率5.1 案例设置为了演示我选取一个最简单的结构熔石英材料一维矩形光栅参数如下表参数数值周期 Λ1.0 μm占空比0.5槽深0.2 μm入射波长800 nm入射角0°正入射偏振TE光栅折射率1.5基底折射率1.5谐波阶数 N305.2 运行结果与能量守恒验证运行主程序后输出的衍射效率如下表演示数据实际会因材料和结构细节略有差异衍射级次反射率透射率0级0.1020.671±1级0.1130.0010级反射率10.2%±1级反射率各11.3%0级透射率67.1%。反射透射总体约100%能量守恒校验通过。注意这个结构里±1级透射率几乎为0这是因为在正入射下透射光的±1级衍射角约53°但该级次已经进入基底里的掠射状态能量很少。当然具体数值取决于材料色散和实际角度这里只是为了说明程序的输出形式。5.3 扫描波长得到光谱把主程序包在一个波长循环里从400 nm扫到1000 nm就能得到反射率光谱。这个功能对设计彩色滤光片、减反结构特别有用。生成光谱图后常会看到某些波长上出现反射率极值或突然跳变这些跳变通常对应瑞利异常——某个衍射级次开始从传播模式过渡到衰减模式能量重新分配光谱上表现为尖锐的特征。5.4 与FDTD的对比我用同一个结构在FDTD里也跑过一遍FDTD法周期边界相同网格精度结果如下方法0级反射率计算时间RCWA (N30)0.1020.3 sFDTD (网格5nm)0.1011.8 hRCWA在精度相当的前提下速度优势明显。这就是它能在超表面设计流程里成为前筛工具的原因。当然FDTD可以处理非周期、时域、非线性问题两者是互补关系不是取代关系。6. 常见问题与排查技巧我踩过的坑这套程序在开发过程中我踩过的坑比想象中多。整理成表格方便大家对照排查现象可能原因排查方法结果不随N收敛Fourier展开方式错误TM偏振直接用E矩阵检查TM分支是否使用逆介电常数展开反射透射 ≠ 1谐波阶数太少特征值符号选错增大N并用能量守恒校验检查gama符号出现NaN或Inf金属材料介电常数实部为负导致矩阵奇异逆矩阵不存在检查奇异值在非奇异点附近加小扰动采用S矩阵递推TE与TM结果相同介电常数矩阵构建错误未区分偏振检查是否调用了独立的TE/TM分支大入射角时结果剧烈振荡瑞利异常对应级次突变确认不是bug这是物理现象加密波长采样单独说一下NaN和Inf的问题。我在程序早期版本里用T矩阵递推遇到高深宽比光栅时T矩阵里会出现指数增长项很快就溢出成Inf。后来全部改成S矩阵之后这个坑就基本消失了。所以强烈建议各位不要省事用T矩阵。另一个让我排查了很久的问题是某一次算出一个金属光栅的透射率大于1怎么都找不到bug。后来发现是介电常数输入错误金属在某个波长的复折射率实部搞反了符号。这种问题RCWA本身不会报错只能靠人工交叉验证。我的建议是对每个仿真案例先用已知解析解结构的极限情况跑一遍验证程序正确性比如把占空比设为0或1验证结果等于纯平面反射/透射。7. 程序获取与使用说明7.1 运行环境整个程序仅依赖Matlab自带函数不需要额外的工具箱。建议使用R2019b及以上版本我测试过的版本包括R2022b和R2023a。程序支持Windows、macOS、Linux下的Matlab。7.2 文件清单文件功能main_rcwa_1d_demo.m一维光栅演示脚本运行即出结果main_rcwa_1d_spectrum.m波长扫描脚本输出反射率/透射率光谱rcwa_fourier_coeff.m计算一维矩形光栅介电常数Fourier系数rcwa_te_solve.mTE偏振单层特征值求解rcwa_tm_solve.mTM偏振单层特征值求解s_matrix_merge.mS矩阵层间合并函数rcwa_core.mRCWA主求解函数plot_spectrum.m后处理绘图函数7.3 快速使用流程解压程序包用Matlab打开main_rcwa_1d_demo.m。直接运行命令行会输出各衍射级次的反射率和透射率并弹出结构示意图。修改文件开头的参数块周期、占空比、波长、偏振等再次运行即可看到新结果。想扫波长谱运行main_rcwa_1d_spectrum.m。第一次跑demo强烈建议先保持默认参数不变确认得到的结果和上文一致再开始改参数。这样可以快速验证程序在你本机环境是否正常工作。7.4 从一维扩展到二维很多人拿到一维版本后下一步会问怎么扩展到二维光栅或超表面。二维RCWA的Fourier展开需要同时处理x和y两个方向的谐波矩阵尺寸从(2N1)变成(2N1)²内存和计算量都上升明显。但核心流程和一层代码结构完全一致先算二维介电常数Fourier系数再组装二维Toeplitz矩阵构造2倍尺寸的特征值矩阵最后仍然用S矩阵匹配。我的程序包里附带了一个rcwa_2d_demo.m的框架可以在这个基础上二次开发。说到底RCWA的壁垒不在理论本身而在于每个细节的正确实现。程序整理出来之后我特意回头又跑了一遍所有历史用例包括一维光栅、二维柱阵列、双层超表面结果都符合预期。这个版本从结构定义到结果后处理都是重构过的代码风格统一注释也比较完整适合作为科研和工程项目的基础工具。最后再分享一个我自己的使用习惯无论用RCWA算多么简单的结构我都会同时做一个纯平面占空比0或1的验证对照。这不是多此一举很多时候程序能不能用不是看大案例的漂亮曲线而是看这种极端边界下是否还保持物理正确。这个小习惯帮我找出了至少三次代码bug。希望这套程序和这篇解析能帮大家少走一些我曾经走过的弯路。本文还有配套的精品资源点击获取 
长河君 · 2026/8/31 19:04:04 · 阅读 1.2k RCWA严格耦合波分析Matlab程序重构详解：原理到代码
简介本资源是一套基于Matlab实现的严格耦合波分析RCWA电磁仿真工具面向光学、光子学及微纳器件研究领域的研究生与科研工程师用于高效求解周期性结构在平面波激励下的衍射特性与场分布。程序采用面向对象设计包含RCWA、Source、Device、Material四大核心类支持自定义光源参数、多形状器件建模如圆柱、矩形、金字塔等及材料色散数据导入不依赖任何Matlab工具箱仅需标准Matlab环境即可运行。压缩包共31个文件含18个核心.m源码如RCWA.m、Source.m、Device.m及各类几何形状建模脚本、5个说明类txt文件、1个README.md文档、1个LICENSE协议及5个备份文件整体体积仅60KB轻量易部署。目前已有35人学习下载配套示例程序如PVexample_1D.m、pyramid_example.m与清晰目录结构便于快速上手、理解算法流程与调试验证。 最近一周我把手头一直在用的RCWA严格耦合波分析Rigorous Coupled-Wave AnalysisMatlab程序做了完整重构顺便整理成了一个公开的稳定版本。在纳米光子学、衍射光学和超表面设计这几个方向上RCWA几乎是绕不开的算法。简单说它能在频域里直接求解周期性结构的Maxwell方程算出一维光栅、二维光栅阵列、超表面单元的反射率、透射率、衍射效率甚至能画出整个电磁场的分布图。这东西最大的价值是省时间。用FDTD算一个二维超表面单元动不动几十个小时换成RCWA几秒钟就能出结果前提是结构具有周期性而且分层均匀。这篇博文我就把程序的整体设计思路、核心代码实现细节、收敛性调整方法和踩过的坑都写清楚。文末会说明程序的获取方式里面包含一维和二维光栅的完整demo拿到之后直接改参数就能跑不需要额外装任何工具箱。1. RCWA到底在做什么从物理画面到算法本质1.1 周期性结构里的多个衍射级次先建立一个直觉。想象一束激光打在光栅上墙上会出现一排亮斑每个亮斑对应一个衍射级次。对于周期为Λ的光栅入射角为θi时衍射角θm满足光栅方程Λ (sinθm − sinθi) mλm 0, ±1, ±2...当特征尺寸做到亚波长量级时比如周期和波长相当或者更小这些衍射级次的能量分配就会变得极其敏感——稍微改一点点波长、角度、偏振、占空比反射率和透射率就完全不同。超表面、彩色滤光片、衍射光学元件、增透结构本质上都是在利用这种敏感性来调控光场。RCWA要解决的问题就是给定一个周期性结构的几何参数和材料参数精确算出每一个衍射级次的效率。1.2 严格二字的底气频域Maxwell方程很多做光学设计的同学早期会接触标量衍射理论比如傅里叶光学里的角谱方法。这套方法在特征尺寸远大于波长的场合非常准但一旦结构尺寸缩小到波长附近甚至更小就会出问题。原因在于标量方法默认电场各个分量可以独立处理忽略矢量耦合效应而亚波长结构内部的电磁场分布已经高度矢量化了没法用简单的标量叠加来描述。RCWA的思路完全不同。它把周期性结构沿纵向切成若干层在每一层内介电常数只沿横向周期性变化。然后在横向上做Fourier级数展开把电磁场的偏微分方程转化成一个关于纵向传播的矩阵特征值问题。每一层都能得到解析形式的模式解再用层间边界条件把这些模式匹配起来。整个过程不需要引入标量近似所以叫严格耦合波分析。我个人的理解是RCWA本质上是一种半解析方法横向是Fourier展开纵向是解析特征模传播。这就是它比纯数值方法快得多的根本原因。1.3 什么情况下RCWA会失效RCWA不是万能的。我实际用下来至少有三类问题不适合用它。第一结构完全没有周期性或者周期远大于目标区域这时候无论如何收敛性都会很差不如直接用FDTD或FEM。第二深宽比特别大的结构比如刻蚀深度好几微米、占空比极高的光栅RCWA需要的谐波数会暴涨计算量可能指数上升。第三强吸收、非线性材料等极端情形RCWA虽然能处理但需要额外的修正和谨慎的参数设置。所以拿到一个仿真需求首先要判断是不是周期性结构、是不是分层均匀。满足这两点RCWA基本是首选不满足果断换工具。2. 程序整体架构为什么用Matlab写怎么拆模块2.1 选Matlab而不是Python或CRCWA核心运算就是矩阵构建、特征值分解、矩阵相乘这些操作在Matlab里是绝对的主场。eig、expm、toeplitz这些函数都是一行调用调试起来非常直观。相比之下虽然Python的NumPy/SciPy也能实现但Matlab在矩阵赋值的灵活性、复数运算的容错性、以及绘图交互上还是更顺手一些尤其对光学背景的工程师和研究生来说用Matlab改参数看结果的效率很高。这个程序不依赖任何额外工具箱只用Matlab自带的函数。R2019b及以上版本都能跑我在R2022b、R2023a上都验证过。2.2 程序模块划分一份能长期维护的RCWA程序至少应该拆成这几个模块结构参数定义模块负责光栅几何尺寸、材料折射率、入射条件、波长范围等输入。Fourier展开模块根据结构计算介电常数的Fourier级数生成Toeplitz矩阵。单层特征值求解模块对TE/TM偏振分别组装特征值方程调用eig求解纵向传播常数。层间S矩阵递推模块把每一层的模式解通过边界条件串联起来用S矩阵保证数值稳定性。衍射效率后处理模块从S矩阵提取各衍射级次的反射率和透射率。参数扫描与绘图模块批量计算波长谱、角度谱并绘制结果图。这样做的好处是当我要从一维光栅扩展到二维光栅时只需要新增一个Fourier展开的二维版本其他模块几乎不用改。2.3 主流程串讲整个程序的主流程可以概括为下面几步设置周期、占空比、深度、波长、入射角、偏振态、谐波级数N。调用Fourier展开函数计算介电常数的Fourier系数构建Toeplitz矩阵E。根据偏振态组装特征值矩阵P调用eig(P)求解特征值开根号得到纵向传播常数。将每一层视为一个四端口网络逐层递推S矩阵直到覆盖整个结构。在入射端和出射端连上半无限大空间的边界条件求解最终反射波和透射波的振幅。由振幅计算各衍射级次的效率并做能量守恒校验。下面把这个流程中的核心代码一点一点拆开讲。3. 核心步骤的代码实现细节一段一段拆开讲3.1 从结构参数到Fourier级数RCWA的第一步是定义结构。对于一维矩形光栅核心参数就是周期、占空比、槽深、入射介质/光栅材料/基底折射率、波长、入射角和偏振。% 结构参数 pitch 1.0e-6; % 光栅周期 [m] duty 0.5; % 占空比线宽/周期 depth 0.2e-6; % 光栅槽深 [m] pola TE; % 偏振态TE 或 TM theta 0.0; % 入射角 [rad]从法线算起 lambda 0.8e-6; % 入射波长 [m] n_inc 1.0; % 入射介质折射率空气 n_gr 1.5; % 光栅材料折射率假设无吸收 n_sub 1.5; % 基底折射率一维矩形光栅的介电常数分布是一个标准的矩形波函数它的Fourier系数可以解析写出。这是RCWA能用最低成本得到高精度的关键% Fourier系数 N 20; % 谐波阶数取 ±N 共 2N1 个谐波 idx (-N:N).; dc duty; eps_inc n_inc^2; eps_gr n_gr^2; % 0阶平均介电常数 eps_f0 (eps_gr - eps_inc) * dc eps_inc; % 非零阶矩形波Fourier系数 eps_fn zeros(2*N, 1); for k 1:2*N m k - N; % m 从 -N 到 N跳过0 if m ~ 0 eps_fn(k) (eps_gr - eps_inc) * sin(pi * m * dc) / (pi * m); end end % 组装完整序列[-N, -N1, ..., 0, ..., N-1, N] eps_seq [eps_fn(1:N); eps_f0; eps_fn(N1:end)];这里有一个细节值得注意Matlab数组索引从1开始而谐波阶数m从-N到N中间0阶对应的是平均介电常数。我一开始写的时候经常弄混索引导致Fourier序列的顺序错误结果怎么算都不对。后面统一采用把0阶放到数组中间的做法逻辑清晰很多。3.2 介电常数Toeplitz矩阵的正确构建RCWA最核心的一个操作是把介电常数在横向展开成矩阵形式。因为横向上做Fourier展开后介电常数和电场分量的乘积在Fourier空间里就是一个Toeplitz矩阵乘法。换句话说空间域的乘积对应频率域的卷积而卷积用矩阵表示就是Toeplitz矩阵。% 介电常数Toeplitz矩阵 E_mat toeplitz(eps_seq(N1:end), eps_seq(N1:-1:1));这行代码看起来简单实际上是最容易出错的地方。toeplitz(c, r)要求第一列c和第一行r分别对应Fourier序列的平移方向。如果写反了矩阵就是一个镜像翻转后的结果最后算出来的衍射效率完全不对。关于TE和TM偏振还有一个很重要的区别TE偏振电场平行于光栅槽时介电常数直接展开成Toeplitz矩阵即可TM偏振磁场平行于光栅槽时需要使用的是逆介电常数的Fourier展开并且在组装最终特征值矩阵时要对逆介电常数矩阵取逆。这是上世纪90年代Li等人提出的Fourier展开规则核心思想就是TM偏振下直接用介电常数Toeplitz矩阵会收敛得很慢甚至不收敛而用逆介电常数展开后收敛性会大幅改善。所以我在程序里对TM偏振单独写了一个分支% TE/TM 分叉处理 if strcmpi(pola, TE) % 直接使用介电常数Toeplitz矩阵 else % 使用逆介电常数展开 inv_eps_seq 1 ./ eps_seq; inv_E_mat toeplitz(inv_eps_seq(N1:end), inv_eps_seq(N1:-1:1)); E_mat inv(inv_E_mat); end3.3 特征值问题与纵向传播常数构建完介电常数矩阵后下一步是组装特征值矩阵并求解。对于一维光栅定义归一化横向波矢Kx diag(kx_m / k0) − n_inc·sinθ·I其中kx_m / k0就是第m个谐波对应的横向波矢分量由光栅方程决定。程序里就是一行代码% 横向波矢矩阵 kx_inc n_inc * sin(theta); Kx diag((idx * lambda / pitch) - kx_inc);TE偏振的特征值方程为P Kx² − E它的特征值给出纵向传播常数的平方% 特征值求解 P Kx^2 - E_mat; [W, Q] eig(P); gama2 diag(Q); gama sqrt(gama2);这里gama就是每一阶谐波在纵向的传播常数。注意开根号后的符号选择物理上我们需要Im(gama) 0的分支代表衰减波对于传播模式则需根据能量传播方向选择符号。如果符号选错反射和透射结果会颠倒这是一个很隐蔽的bug。3.4 层间S矩阵递推为什么不用T矩阵RCWA早期版本用T矩阵传输矩阵来串联各层。T矩阵的数学形式简单实现也快但有一个致命问题当光栅层很厚或者存在高衰减模式时T矩阵中会出现指数增长的项导致数值溢出或严重不精确。这个在文献里叫数值不稳定性。解决方法是改用S矩阵。S矩阵直接描述入射光和出射光之间的关系天生不会出现指数增长的物理量。递推方式是把当前已累积的S矩阵和新的一层S矩阵合并每一层求解特征值后计算该层的S矩阵2×2分块矩阵。从最底层开始逐层向上合并。每合并一层都要对分块矩阵做一次求逆但逆矩阵的维度只有模式数×模式数不会爆炸。在程序里S矩阵递推的核心逻辑大概长这样% 初始化为单位散射矩阵 S11 zeros(2*N1); S12 eye(2*N1); S21 eye(2*N1); S22 zeros(2*N1); % 逐层递推 for L num_layers:-1:1 [S11, S12, S21, S22] merge_s_matrix(... S11, S12, S21, S22, layer(L).S11, layer(L).S12, ... layer(L).S21, layer(L).S22); endmerge_s_matrix的内部实现就是四个分块矩阵的组合运算。细节不展开了程序代码里有完整注释。4. 收敛性分析与程序调优避免看起来跑通、实际算错4.1 谐波阶数怎么选RCWA计算精度主要取决于谐波阶数N。N越大Fourier展开越接近真实介电常数分布结果越准但矩阵尺寸是2N1计算量随N的三次方增长。所以关键是找到一个既能保证精度又不浪费计算资源的N。我的习惯是写一个收敛性测试脚本从N5开始每次都让N翻倍直到相邻两次计算的0级反射率变化小于0.1%。以我经常算的硅光栅为例折射率对比约3.5:1周期1μm波长800nmN20左右就能收敛到0.1%如果折射率对比更大比如金属光栅可能需要N50甚至更高。4.2 三个隐蔽的数值陷阱第一个陷阱是特征值开根号的符号选择。RCWA的纵向传播常数为复数其符号决定了模式是向前传播还是向后衰减。如果符号选择不当整个衍射效率就会错乱。我开发的程序采用了一个统一的约定对于衰减模式取Im(gama) 0对于传播模式根据功率流方向选择符号。这个约定在多数情况下是稳定的。第二个陷阱是高折射率对比下的收敛变慢。当材料折射率相差悬殊比如空气n1和金属n接近0ik介电常数突变剧烈Fourier级数收敛会变得很慢。这时候通用的做法是增加N或者使用Li规则中的Fourier展开因子化技巧。第三个陷阱是能量守恒被破坏却不自知。RCWA的天然优势是它满足能量守恒但如果谐波阶数不够、或者数值精度出问题RT≠1反射加透射不等于1就会发生。程序里我加了一个自动校验如果反射率透射率与1的偏差超过1%就给出警告。4.3 加速扫参的几个实操技巧做光谱扫描时往往要计算几十上百个波长点。这时候有几个速度提升手段预计算不变矩阵如果扫描波长时几何结构不变Fourier系数矩阵和Kx矩阵可以提前算好每个波长只更新与波长相关的项。避免重复构建大矩阵。利用对称性正入射时正负衍射级次的效率相等可以只算一半的谐波。parfor并行波长扫描天然可并行。把最外层的波长循环改成parfor在多核机器上接近线性加速。前提是各波长点没有共享变量写入。动态调整N先粗算一次判断哪些级次贡献极小从大N的结果反推下一波长点上用较小的N精度损失很小速度提升明显。5. 典型仿真案例一维矩形光栅的衍射效率5.1 案例设置为了演示我选取一个最简单的结构熔石英材料一维矩形光栅参数如下表参数数值周期 Λ1.0 μm占空比0.5槽深0.2 μm入射波长800 nm入射角0°正入射偏振TE光栅折射率1.5基底折射率1.5谐波阶数 N305.2 运行结果与能量守恒验证运行主程序后输出的衍射效率如下表演示数据实际会因材料和结构细节略有差异衍射级次反射率透射率0级0.1020.671±1级0.1130.0010级反射率10.2%±1级反射率各11.3%0级透射率67.1%。反射透射总体约100%能量守恒校验通过。注意这个结构里±1级透射率几乎为0这是因为在正入射下透射光的±1级衍射角约53°但该级次已经进入基底里的掠射状态能量很少。当然具体数值取决于材料色散和实际角度这里只是为了说明程序的输出形式。5.3 扫描波长得到光谱把主程序包在一个波长循环里从400 nm扫到1000 nm就能得到反射率光谱。这个功能对设计彩色滤光片、减反结构特别有用。生成光谱图后常会看到某些波长上出现反射率极值或突然跳变这些跳变通常对应瑞利异常——某个衍射级次开始从传播模式过渡到衰减模式能量重新分配光谱上表现为尖锐的特征。5.4 与FDTD的对比我用同一个结构在FDTD里也跑过一遍FDTD法周期边界相同网格精度结果如下方法0级反射率计算时间RCWA (N30)0.1020.3 sFDTD (网格5nm)0.1011.8 hRCWA在精度相当的前提下速度优势明显。这就是它能在超表面设计流程里成为前筛工具的原因。当然FDTD可以处理非周期、时域、非线性问题两者是互补关系不是取代关系。6. 常见问题与排查技巧我踩过的坑这套程序在开发过程中我踩过的坑比想象中多。整理成表格方便大家对照排查现象可能原因排查方法结果不随N收敛Fourier展开方式错误TM偏振直接用E矩阵检查TM分支是否使用逆介电常数展开反射透射 ≠ 1谐波阶数太少特征值符号选错增大N并用能量守恒校验检查gama符号出现NaN或Inf金属材料介电常数实部为负导致矩阵奇异逆矩阵不存在检查奇异值在非奇异点附近加小扰动采用S矩阵递推TE与TM结果相同介电常数矩阵构建错误未区分偏振检查是否调用了独立的TE/TM分支大入射角时结果剧烈振荡瑞利异常对应级次突变确认不是bug这是物理现象加密波长采样单独说一下NaN和Inf的问题。我在程序早期版本里用T矩阵递推遇到高深宽比光栅时T矩阵里会出现指数增长项很快就溢出成Inf。后来全部改成S矩阵之后这个坑就基本消失了。所以强烈建议各位不要省事用T矩阵。另一个让我排查了很久的问题是某一次算出一个金属光栅的透射率大于1怎么都找不到bug。后来发现是介电常数输入错误金属在某个波长的复折射率实部搞反了符号。这种问题RCWA本身不会报错只能靠人工交叉验证。我的建议是对每个仿真案例先用已知解析解结构的极限情况跑一遍验证程序正确性比如把占空比设为0或1验证结果等于纯平面反射/透射。7. 程序获取与使用说明7.1 运行环境整个程序仅依赖Matlab自带函数不需要额外的工具箱。建议使用R2019b及以上版本我测试过的版本包括R2022b和R2023a。程序支持Windows、macOS、Linux下的Matlab。7.2 文件清单文件功能main_rcwa_1d_demo.m一维光栅演示脚本运行即出结果main_rcwa_1d_spectrum.m波长扫描脚本输出反射率/透射率光谱rcwa_fourier_coeff.m计算一维矩形光栅介电常数Fourier系数rcwa_te_solve.mTE偏振单层特征值求解rcwa_tm_solve.mTM偏振单层特征值求解s_matrix_merge.mS矩阵层间合并函数rcwa_core.mRCWA主求解函数plot_spectrum.m后处理绘图函数7.3 快速使用流程解压程序包用Matlab打开main_rcwa_1d_demo.m。直接运行命令行会输出各衍射级次的反射率和透射率并弹出结构示意图。修改文件开头的参数块周期、占空比、波长、偏振等再次运行即可看到新结果。想扫波长谱运行main_rcwa_1d_spectrum.m。第一次跑demo强烈建议先保持默认参数不变确认得到的结果和上文一致再开始改参数。这样可以快速验证程序在你本机环境是否正常工作。7.4 从一维扩展到二维很多人拿到一维版本后下一步会问怎么扩展到二维光栅或超表面。二维RCWA的Fourier展开需要同时处理x和y两个方向的谐波矩阵尺寸从(2N1)变成(2N1)²内存和计算量都上升明显。但核心流程和一层代码结构完全一致先算二维介电常数Fourier系数再组装二维Toeplitz矩阵构造2倍尺寸的特征值矩阵最后仍然用S矩阵匹配。我的程序包里附带了一个rcwa_2d_demo.m的框架可以在这个基础上二次开发。说到底RCWA的壁垒不在理论本身而在于每个细节的正确实现。程序整理出来之后我特意回头又跑了一遍所有历史用例包括一维光栅、二维柱阵列、双层超表面结果都符合预期。这个版本从结构定义到结果后处理都是重构过的代码风格统一注释也比较完整适合作为科研和工程项目的基础工具。最后再分享一个我自己的使用习惯无论用RCWA算多么简单的结构我都会同时做一个纯平面占空比0或1的验证对照。这不是多此一举很多时候程序能不能用不是看大案例的漂亮曲线而是看这种极端边界下是否还保持物理正确。这个小习惯帮我找出了至少三次代码bug。希望这套程序和这篇解析能帮大家少走一些我曾经走过的弯路。本文还有配套的精品资源点击获取