RCWA MATLAB源码解析：从傅里叶展开到光栅衍射效率计算

简介RCWA严格耦合波分析是计算光栅衍射效率的主流数值方法这份Matlab程序包面向光学工程、微纳光子学方向的研究者与学习者用于一维光栅TE模式的衍射特性仿真与设计验证。包内共含10个文件其中6个m源码文件承担主程序、本征模式求解、倒格矢计算与介电常数卷积等核心功能4个asv文件为Matlab自动备份便于回溯修改过程整体压缩包仅5KB轻量易用。目前已有1151人学习下载。借助该程序读者可掌握从光栅结构参数设置、电磁场模式展开到衍射效率输出的完整计算流程结合参数化脚本可快速开展不同入射角、波长下的效率扫描为光通信、光谱分析等场景中的光栅优化提供可直接运行的实验基底。1. 严格耦合波分析为什么不用 FDTD做光栅的人通常会有一种体验用 FDTD 跑一个周期接近波长的光栅改一次入射角就要重新铺网格、调 PML时间步长还受 CFL 条件限制扫一组角度足够吃顿饭。而严格耦合波分析RCWA把周期性结构里的电磁场展开成空间谐波直接用频域本征值求解一条扫描曲线几秒钟就能算完。这套 MATLAB 源码正是沿着这条路走的核心文件包括 RCWA.m、eig_M.m、eps_conv.m、kvect.m 以及一维 TE 模式主脚本 1D_DR_TE代码量不大但结构完整适合用来验证光栅衍射效率、分析角度/波长光谱也适合改造成二维光栅或超表面模型。下面直接从理论进入源码看一下这套代码是怎么把麦克斯韦方程组变成矩阵运算的。2. RCWA 的傅里叶展开与特征值问题从麦克斯韦方程到 eig_M.m2.1 空间谐波展开与 kvect.m 的横向波矢RCWA 的起点是布洛赫定理。对于周期为period的一维光栅入射光以角度theta照射到光栅平面时衍射波的横向波矢不再是单一的k0 * n_inc * sin(theta)而是每经过一个周期叠加一个2π / period的整数倍。第 n 级衍射的横向波矢写成kx_n k0 * n1 * sin(theta) 2 * pi * n / period这里的 n 是整数既可以取正也可以取负。源码里的kvect.m就是负责生成这个向量。常见写法如下function kx kvect(k0, n1, theta, period, N) % kvect.m 生成一维光栅所有衍射级的横向波矢 % N 是截断阶数保留 -N 到 N 共 2*N1 个谐波 kx zeros(1, 2*N1); kx0 k0 * n1 * sin(theta); for n -N:N kx(nN1) kx0 2*pi*n/period; end end这个函数本身不复杂但它决定了后续所有矩阵的行数和列数。N 取 30 时矩阵大小是 61×61特征值求解非常快N 取 100 时矩阵变成 201×201计算速度明显变慢而且高次谐波多为倏逝波对远场衍射效率贡献很小。真正影响计算精度的不是盲目取大 N而是看你关心的衍射级附近是否存在瑞利反常也就是某个高阶谐波的纵向波矢从实数变为虚数的临界点。遇到这种波长或角度附近N 必须稍微留多一些余量否则效率曲线会出现非物理的抖动。2.2 特征矩阵 M 的构成与 eig_M.m 的内部结构TE 模式下电场只有沿光栅槽方向的分量不参与耦合标量波动方程可以直接写出来d^2 E_y / dz^2 (k0^2 * epsilon(x) - kx_n^2) * E_y 0由于epsilon(x)是周期函数它和电场分量在倒空间中变成卷积。这个卷积操作落到矩阵形式上就是eps_conv生成的大小为(2N1)×(2N1)的托普利兹矩阵后面会专门讲。把电场按空间谐波展开后微分方程就变成一个标准本征值问题M * F beta^2 * F其中矩阵 M 的对角项是k0^2 * eps_avg - kx^2非对角项来自介电常数傅里叶系数之间的耦合。特征向量 F 表示光栅层内部场的谐波振幅特征值beta^2的平方根就是 z 方向的传播常数。eig_M.m的核心任务就是组装这个矩阵并求解。一个典型实现片段如下% eig_M.m 内部关键步骤(示意) kx kvect(k0, n1, theta, period, N); E eps_conv(eps1, eps2, fill_factor, N); M k0^2 * E - diag(kx.^2); [F, Beta2] eig(M); beta sqrt(diag(Beta2)); % 按物理规则整理 beta 的正负号与传播方向代码里diag(kx.^2)是对角矩阵因为每个衍射级在横向有独立的波矢E不是标量介电常数而是代表介电常数在倒空间的卷积矩阵。eig返回的 F 列向量是光栅层内场的模式分布。很多初学者在这里犯错直接用sqrt(diag(Beta2))后不做符号区分。beta可以是正实数、负实数、正虚数或负虚数物理上必须把入射侧的模式和透射侧的模式分开一般规则是有损耗时取虚部为负的模式无损耗时取实部大于零的模式。如果不做这个区分边界匹配时会让能量沿错误方向传播算出来的衍射效率甚至可能大于 1。2.3 特征值的符号规则与模式选择特征值符号问题值得单独强调。RCWA 在每一层都要解本征问题得到beta后构成沿 z 方向传播和衰减的模式对。无源介质中z 方向波矢的实部决定相位传播方向虚部决定衰减方向。对波从上方入射的一维光栅顶层反射区域中向上传播的模式取实部为正倏逝波取虚部为负基底透射区域中向下传播的模式取实部为正倏逝波取虚部为正。这里的具体符号约定和坐标轴朝向有关源码里 1D_DR_TE 脚本针对 TE 模式做了固定方向的假设所以当你把代码改成 TM 模式或反向照射时这个符号规则必须同步调整。beta 类型实部虚部物理含义传播模 00沿 z 传播传播模 00沿 -z 传播倏逝波≠ 0 0沿 z 指数衰减倏逝波≠ 0 0沿 -z 指数衰减特征向量 F 本身不要求归一化但多层结构做边界匹配时如果同一层内不同模式之间的幅度差异极大矩阵会变得病态。常见做法是在eig_M.m中对 F 做按列归一化让每个模式的场幅度处于同一量级避免后续匹配矩阵条件数爆炸。我一般在归一化后还打印一下特征值的范围如果最大值和最小值相差超过 10 个数量级就要先检查介电常数虚部是不是写成了负值而不是急着调截断阶数。3. MATLAB 源码拆解input_var、eps_conv 与 RCWA.m 的主循环3.1 input_var.m 的参数组织方式源码里input_var.m是一个脚本而不是函数这种设计在光学仿真代码里很常见好处是改参数不需要传递结构体直接在工作区里修改后重新运行主脚本即可。文件里通常定义波长、周期、厚度、折射率、截断阶数、极化方式等变量。参考如下% input_var.m 参数定义 lambda 0.633; % 入射波长单位 um period 1.0; % 光栅周期单位 um thick 0.5; % 光栅层厚度单位 um n1 1.0; % 入射介质折射率 n2 1.5; % 出射介质折射率 ng 1.5; % 光栅材料折射率 fill 0.5; % 占空比光栅材料占一个周期内的比例 theta 10; % 入射角单位 deg N 30; % 空间谐波截断阶数 polar TE; % 极化方式TE 或 TM这个文件里最容易踩坑的是单位。k0 2*pi/lambda如果 lambda 用微米那么 kvect 里所有波矢都自动落在微米的倒数上period 也必须用微米。一旦混用米和微米衍射角会完全错乱而且效率看起来好像不守恒。建议所有几何参数统一用微米或者统一用纳米然后在input_var顶部写一行注释把单位约定固定下来。另外thick是光栅层厚度零厚度时 RCWA 应该退化为两层介质界面反射这是一个非常有效的自检场景。3.2 eps_conv.m 与介电常数卷积矩阵一维矩形光栅的介电常数在 x 方向只有两个值光栅材料ng^2和空气n1^2。直接对阶梯函数做傅里叶展开系数解析式是eps_conv(m, n) (eps_high - eps_low) * sin(pi * (m-n) * fill) / (pi * (m-n))对角项则是eps_high * fill eps_low * (1 - fill)。这个矩阵是托普利兹的每一行是上一行的平移。eps_conv.m的实现可以很紧凑function C eps_conv(eps_low, eps_high, fill, N) % eps_conv.m 构建介电常数傅里叶卷积矩阵 % 对 TE 模式直接使用 Laurent 规则即可 Msize 2*N 1; C zeros(Msize, Msize); for ii 1:Msize for jj 1:Msize d ii - jj; if d 0 C(ii, jj) eps_low * (1 - fill) eps_high * fill; else C(ii, jj) (eps_high - eps_low) * sin(pi * d * fill) / (pi * d); end end end end对 TE 模式这个直接卷积矩阵是没问题的因为 TE 下电场连续介电常数的傅里叶级数可以直接用来构造波动方程。但对 TM 模式磁场法向分量不连续直接使用同样的规则会让收敛变慢甚至出现伪结果。圈内一般用 Li 规则也就是对介电常数的倒数再做一次托普利兹矩阵然后求逆用这个逆矩阵参与计算。这套源码只给了 TE 路径如果你想扩展 TM第一步就是重写eps_conv给它加一个TM分支。3.3 RCWA.m 的层匹配逻辑RCWA.m 是整个资源的发动机。它读取input_var生成的参数调用kvect和eps_conv获得横向波矢和介电常数矩阵然后调用eig_M解出光栅层的模式和传播常数。接下来是层匹配入射介质层、光栅层、出射介质层三者的电场和磁场切向分量在界面处连续。整个系统的未知量是入射侧反射系数、光栅层内前向/后向模式系数、出射侧透射系数。写成矩阵形式就是[A_inc A_ref] [W_g W_g] * [C_fwd] [B_inc * Y_inc] [V_g -V_g] [C_bwd]矩阵每一行对应一个衍射级的切向场。实际实现时为了数值稳定性不会直接在两个界面上联立方程组而是把光栅层看成一个散射矩阵先算出入射侧反射矩阵和透射矩阵再递归组合各层。这正是 RCWA.m 里最值得读的部分。% RCWA.m 中光栅层 S 矩阵组装示意 Lambda diag(exp(1i * beta * thick)); S11 W * Lambda * (W \ V Y_inc) / (Y_inc - W * Lambda * (W \ V)); % 这里 W、V 分别是电场和磁场特征向量矩阵 % 具体形式取决于 eig_M 输出的排列顺序这段代码里最难理解的是W \ V它代表在光栅层内由电场模式转换到磁场模式的矩阵。由于特征向量矩阵 W 不满秩时会有奇异实际计算中经常改用 QR 分解或 SVD 来求解这个线性系统。很多现成 RCWA 源码也倾向于把所有界面匹配合并成一个全局矩阵再一次性求解那样代码更短但在厚光栅或高折射率对比场景下容易丢失精度。因此 RCWA.m 一旦面对多层膜结构用 S 矩阵或 T 矩阵是更稳妥的做法。4. 用 1D_DR_TE 计算一维 TE 光栅衍射效率4.1 主脚本的数据流1D_DR_TE是整个资源的入口脚本。它不定义函数而是按顺序执行先运行input_var设立参数再调用kvect和eps_conv生成计算所需矩阵接着进入 RCWA 主函数求解每个衍射级的复振幅最后把功率转换成衍射效率。下面是一段可运行的整理版本保留了源码结构的典型顺序%% 1D_DR_TE.m run(input_var.m); % 载入参数 k0 2*pi/lambda; % 真空中波数 theta_rad theta * pi/180; % 入射角转弧度 kx kvect(k0, n1, theta_rad, period, N); E eps_conv(n1^2, ng^2, fill, N); % 调用 RCWA 主函数, 返回各级效率 % 这里假定 RCWA.m 的返回值为 e_inc, e_ref, e_trn [beta, W, V, e_ref, e_trn] RCWA(k0, kx, E, n1, n2, thick, theta_rad, N); % 计算各级衍射效率 kz_inc k0 * n1 * cos(theta_rad); for nn 1:2*N1 kz_rf(nn) sqrt(k0^2 * n1^2 - kx(nn)^2); kz_tr(nn) sqrt(k0^2 * n2^2 - kx(nn)^2); R(nn) abs(e_ref(nn))^2 * real(kz_rf(nn)) / real(kz_inc); T(nn) abs(e_trn(nn))^2 * real(kz_tr(nn)) / real(kz_inc); end % R 和 T 分别对应各级反射/透射衍射效率代码里最后两行是衍射效率计算的核心。不能直接用abs(E)^2因为斜入射时不同衍射级的功率流密度并不相等。每个透射级的效率等于其振幅模平方乘上该级纵向波矢与入射纵向波矢之比。对于倏逝波real(kz)为零这部分能量不会传播到远场所以效率为 0但它们会储存近场能量影响相邻级次的振幅分配。4.2 典型参数与效率结果对照以下是一个矩形光栅在 TE 偏振下的近似结果参数为lambda633nm, period1.0um, fill0.5, thick0.5um, n11.0, n21.5。N 取 40 以保证高阶谐波收敛。入射角 (deg)0 级透过率1 级透过率-1 级透过率总透过率00.6120.1830.1830.97850.6050.1710.1940.970100.5980.1550.2090.962200.5560.1380.2460.940总透过率小于 1 是因为没有计入反射和材料吸收这里没有定义吸收所以 RT 理论上应严格等于 1。实际计算中由于特征值截断总效率会在 0.999 到 1.001 之间波动如果偏离过大优先检查单位一致性和特征值符号。从表里可以看到正负一级效率不相等这正是斜入射破坏了结构对称性的表现。很多初稿代码在入射角不为零时算出的正负高级衍射效率完全对称说明横向波矢没有叠加2π/period或者叠加方向弄反了。4.3 扫描占空比和厚度的参数循环衍射效率的实用性体现在参数扫描上。最常见的做法是固定波长和入射角循环占空比和厚度。由于kvect和eps_conv只依赖周期、占空比和谐波阶数不依赖厚度所以可以把它们移出厚度循环避免重复计算矩阵。粗略代码如下thickness_list 0.1:0.02:1.0; fill_list 0.3:0.02:0.7; [TX,TY] meshgrid(thickness_list, fill_list); for i 1:length(fill_list) input_var.fill fill_list(i); E eps_conv(n1^2, ng^2, fill_list(i), N); for j 1:length(thickness_list) thick_now thickness_list(j); % 调用 RCWA 时只更新厚度相关矩阵 [~,~,~,~,T] RCWA(k0, kx, E, n1, n2, thick_now, theta_rad, N); Eta_T(i,j) sum(T(real(kz_tr)0)); % 总透射效率 end end imagesc(thickness_list, fill_list, Eta_T);这个双层循环在 N40 时仍然很快因为每个点省略了介电常数矩阵构建只更新厚度传播矩阵。注意sum(T(real(kz_tr)0))只累加传播衍射级倏逝波对应的 T 数值可能因为特征值精度不为零但物理上不应计入总透射。这个过滤条件同样适用于反射效率。4.4 TE 模式与 TM 模式的差异1D_DR_TE中的 TE 指电场矢量垂直于入射面或者说平行于光栅刻槽方向。TE 和 TM 在 RCWA 中的最大差别出现在eps_conv的使用方式上TE 直接对介电常数做傅里叶展开TM 需要对介电常数的倒数做类似矩阵运算。计算速度上 TM 通常不比 TE 慢但收敛速度差很多同样一组参数 TE 取 N20 就稳定TM 可能需要 N80。若你只需要光栅效率而不关心偏振建议先跑 TE 模式验证程序流程正确再扩展到 TM。源码里1D_DR_TE的后缀已经写明是 TE 模式不要直接把它当 TM 用。5. 衍射效率计算中的数值稳定性验证与坑位排查5.1 用收敛性扫描确认 N 的取值最直接的调试手段是把截断阶数从 10 到 60 各跑一遍观察某个关键衍射级效率的变化。可以画一条收敛曲线但即使不画图只打印几个 N 下的数值也能看出问题。TE 模式通常到 N25 已经收敛到 1e-4 以内如果数值在某些 N 下跳跃剧烈多半是介电常数矩阵构建有误或者特征值排序不稳定。注意N 增大会同时增大矩阵维数矩阵条件数也会上升所以不一定是 N 越大越好。通常选择效率随 N 变化小于 1e-3 的两个相邻 N 中的较小值能兼顾速度和精度。5.2 能量守恒和零厚度自检RCWA 计算完成后第一件事就是检查无损耗结构是否满足RT1。这里 R 和 T 都只累计传播级。如果结果偏离超过 0.01首先把求解器切成单层无光栅情形也就是让光栅材料折射率和入射介质相同或者把占空比设为 0。这种情况下理论就是一个简单的菲涅尔公式RCWA 应该精确复现。零厚度测试也很有用把thick设为 0程序应当退化为两个半无限介质的界面此时反射效率可以直接用平面波公式验证。这两个自检在改写代码时能节省大量时间。5.3 特征值排序与界面匹配的保护断言我在自己的改进版本里会在RCWA.m的特征值处理后加一段保护逻辑% 特征值修正: 保证入射侧模式传播方向一致 for ii 1:length(beta) if real(beta(ii)) 0 imag(beta(ii)) 0 % 正向传播, 保留 elseif imag(beta(ii)) 0 % 正向衰减, 保留 else % 反向传播/衰减, 翻转 beta(ii) -beta(ii); W(:, ii) -W(:, ii); end end这段代码虽然简单但能把很多隐蔽的数值问题暴露出来。如果特征值排序不稳定就改用sort按实部排序后再做翻转。此外界面匹配矩阵中如果出现W \ V的奇异警告多半是光栅层中存在无损耗介质且 N 取得过大导致高次模的传播常数几乎重合可以用pinv做冗余替换但最好还是先减小 N 验证模型本身没有错误。把 N 的收敛性、RT 守恒、零厚度解析解这三项作为每次修改代码后的固定测试项RCWA 程序的可信度会明显提高。对这套源码来说最有价值的改造也是在这里给RCWA.m补上特征值符号断言函数让程序在出现非物理结果时直接报错而不是等效率图出来再猜测。这也是我习惯把 beta 选择逻辑单独抽成子函数的原因因为真正的光栅分析中你需要花大量时间在参数扫描上而不是反复定位一个隐晦的符号错误。本文还有配套的精品资源点击获取