简介：本资源是面向光学仿真与微纳光子器件设计领域的研究者及工程师的RCWA（严格耦合波分析）1D亚波长光栅建模与设计工具包，聚焦非周期性偏转/汇聚型光栅的参数化建模与性能优化。资源提供完整的MATLAB实现框架，涵盖光栅结构建模、傅里叶级数展开、特征矩阵求解、衍射效率与相位响应计算等核心模块，适用于微纳光学、超构表面预研、光纤传感及太阳能陷光结构设计等场景。压缩包共179个文件，含105个.m主程序与函数脚本（实现RCWA算法全流程）、25个.mat参数与仿真结果数据、31个.txt说明与配置文件、8个.fig可视化图例（如波长-占空比响应曲线、偏转角分布图等），另有PDF文档与license文件，整体9.14MB，结构层次分明，便于按功能模块调用与二次开发。目前已有186人学习下载，用户可直接复现亚波长光栅在不同入射角、波长、深度与占空比下的偏转/聚焦特性，并基于现有代码快速拓展至非周期梯度光栅或CCSWG（啁啾耦合亚波长光栅）设计。
1. RCWA 1D 工具包不是仿真软件，而是亚波长光栅设计的“参数化引擎”
你下载了rcwa-1d-02保留版本.zip，解压后看到一堆.m文件、README.txt和几个示例脚本——但没有图形界面，也没有“运行”按钮。这不是 bug，而是 RCWA 1D 的本质：它不渲染电场分布，也不自动生成结构图；它是一套用 MATLAB 实现的严格耦合波分析（RCWA）一维算法封装，专为快速迭代非周期亚波长光栅（sub-wavelength grating）的衍射响应而生。典型场景是：光学工程师在设计偏转型（beam-steering）或聚焦型（focusing）超构表面时，需要在几十微秒内评估上千种单元结构的 0 级与 ±1 级衍射效率、相位响应、偏振转化率。此时，调用rcwa_1d函数比启动全波电磁仿真器快两个数量级。它适合两类人：一是熟悉 MATLAB 的光学设计者，需将光栅参数（占空比、高度、材料折射率、周期序列）直接映射为衍射谱；二是算法开发者，需在其反演优化流程中嵌入高精度、低开销的正向模型。标题中“非周期偏转或者汇聚”正是其核心能力边界——它不假设光栅具有平移对称性，而是通过分段常数近似（piecewise-constant approximation）处理任意变化的单元轮廓，从而支撑 chirped、apodized、渐变周期等真实器件建模。
2. 用 rcwa_1d 在本地跑通 sub-wavelength 光栅的最小命令链
RCWA 1D 的入口函数是rcwa_1d，但它不接受原始几何描述。必须先将物理结构离散为“层-切片”模型：每一层由厚度、介电常数和横向分段数定义。整个流程分三步：构建结构描述体（structure）、设置入射条件（incident）、调用求解器（solve）。下面以一个最简非周期光栅为例——5 个单元，每个单元周期线性递减（chirped），用于实现角度偏转——展示从零开始的可复现命令链。
2.1 构建非周期光栅的 structure 结构体
非周期性体现在period字段为向量而非标量。注意：rcwa-1d-02版本要求所有单元在同一层内定义，因此需将光栅视为单层介质柱阵列，其轮廓由eps_r（相对介电常数）矩阵按 x 方向分段给出：
% 定义 5 个非周期单元，周期从 800nm 递减至 600nm（单位：nm） lambda0 = 1550; % 设计波长 1550 nm periods = [800, 750, 700, 650, 600] * 1e-9; % 转为米 height = 300e-9; % 光栅高度 300 nm n_Si = 3.48; % 硅在 1550nm 的折射率 n_air = 1.0; Ncell = 5; % 单元总数 Nslice = 32; % 每单元横向切片数（分辨率控制） % 初始化 structure 结构体 s = struct(); s.period = periods; % 关键：非周期 → 向量 s.height = height; s.Nslice = Nslice; s.eps_r = zeros(Ncell, Nslice); % 每单元每切片的 eps_r % 为每个单元填充介电常数：空气基底 + 硅柱（占空比 0.6） for i = 1:Ncell s.eps_r(i, 1:round(0.6*Nslice)) = n_Si^2; % 硅区域 s.eps_r(i, round(0.6*Nslice)+1:end) = n_air^2; % 空气区域 end % 设置背景介质（基底和覆盖层） s.eps_bkg_up = n_air^2; % 上方覆盖层（空气） s.eps_bkg_down = n_Si^2; % 下方基底（硅衬底）提示：
s.eps_r是Ncell × Nslice矩阵，行对应单元，列对应该单元内横向位置。RCWA 1D 不解析连续函数，只认这个“像素化”的介电常数分布。若要实现平滑渐变轮廓（如抛物线型汇聚光栅），需增大Nslice并用插值填充eps_r。
2.2 设置入射光与求解参数
入射条件决定衍射级次的物理意义。对于偏转/汇聚设计，必须启用斜入射（theta_i ≠ 0）并指定偏振态。rcwa-1d-02支持 TE/TM 及任意椭圆偏振，但最常用的是线偏振：
% 入射参数 inc = struct(); inc.lambda0 = lambda0; % 波长（米） inc.theta_i = 0; % 入射角（弧度），0 表示正入射 inc.phi_i = 0; % 方位角（弧度），固定为 0（1D 模型仅 y-z 平面有效） inc.pol = 'TE'; % 偏振：'TE'（电场平行于光栅条纹）或 'TM' % 求解参数 opt = struct(); opt.NG = 11; % 傅里叶阶数（必须为奇数），控制衍射级截断 opt.tol = 1e-8; % S 矩阵收敛容差 opt.maxit = 100; % 最大迭代次数 opt.verbose = 1; % 显示进度注意：
opt.NG直接决定能解析的最高衍射级。若设计目标是将能量集中到 ±1 级实现偏转，则NG=11可覆盖 -5 到 +5 级，足够冗余；但若NG过小（如 3），高阶衍射被截断，结果严重失真。经验法则是：NG ≥ 2 × ceil(最大期望衍射级 × 1.5)。
2.3 调用 rcwa_1d 并提取关键输出
执行求解后，返回值包含各衍射级的复振幅反射/透射系数，以及能量守恒验证信息：
% 执行 RCWA 计算 [refl, tran, info] = rcwa_1d(s, inc, opt); % 提取 0 级和 ±1 级的衍射效率（功率占比） R0 = abs(refl(1))^2; % 0 级反射效率 R1p = abs(refl(2))^2; % +1 级反射效率 R1m = abs(refl(3))^2; % -1 级反射效率 T0 = abs(tran(1))^2; % 0 级透射效率 T1p = abs(tran(2))^2; % +1 级透射效率 T1m = abs(tran(3))^2; % -1 级透射效率 % 验证能量守恒：R_total + T_total ≈ 1 R_total = sum(abs(refl).^2); T_total = sum(abs(tran).^2); fprintf('能量守恒误差: %.2e\n', abs(1 - (R_total + T_total)));逻辑说明：
refl和tran是长度为NG的复数向量，索引1对应 0 级，2对应 +1 级，3对应 -1 级，依此类推。rcwa-1d-02默认按k_x = k0×sin(theta_i) + m×2π/period_avg排序，但由于输入period是向量，实际采用的是等效平均周期（mean(periods)）作为基准。这意味着：非周期结构的衍射角计算需后处理校正——见第 4 章。
3. 设计非周期偏转光栅：从衍射效率到相位梯度的三步映射
标题中“非周期偏转或者汇聚”不是靠改变单个单元周期就能实现，而是通过空间变化的相位响应构造波前调控。RCWA 1D 本身不输出相位，但可通过angle(refl(m))或angle(tran(m))提取复振幅的相位。关键在于：如何将“每个单元的相位响应”与“整体偏转角”关联？答案是广义斯涅尔定律（Generalized Snell’s Law）的离散形式。
3.1 提取每个单元的局部相位响应
由于rcwa_1d对整个非周期结构一次性求解，无法直接返回每个单元的独立响应。但rcwa-1d-02提供了rcwa_1d_cell函数，可对单个单元单独建模。这是设计流程中不可或缺的预处理步骤：
% 预计算每个单元的透射相位（固定入射角 theta_i = 0，TE 偏振） phi_T = zeros(1, Ncell); eta_T = zeros(1, Ncell); % 对应透射效率 for i = 1:Ncell s_cell = struct(); s_cell.period = periods(i); s_cell.height = height; s_cell.Nslice = Nslice; s_cell.eps_r = repmat(s.eps_r(i,:), 1, Nslice); % 复制为单单元矩阵 s_cell.eps_bkg_up = n_air^2; s_cell.eps_bkg_down = n_Si^2; inc_cell = inc; inc_cell.theta_i = 0; [~, tran_cell, ~] = rcwa_1d(s_cell, inc_cell, opt); phi_T(i) = angle(tran_cell(1)); % 取 0 级透射相位 eta_T(i) = abs(tran_cell(1))^2; % 0 级透射效率 end参数说明：此处
tran_cell(1)是单单元 0 级透射复振幅。相位phi_T(i)即该单元引入的透射波前延迟。注意：rcwa_1d_cell需手动构造单单元s_cell，不能复用原s结构体。
3.2 构建相位梯度与目标偏转角的约束方程
对线性 chirped 光栅，期望相位分布为phi(x) = (2π/λ) × Δn × x，其中Δn是等效折射率差，x是沿光栅方向的位置。但更实用的是离散形式：设第i个单元中心坐标为x_i，则要求phi_T(i) ≈ phi_target(x_i)。目标相位函数由广义斯涅尔定律导出：
$$ \phi_{\text{target}}(x_i) = \frac{2\pi}{\lambda_0} \left( \sin\theta_t - \sin\theta_i \right) x_i $$
其中θ_t是目标透射角（如偏转 10°），θ_i是入射角（通常为 0）。在 MATLAB 中实现：
theta_t_target = deg2rad(10); % 目标偏转角 10 度 x_pos = cumsum([0, periods(1:end-1)/2 + periods(2:end)/2]); % 近似单元中心位置 phi_target = (2*pi/lambda0) * (sin(theta_t_target) - sin(inc.theta_i)) * x_pos; % 计算相位误差（用于后续优化） phase_error = mod(phi_T - phi_target + pi, 2*pi) - pi; % wrap to [-π, π]关键点：
phi_T和phi_target必须在同一参考系下。rcwa_1d计算的相位是相对于入射波前的绝对相位，而phi_target是理想波前所需的相位，二者可直接相减。mod(... + pi, 2*pi) - pi是标准相位解包裹操作，避免±π跳变干扰误差统计。
3.3 用参数扫描逼近最优非周期序列
rcwa-1d-02不内置优化器，但可结合 MATLAB 的fmincon或网格搜索实现自动设计。以下是最小可行方案：固定高度和材料，扫描占空比duty和周期period组合，寻找使mean(abs(phase_error))最小的序列：
% 定义搜索空间（示例：5 个单元，每个单元 duty ∈ [0.4, 0.8], period ∈ [600, 800] nm） duty_grid = linspace(0.4, 0.8, 5); period_grid = linspace(600e-9, 800e-9, 5); best_error = Inf; best_params = []; for i1 = 1:5 for i2 = 1:5 for i3 = 1:5 for i4 = 1:5 for i5 = 1:5 duties = [duty_grid(i1), duty_grid(i2), duty_grid(i3), duty_grid(i4), duty_grid(i5)]; periods_test = [period_grid(i1), period_grid(i2), period_grid(i3), period_grid(i4), period_grid(i5)]; % 构建新 structure 并计算 phi_T s_test = construct_structure(periods_test, duties, height, n_Si, n_air, Nslice); phi_T_test = compute_phi_T_per_cell(s_test, inc, opt); % 计算 phase_error phase_error_test = mod(phi_T_test - phi_target + pi, 2*pi) - pi; error_test = mean(abs(phase_error_test)); if error_test < best_error best_error = error_test; best_params = {duties, periods_test}; end end; end; end; end; end工程权衡：全组合扫描计算量大（5⁵=3125 次 RCWA），但
rcwa_1d单次耗时约 20ms（i7 CPU），总耗时可控。生产环境建议改用patternsearch或代理模型加速。
4. RCWA 1D 的三个必调参数与两个典型失效模式
rcwa-1d-02的鲁棒性高度依赖三个核心参数的协同设置。它们不是孤立选项，而是构成数值稳定的三角约束。同时，两类常见失效（伪收敛与相位跳变）有明确诊断路径。
4.1 NG、Nslice、tol 的耦合关系表
| 参数 | 物理意义 | 过小后果 | 过大后果 | 推荐初始值 | 调整策略 |
|---|---|---|---|---|---|
NG | 傅里叶展开阶数 | 高阶衍射被截断，衍射角计算错误；能量不守恒 >5% | 内存占用激增，求解时间平方增长 | 11（覆盖 ±5 级） | 若R_total + T_total < 0.95，先增NG |
Nslice | 单元横向离散数 | 无法分辨精细轮廓，相位响应失真；出现虚假谐振峰 | eps_r矩阵过大，rcwa_1d内部矩阵求逆失败 | 32（亚波长光栅） | 若phi_T随Nslice变化 >0.1 rad，需增大 |
tol | S 矩阵迭代收敛容差 | 未真正收敛，refl/tran值随机波动 | 迭代超时（info.flag = 0），无结果返回 | 1e-8 | 若info.flag == 0，先放宽至1e-6再收紧 |
提示：三者需同步验证。例如，当
Nslice=16时NG=11可收敛，但Nslice=64时需NG=15才能保持精度。无脑增大所有参数会显著拖慢设计迭代速度。
4.2 伪收敛：识别与修复
伪收敛指info.flag == 1（声称收敛），但结果违反物理常识。典型表现：
- 正入射下
R1p和R1m效率差异超过 10⁻³（对称结构应相等）； theta_i = 0时phi_T随period单调增加，但在某点突变 2π（相位跳变）。
诊断代码：
% 检查对称性（适用于周期性或近似对称非周期结构） if max(abs(periods - mean(periods))) < 10e-9 % 近似周期 fprintf('对称性检查: |R1p - R1m| = %.2e\n', abs(R1p - R1m)); if abs(R1p - R1m) > 1e-3 warning('检测到伪收敛：反射对称性破坏，请增大 NG 或 tol'); end end % 检查相位连续性 dphi = diff(phi_T); if any(abs(dphi) > 2.5) % 跳变 > 2.5 rad 即可疑 fprintf('相位跳变位置: %d\n', find(abs(dphi) > 2.5)); warning('存在相位跳变，尝试增大 Nslice 或减小 tol'); end4.3 非周期结构的衍射角修正公式
rcwa_1d输出的refl/tran索引对应的是以平均周期P_avg = mean(periods)为基准的衍射级。但真实衍射角θ_m应满足：
$$ k_0 \sin\theta_m = k_0 \sin\theta_i + \frac{2\pi m}{P_{\text{local}}} $$
其中P_local是局部周期。由于非周期结构无单一P_local，工程上采用等效动量匹配：将m级衍射能量重心位置x_m与目标偏转角关联。rcwa-1d-02不提供此功能，需后处理：
% 计算各衍射级能量重心（单位：米） x_m = zeros(1, length(refl)); P_avg = mean(periods); for m = -(length(refl)-1)/2:(length(refl)-1)/2 idx = m + (length(refl)+1)/2; % 索引转换 x_m(idx) = m * P_avg; % 线性映射到空间域 end % 加权重心位置（用于偏转光栅校准） x_centroid = sum(x_m .* abs(refl).^2) / sum(abs(refl).^2); theta_est = asin(x_centroid * 2*pi / (lambda0 * P_avg)); % 估算偏转角 fprintf('估算偏转角: %.2f deg\n', rad2deg(theta_est));注意：此公式假设光栅长度远大于周期，且
x_m线性分布。对短光栅（<10 单元），需用傅里叶变换直接计算远场分布，而非依赖rcwa_1d的级次输出。
本文还有配套的精品资源，点击获取
