简介：面向光学与电磁学领域研究者的MATLAB实现资源包，聚焦严格耦合波分析（RCWA）与平面波展开法在周期性介质结构中的数值模拟。资源适合正在学习衍射光学、光栅设计或光子晶体分析的学生与工程师，可帮助理解从网格设置、PML边界处理到透反射计算的核心流程。资源包共24个文件，以13个m源码文件为主体，配合8个dat数据文件及说明文档（txt、md）组成，压缩包约17KB，结构紧凑，便于对照学习。已有82人学习。通过阅读RCWA.m、setupGrid.m、convmat.m等关键程序，读者能掌握傅里叶展开阶次选择、耦合波方程组构建与求解的实现细节；同时，附带的TransferMatrixMethod与RCWA-master目录为不同方法对比提供了直接范例。整套代码可作为光学数值仿真入门与进阶的参考资料，也可作为开发自定义计算工具的起点。 做光子晶体和亚波长光栅仿真的人，大概率绕不开两个名字：严格耦合波分析（RCWA）和平面波展开法（PWEM）。我第一次用MATLAB把这两套方法落地，是为了算一个硅波导上的周期光栅耦合器——既想知道正入射光在哪个波长会被反射掉多少，又想知道周期性结构本身的色散关系长什么样。RCWA帮我拿到了各级衍射效率，PWEM帮我画出了带隙位置。两个方法一个解决“周期结构对外如何散射”，一个解决“周期介质内部支持什么模式”，合起来基本覆盖了周期微纳结构仿真的大部分需求。
这篇文章把我自己在MATLAB里跑通整个流程的思路、核心代码片段和踩过的坑整理出来。不是完整工具箱，是能让你从零开始理解并动手实现两套数值方法的一条清晰路径。适合正在学光子晶体、亚波长光栅、超表面或者周期性滤波器设计的人参考，也适合想搞懂数值算法而不是只会点软件按钮的同学。
1. 先分清两个方法在解什么问题
很多教材把RCWA和PWEM放在一起讲，但初学者最容易迷糊的是：这两个东西明明都叫“展开法”，到底有什么区别？什么时候用哪个？
1.1 严格耦合波分析：周期性结构对外的“反射/透射”
RCWA解决的是“一束平面波打到周期性结构上，反射和透射各是多少”。它把结构的介电常数按空间周期做傅里叶展开，把电磁场也按衍射级次展开，然后通过每一层内部的模式场匹配，把麦克斯韦方程组变成一组线性方程。最终输出的是各衍射级的反射率R_m和透射率T_m，有了这些就能画反射谱、透射谱、衍射效率，还能看能量守恒是否满足。
RCWA在工程上最常用的场景是光栅、光子晶体平板、超表面单元结构。它的核心假设是结构在一个或两个方向上严格周期，但在第三个方向上可以分层。换句话说，RCWA把任意复杂截面的周期结构沿纵深切成多层，每一层被当作均匀的周期介质，再逐层传递边界条件。
1.2 平面波展开法：周期介质里到底藏着哪些模式
PWEM解决的是“无限大周期结构里，对应某个波矢k，允许哪些频率的模式存在”。它的出发点是布洛赫定理：在周期介质中，电磁场可以写成布洛赫波形式，也就是平面波包络乘以周期函数。把这个周期函数再按倒格子矢量展开，麦克斯韦方程组就变成了一个本征值问题。求解这个本征值问题，就得到了频带结构（能带图）。
PWEM的典型输出是ω-k关系，也就是能带图。从能带图上你能直观地看到带隙在哪儿、慢光区域在哪儿、某个频率的光能不能在结构里传播。它是设计光子晶体、研究带隙结构最常用的分析工具。
1.3 为什么我把两者都放在MATLAB里实现
商业电磁仿真软件（比如FDTD、FEM类工具）当然能做这两个计算，但它们的参数控制逻辑往往封装在GUI里，改算法细节很麻烦。RCWA和PWEM的核心计算本质都是矩阵构建加特征值/线性方程组求解，MATLAB在这件事上有天然优势：矩阵操作语法简洁、eig指令强大、画图方便。
更重要的是，自己用MATLAB实现会让你真正理解每一个数值步骤的含义。比如截断阶数为什么取奇数、为什么TM偏振的介电常数矩阵要用逆介电常数做Toeplitz、为什么高频能带需要更多平面波才能收敛——这些知识在工程调整参数时非常有用，不是给个结果就完了的商业软件能教的。
2. 写代码前的数学准备和参数选型
代码能跑起来之前，有两个参数选型你必须搞清楚，否则后面排查问题时会很痛苦。这两个参数分别对应两个方法各自的“截断”操作。
2.1 级次截断和收敛性之间的关系
RCWA里我们要保留多少个衍射级次？这句话听起来简单，实际是RCWA最核心的数值技巧。衍射级次理论上从负无穷到正无穷，但实际计算只能截断。截断级数N直接决定计算精度和速度。
我自己的经验是：对于常见的硅光栅，折射率对比在3.5比1左右的场景，N取21到41（即-10到10或-20到20）基本能收敛；如果结构里有金属或者折射率对比特别大，N可能需要取到81甚至更多。截断阶数太少会导致反射率出现负值或大于1，这是典型的高阶衍射贡献被强行砍掉造成的伪影。
还有一个容易忽略的细节：N一定要取奇数。因为衍射级次是在0级两侧对称分布的，取奇数才能保证级次集合是-（N-1）/2到（N-1）/2的完整对称集合。如果取偶数，某个一级衍射要么多了正侧要么多了负侧，能量分配就会不守恒。
PWEM的截断逻辑类似，只是截断的是倒格子矢量集合。一个常见做法是以原点为中心，截取|G|<G_max范围内的所有倒格矢。G_max越大，低频带越收敛，高频带则需要更多的平面波才能让结果可信。
2.2 布里渊区路径：能带图的横轴怎么打点
PWEM算能带图时，横轴是倒空间里的波矢k。通常的做法不是遍历整个第一布里渊区，而是沿着高对称点连线扫描。因为能带的极值一般出现在高对称点上，扫描Γ-X-M-Γ这种路径就能把握带隙的全貌。
扫描时需要注意分辨率问题。每个高对称段之间至少要取20到30个点，否则带边位置看不清楚。如果要做慢光或者群速度分析，点数还要加密。我在实际代码里一般用linspace在每段路径上取50个点，既不会让矩阵反复求解太多次，又能看清能带细节。
2.3 结构参数怎么设置，算出来的结果才好看
RCWA和PWEM的输入参数有重叠也有差异。比如周期、填充率、材料折射率都是公共参数，但RCWA还需要入射角度、偏振和结构深度，PWEM则需要确定晶格类型和高对称路径。
参数设不好最直接的后果是结果“看起来很奇怪”。比如RCWA算反射谱，在短波段看到周期性振荡，大概率是结构深度太大、模式干涉效应叠加所致，这不一定是代码bug，而是真实的Fabry-Perot共振。PWEM算能带图，如果带隙出现在太高的频率，说明填充率或者折射率对比选得不合适。
一个快速判断参数合理性的方法是：先用极值情况做验证。RCWA里把光栅深度设为0或让n1=n2，应该得到反射率接近0的结果；PWEM里把介电常数对比去掉（ε_a=ε_b），能带应该退化成自由空间色散直线ω=ck。如果这两个最简单的极限测试都过不了，说明代码里有问题，先不要碰复杂结构。
3. 一维矩形光栅RCWA的MATLAB落地过程
RCWA的完整实现代码量不算小，但核心思想非常清晰：第一步，把介电常数展开成傅里叶级数；第二步，构建本征矩阵并求解传播常数；第三步，用边界条件拼接各层场；第四步，提取衍射效率。下面按步骤拆开讲我用到的写法。
3.1 从介电常数到Toeplitz矩阵
以最常见的一维矩形光栅为例，结构在x方向周期排列，入射面是x-z平面，光栅脊宽为fill*Λ。介电常数ε(x)是x的周期函数，把它展开成傅里叶级数时，矩形光栅的展开系数有解析表达式，不需要做FFT数值积分，直接写公式即可。
% 一维矩形光栅RCWA，TE偏振（电场沿y方向） lambda = 1.0; % 自由空间波长，单位um Lambda = 0.8; % 光栅周期 theta = 30; % 入射角，单位度 h = 0.5; % 光栅深度 fill = 0.5; % 占空比：脊宽 / 周期 n1 = 1.0; % 入射区折射率 n2 = 1.5; % 光栅介质折射率 Norder = 21; % 截断级次，必须取奇数 % 衍射级次索引 mid = (Norder+1)/2; m = (-(Norder-1)/2:(Norder-1)/2).'; % 介电常数傅里叶展开 eps1 = n1^2; eps2 = n2^2; f = fill; eps_coeff = f*eps2 + (1-f)*eps1; % 占位 for idx = 1:Norder if m(idx) == 0 eps_coeff(idx) = f*eps2 + (1-f)*eps1; else eps_coeff(idx) = (eps2 - eps1) * sin(pi*f*m(idx)) / (pi*m(idx)); end end % 注意：用 Toeplitz 矩阵表达介电常数在倒空间的卷积 Emat = toeplitz(eps_coeff);这段代码里最值得留意的是toeplitz(eps_coeff)。RCWA里介电常数在倒空间的卷积运算，恰好等价于一个Toeplitz矩阵左乘场展开系数向量。利用MATLAB的toeplitz指令可以三行写完，手写循环反而容易出错。
3.2 本征求解与层内场展开
有了介电常数矩阵，紧接着构建无量纲化的波矢矩阵和本征矩阵。TE偏振的公式相对简单：Q = Kx² - E，其中Kx是各衍射级横向波矢组成的对角矩阵，E就是刚才的Toeplitz矩阵。
% 各衍射级横向波矢（归一化到k0） k0 = 2*pi/lambda; kx = k0 * ( n1*sind(theta) - m*lambda/Lambda ); Kx = diag(kx / k0); % TE偏振的本征矩阵 Q = Kx*Kx - Emat; % 求解本征值 q^2 和本征向量 W [W, q2] = eig(Q); q = sqrt(diag(q2));到这里已经拿到了光栅层内部的模式解：W是各模式的场分布本征向量，q是归一化传播常数。这里的q可能出现纯虚数或复数，对应的是渐逝模。物理上要求模式在传播方向上衰减而不增长，所以取本征值时要保证虚部符号正确，否则边界条件组装后会出现数值发散。
TM偏振的差异在于介电常数矩阵的处理。按照Li规则，TM偏振下应该用逆介电常数1/ε的傅里叶展开系数做Toeplitz矩阵，再取逆，而不是直接把ε的Toeplitz矩阵求逆。这个细节很多早期论文都踩过坑。我之前因为偷懒直接用了Emat的逆，结果TM偏振的反射谱在高阶衍射阈值附近出现明显的数值振荡，换成1/ε展开之后立刻干净了。
从本征解到反射率/透射率的最后一步，是用入射区和出射区的平面波展开，把电场和磁场的切向分量在z=0和z=h两个界面上做匹配。这一步本质上是组装一个稠密线性方程组，方程组右端是入射波在0级上的归一化振幅。解出反射系数R_m和透射系数T_m之后，衍射效率由如下关系给出：
- 反射效率：
DE_r(m) = R_m² · Re(γ_I,m) / γ_I,0 - 透射效率：
DE_t(m) = T_m² · Re(γ_III,m) / γ_I,0
其中γ是各衍射级在入射区/出射区的纵向波矢分量。注意这里分子分母都要用纵向波矢的实部，因为渐逝模虽然在场展开中存在，但不携带能量，不能算进能量守恒。
3.3 算衍射效率的时候最容易犯的错
我自己在写这一步时踩过两个坑，说出来给大家提个醒。
第一个坑是波矢符号。很多教科书在推导时采用时谐因子exp(iωt)，另一些用exp(-iωt)。如果混着抄公式，横向波矢的符号会整体反号，结果看起来似乎对称，但斜入射时某几个衍射级的能量就完全错了。我的做法是统一用exp(-iωt)，所有公式从S矩阵推导开始就保持一致，中途不混用。
第二个坑是渐逝模的能量计算。前面提到衍射级进入渐逝状态时，纵向波矢为虚数，此时反射率公式里的Re(γ)为0，能量贡献为0。但如果你忘了取实部而是直接用模值，就会出现R+T明显大于1的异常结果。这种bug特别隐蔽，因为它不会报错，画出来的反射谱也大体合理，只有检查能量守恒时才会暴露。
4. 二维光子晶体PWEM的实现与能带计算
PWEM的实现比RCWA更简洁，核心就三步：建倒格矢、构矩阵、扫K点。这里我以一维光子晶体为例给出完整可运行代码，再说明二维扩展的做法，这样代码逻辑最清晰，也不容易踩坐标系的坑。
4.1 倒格子基矢与平面波集合生成
一维光子晶体周期为a，介质A的占空比为f，介电常数ε_a，背景ε_b。倒格子矢量G = m * (2π/a)。平面波截断集合就取m从-Ng到Ng的所有整数。这个集合的大小直接决定矩阵尺寸。Ng取8对应17个平面波，能带低频部分基本够看；要精确算高能带或带隙边界，Ng至少取15以上。
% 一维光子晶体的平面波展开法 a = 1.0; % 晶格常数（归一化） f = 0.3; % 介质A占空比 eps_a = 11.0; % 介质A介电常数（高折射率） eps_b = 1.0; % 背景介电常数 Ng = 8; % 每个方向截断范围 G = (-Ng:Ng)' * (2*pi/a); N = length(G); % 逆介电常数1/eps的傅里叶展开系数 inv_eps = zeros(N,1); inv_eps(mid_index) = f/eps_a + (1-f)/eps_b; mm = (-Ng:Ng)' ; nonzero = mm ~= 0; inv_eps(nonzero) = (1/eps_a - 1/eps_b) .* sin(pi*f*mm(nonzero)) ./ (pi*mm(nonzero));这里最核心的细节是：展开的是逆介电常数，不是介电常数本身。原因在于TM偏振下我们求解的方程是：
∇ × (1/ε) ∇ × E = (ω/c)² E
算符里自带1/ε，所以展开它才是符合物理的。如果误把ε的傅里叶展开拿来用，会在介电常数突变界面处产生严重的Gibbs振荡，带隙位置算不准。
4.2 哈密顿量矩阵构建的关键写法
一维情况下本征方程退化成标量积分方程：
∑_G' κ_{G-G'} (k+G)(k+G') E_{G'} = (ω/c)² E_G
矩阵元H_{G,G'} = (k+G)(k+G') * κ_{G-G'}。这个矩阵是稠密的，直接用MATLAB的eig求全谱即可。对于一维问题，平面波数量通常不超过31个，矩阵也就31×31，求解速度可以忽略不计。
% 中间代码：以k=0.3*(2*pi/a)为例构建矩阵 k = 0.3 * 2*pi/a; H = zeros(N,N); for i = 1:N for j = 1:N idx_diff = (i - j) + (Ng+1); % G_i - G_j 在 inv_eps 中的索引 H(i,j) = (k + G(i)) * (k + G(j)) * inv_eps(idx_diff); end end omega_sq = eig(H); omega = sqrt(max(0, real(omega_sq))); % 频率取正根构建索引idx_diff时，我直接用i-j的偏移量换算到inv_eps数组的索引。因为G_i和G_j的差依然在截断集合内（一维截断集合是完整的），这样做准确且高效。二维情况略复杂，G_i-G_j可能超出截断范围，超出就置0，不影响低频结果的正确性。
4.3 一维光子晶体案例：带隙在哪里
下面这个例子结构是空气孔型光子晶体的一维对应：高折射率层厚度占比30%，背景为空气。扫描k从0到π/a，取每个k点的前几条频率解，就能画出标准的ω-k图。
% 扫描第一布里渊区边界 0 到 pi/a klist = linspace(0, pi/a, 60); bands = zeros(N, length(klist)); for nk = 1:length(klist) k = klist(nk); H = zeros(N,N); for i = 1:N for j = 1:N idx_diff = (i - j) + (Ng+1); H(i,j) = (k + G(i)) * (k + G(j)) * inv_eps(idx_diff); end end omega_sq = sort(eig(H)); % 注意排序 bands(:, nk) = sqrt(max(0, real(omega_sq))); end % 画图：频率作为色散纵轴，只显示前几条 figure; plot(klist/(2*pi/a), bands(1:6,:)*a/(2*pi), 'b-', 'LineWidth', 1.2); xlabel('波矢 k (2\pi/a)'); ylabel('归一化频率 \omega a / 2\pi c');跑完这个代码，大概在归一化频率0.2附近能看到明显的带隙。如果ε_a和ε_b的对比越大，带隙越宽；占空比f越接近0.5，带隙宽度一般也越大。这是设计一维光子晶体反射镜最基本的结论，拿这个脚本调参数非常直观。
扩展到二维时，改动点其实只有三处：倒格矢变成二维向量、矩阵元里的(k+G)(k+G')换成二维点乘、逆介电常数傅里叶系数换成圆域/方域的结构因子。其他逻辑完全一致。我第一次写二维三角晶格能带，就是把一维代码里的标量改成向量，半小时就调通了。
5. 常见问题与排查技巧速查
两个方法跑起来之后，一定会遇到各种“看起来像bug但其实是数值问题”的情况。我整理了一份速查表，都是我实际调试过程中遇到过的典型问题。
5.1 收敛性差：先怀疑这两个参数
RCWA反射谱出现负值、或R+T总和远离1，先别急着查边界条件，大概率是截断级次不够。尤其当入射角接近衍射级次“临界角”时，某个级次正处于传播到渐逝的临界状态，截断级次少一个都会引起明显的能量泄漏。
PWEM能带低频部分收敛快、高频部分慢，这是物理事实。但如果最低几条能带在同一个k点处出现明显的折线转弯，通常是平面波数量太少，导致带边位置偏移。最简单的方法是把Ng从8加到16，看目标频段结果是否变化明显——如果变化超过1%，继续加。
5.2 “伪模式”与数值杂波的识别
PWEM能带图里有时会出现一些孤立点或很陡的“鬼线”，这些不是真实物理模式，而是数值伪模式。常见来源有两个：一是矩阵存在零特征值附近的小扰动，平方根后产生虚假低频模式；二是在某些k点上出现了简并模式的数值分裂。
我的排查思路是：把能带数据按k排序后，看这些异常点是否在其他k点呈连续分布。如果只是孤立的、不连续的跳点，就先用后处理把它们剔除；如果连续出现，则要怀疑是Ng截断不够或者介电常数展开错误。这个排查过程没有捷径，只能多看几个k点、多对比不同截断下的结果。
还有一个实用技巧：RCWA计算完成后，把各衍射级效率数组重新相加，务必做到|R+T-1|<1e-8。如果这个容差达不到，代码里一定有bug，不要试图用更大的截断去掩盖。
5.3 几个用着顺手的MATLAB写法
最后分享几个我实际用下来效率很高的MATLAB小技巧。
eig对大矩阵是稠密求解，PWEM矩阵如果超过500×500会很慢。这时候可以改用eigs只求最低的10条能带，速度能快一个数量级。但注意eigs对参数比较敏感，需要给出合理的初始向量和容差，否则可能出现漏解。
RCWA里构建Toeplitz矩阵时，尽量用toeplitz(eps_coeff)而不是循环赋值。TOEPLITZ指令是内置的，支持向量化，几百阶矩阵瞬间生成。我自己早期的循环版本跑一次反射谱需要几十秒，改成Toeplitz后几乎可以实时扫描参数。
另外，所有涉及波矢和纵向传播常数的计算，最好统一归一化到k0，即所有波矢都除以k0再参与矩阵构建。这样做的优势是矩阵元都在0.1到100这个量级，数值稳定性好，不容易出现小数目除大数目导致的信息丢失。这个问题在深度很大或者入射角接近掠入射时特别明显。
最后的实操体会
我自己做完这套MATLAB实现之后最深的感触是：RCWA和PWEM虽然数学形式不同，但都是“周期性假设+傅里叶展开+截断求解”这个框架下的产物。把一维问题彻底想清楚，二维三维只是向量维度增加，不会引入本质困难。
如果现在有人让我快速上手这两个方法，我会建议他先跑通一维PWEM的能带代码，再做一维RCWA的反射谱。两个代码加起来不到200行，但跑通之后再去看文献里的S矩阵、Li规则、收敛性分析这些话题，就有了具象的锚点。数值方法这东西，光看书很容易停在概念层面，真正把矩阵构造、特征求解、效率提取这几个环节亲手敲一遍，才算真正掌握。
本文还有配套的精品资源，点击获取
