简介基于MATLAB实现的平面波扩展法PWE仿真代码可用于模拟具有平面外传播的无限二维光子晶体能带特性。主要面向光子晶体研究、光电器件设计或数值电磁计算的学生与科研人员尤其适合需要快速获取可运行PWE算法框架的中高级MATLAB使用者。压缩包共16个文件核心代码为12个m脚本涵盖参数设置、倒格矢生成、介电常数傅里叶展开、特征值求解及能带绘制等全流程另含2个mat数据文件保存特定结构参数下的工作区结果以及1个md格式说明文档和1个txt许可证文件包体仅32KB轻量易部署。已有222人学习下载。通过调整晶体周期、填充比、背景折射率与平面外波矢等参数可观察不同Kz下的能带结构与带隙变化为理解光子晶体色散关系、设计波导或滤波器件提供直观计算工具。代码结构模块化并附带多种绘制函数便于二次开发。1. 平面外传播的二维光子晶体不是把二维能带平移就能得到如果只在二维布里渊区内扫 k 点TE 和 TM 两个偏振可以各算各的本征问题矩阵规模可以压得很小。可一旦把平面外波矢 Kz 加进来电磁场的纵向分量不再是零原本解耦的偏振会混成一体若仍按二维问题处理能带图在 Kz 附近会出现连续的伪泄漏模。泄漏模恰好是光子晶体光纤、光栅耦合器设计里最关心的部分所以这个场景必须用带平面外传播项的完整 PWE。PhotonicCrystalSim 这套代码就是为这个问题准备的用平面波扩展法求解无限二维光子晶体在任意平面外传播常数下的三维本征问题支持正方、三角两种晶格附带 Kz0 到 6g、R0.30 到 0.50a、nb1.6 的预计算能带数据。对要做光源耦合、泄漏模式分析或带隙扫描的 MATLAB 使用者来说它比从零搭矩阵环境省很多事适合先跑通、再改参数、最后替换成自己的结构。2. PWE 的倒空间矩阵组装epsgg.m 与傅里叶系数2.1 麦克斯韦方程如何退化为代数本征问题对于无磁、非色散介质频率域的麦克斯韦方程可以消去磁场得到电场的旋度方程[ abla imes abla imes \mathbf{E} \frac{\omega^2}{c^2} \varepsilon(\mathbf{r}) \mathbf{E} ]周期结构下介电常数满足 (\varepsilon(\mathbf{r}\mathbf{R}) \varepsilon(\mathbf{r}))。平面波展开的做法是把 (\varepsilon^{-1}(\mathbf{r})) 在倒格矢上展开成傅里叶级数同时把电场写成 Bloch 波形式即每个平面波分量携带一个未知复振幅。代回方程后对每个倒格矢 (\mathbf{G}) 配平指数项旋度运算就变成 ((i\mathbf{k}i\mathbf{G}) \times) 的代数操作。最终所有 (\mathbf{G}) 耦合起来形成一个广义特征值问题矩阵 (M) 里除了波矢与倒格矢的叉乘项其余部分就是介电常数傅里叶系数。PhotonicCrystalSim 中的 epsgg.m 负责构造这些系数它相当于整个求解器的骨架。2.2 介电常数傅里叶系数的圆孔积分对圆孔或圆柱形散射单元介电常数傅里叶系数有解析积分圆形截面上的指数积分可以写成贝塞尔函数 (J_1) 的形式。下面是还原 epsgg.m 计算逻辑的核心片段它和我独立实现过的一段代码思路一致可直接对照function kap epsgg_circle(R, nb, nr, Gx, Gy, cell) % R : 圆截面半径 % nb : 背景折射率 % nr : 孔内/柱内折射率 % Gx : 展平的倒格矢 x 分量N*1 % Gy : 展平的倒格矢 y 分量N*1 % cell: 晶格基包含 ax, ay N numel(Gx); A_cell cell.ax * cell.ay; inv_bg 1/nb^2; inv_in 1/nr^2; kap zeros(N); for i 1:N for j 1:N dGx Gx(i) - Gx(j); dGy Gy(i) - Gy(j); G sqrt(dGx^2 dGy^2); if G 0 int_s pi * R^2; else % 圆形截面的二维傅里叶积分 int_s 2 * pi * R / G * besselj(1, G * R); end kap(i, j) inv_bg * (i j) (inv_in - inv_bg) * int_s / A_cell; end end end代码的核心是把每个倒格矢差 (\mathbf{G}-\mathbf{G}) 对应的面积分算出来。(\mathbf{G}0) 时取圆面积代表平均介电常数部分非零 (\mathbf{G}) 时用贝塞尔函数表示圆对称结构的影响。实际 PhotonicCrystalSim 中的 epsgg.m 还会对三角晶格单元做类似的形状因子修正但圆形积分的数学形式相同。参数说明R控制填充比nb、nr分别对应文件名里的背景折射率和散射体折射率Gx、Gy由倒空间网格生成矩阵维度等于平面波个数的平方。截断越大系数矩阵携带的高频傅里叶成分越完整能带计算越接近真实周期结构。2.3 倒空间截断与矩阵尺寸的伸缩PWE 的精度取决于截断后的平面波数量 (M(2N1)^2)。对于二维结构但带 Kz 的矢量本征问题每个倒格矢点有三分量电场直接组装矩阵阶数是 (3M)。下表给出典型截断下的规模每维截断 N平面波数 M三维全分量矩阵阶数 3M适用场景349147快速试探只看趋势5121363初步对比带底误差约 1%7225675折射率对比 1.6 时的收敛结果93611083发布带隙边界前精算当 Kz 不为零时矩阵计算量比 Kz0 时明显增加因为偏振不再解耦。实际运行中直接对 1000 阶稠密矩阵做eigs单个 Kz 点需要数秒到十几秒。PhotonicCrystalSim 在pwem3DIterKzR.m中采用逐 Kz 点组装、逐点求解的方式把结果写入.mat文件避免让多个 Kz 的矩阵同时驻留内存。这也是我推荐的方式不要对 Kz 连续体构造超矩阵一次性求解分点迭代是标准做法。3. pwem3DIterKzR.m平面外波矢 Kz 的迭代与能带生成3.1 主函数的三层循环结构打开pwem3DIterKzR.m后可以看到它的执行逻辑是三层循环最外层遍历 Kz中层由bz_irr_sqr.m或bz_irr_tri.m给出不可约布里渊区边界点并插值成连续 k 路径最内层通过kvect3D.m把二维 k 向量与当前 Kz 组合交给eigs3D.m求解本征频率。理解这个循环很重要因为每个 Kz 都要重新组装矩阵不能把 Kz0 的结果通过平移简单复用。3.2 运行入口和参数设置下载解压后建议先用addpath把整个目录加入搜索路径再运行主函数。README 中给出的参数定义一般可以改成下面的形式% PhotonicCrystalSim-master 根目录 addpath(fullfile(pwd, Helper Files)); R 0.40; % 圆孔半径单位取晶格常数 a nb 1.6; % 背景折射率 nr 1.0; % 孔内折射率空气孔取 1.0 Nx 7; Ny 7; % 两个方向的倒空间截断 KzRange linspace(0, 6, 7); % 平面外波矢单位 2*pi/a % 运行主函数并保存能带数据 [bandData, kpath, kzAxis] pwem3DIterKzR(R, nb, nr, Nx, Ny, KzRange);pwem3DIterKzR的函数名可以拆成三部分pwem3D表示三维平面波展开IterKz指对 Kz 迭代R则是把半径作为主扫描变量。运行后返回的kpath是归一化 k 路径坐标bandData的每一页对应一个 Kz 值下的能带集合。参数说明R直接决定带隙位置原包保存的预计算数据覆盖 0.30 到 0.50a说明这个区间是带隙容易打开的填充比范围nb1.6对应玻璃或聚合物这类低折射率背景需要较高的截断数才能看清带隙边界KzRange中的 6 本质是 (6 \cdot 2\pi/a)即波矢延伸到布里渊区边界以外用于观察泄漏模。ipanel.m还提供了一个简单滑条界面拖动 R 或 Kz 可以实时刷新能带适合直观感受参数移动对带边的影响。3.3 eigs3D.m 的矢量本征问题eigs3D.m不是把二维矩阵简单复制成三维而是重新构造包含三分量电场耦合的矩阵。Kz 加入后TE 和 TM 的解耦关系消失叉乘项 ((\mathbf{k}\mathbf{G}) \times (\mathbf{k}\mathbf{G}) \times) 会把 (E_x, E_y, E_z) 全部耦合起来。代码先形成 (M\times 3) 的完整倒空间波矢表再组装成块矩阵。若自己写替代品最容易出错的是叉乘项符号和顺序。建议先用 (N_xN_y3) 的小模型把所有矩阵元素打印出来验证 Kz0 时能严格退化为 TE/TM 两个独立本征问题再做大规模扫描。3.4 直接读取保存的 workspace预计算好的数据在SavedWorkspaces目录下文件名直接标明了参数范围Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat就是正方晶格、Kz 从 0 到 6g、半径 0.30 到 0.50a、折射率 1.6 的结果。加载方式如下S load(SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat); whos(-file, SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat)whos会列出文件内的变量名之后用S.变量名取能带值。注意文件名里有逗号和等号load必须使用完整字符串且不能省略扩展名。这些数据不需要先跑主循环就能直接用来后处理和复现论文图比重新扫描快得多。4. 正方/三角晶格的 k 路径与能带后处理bz_irr_sqr.m 到 plotBandStruct.m4.1 不可约布里渊区的高对称点差异正方晶格的二维布里渊区高对称点是 (\Gamma(0,0))、(X(\pi/a,0))、(M(\pi/a,\pi/a))三角晶格则是 (\Gamma)、(M)、(K) 三个点构成的三角形。加入 Kz 后不可约布里渊区从二维平面变成三维棱柱路径通常会变成 (\Gamma-X-M-\Gamma-Z-X) 这种带 Z 轴的组合。bz_irr_sqr.m和bz_irr_tri.m分别生成对应的 k 点连线坐标供主循环和绘图脚本共用。手动输入三角晶格 K 点坐标时最容易出错坐标偏差会让能带图出现假带隙这是移植到六角结构时最大的坑。4.2 用 plotBandStruct 系列出能带图plotBandStruct1.m和plotBandStruct2.m是两种不同风格的绘图入口。前者适合把同一 Kz 下的多条能带叠在一张图上横轴是归一化 k 路径纵轴是归一化频率 (\omega a / 2\pi c)后者适合把多个 Kz 的能带叠加到同一个坐标系观察泄漏模随 Kz 的连续变化。直接运行% 用预计算数据出图 plotBandStruct1(SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat);如果不喜欢默认配色可以自己解析数据。先加载文件再用下面的代码把所有 Kz 的能带画成灰阶线S load(SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat); figure; hold on; KzValues S.KzList; % 实际变量名以 whos 为准 for iz 1:numel(KzValues) plot(S.kpath, S.freq(:,:,iz), Color, [0.6 0.6 0.6]); end xlabel(Wave vector along IBS path); ylabel(Frequency \omega a / 2\pi c);这里假定能带数据按 Kz 存在第三维如果whos看到的变量名是bandData或omega把字段名换掉即可。绘图参数说明kpath是插值后的横坐标不是真实倒空间坐标KzList存储扫描的 Kz 值把线宽设成 0.5 并按 Kz 调整透明度可以更清楚地看出带隙随 Kz 的闭合过程。plotAlpha.m则是把半径 R 当作扫描量快速画出带隙宽度随 R 变化的曲线适合在参数扫描后选择最佳填充比。4.3 用 saved workspace 验证收敛性预计算数据还能用来检查收敛。取同一文件里 R0.40a 的能带对比 Nx5 和 Nx7 的结果导带和价带边频率差通常在 2% 以内。如果带隙边界随截断变化很小说明该模式已收敛如果还在抖则要增加 Nx、Ny同时留意内存占用。常见做法是先用 Nx5 扫一遍 R/Kz 参数网格把带隙最宽的候选点拎出来再用 Nx9 精算这样比一开始就上大矩阵省很多时间。5. bandGaps.m 扫描带隙时的三个经验网格、对称性与内存5.1 只扫布里渊区边界bandGaps.m不是扫描所有 k 点而是只扫描不可约布里渊区的外边界。对二维光子晶体带隙极小值只会出现在高对称点或边界中点因此沿 (\Gamma-X-M-\Gamma) 扫描即可。加入 Kz 后还要补一条沿 Kz 方向的 (\Gamma-Z) 线。实际操作中我在路径转折点附近会密集插值因为带边经常出现在高对称点附近均匀取点容易漏掉最窄带隙。5.2 用 Nx5 粗扫、Nx9 精算下表是不同截断下带隙边界误差和耗时的典型经验值对应 R0.40a、nb1.6 的低折射率对比截断 Nx带隙上边界误差单 Kz 点耗时推荐用途5约 2%0.4 秒粗扫参数域7约 0.8%2 秒带隙趋势研究9约 0.3%8 秒最终精算如果bandGaps.m里没有并行可以在外层 Kz 循环加parforparfor iz 1:numel(KzRange) [bandData(:, :, iz), kpath, kzval] pwem3DIterKzR(R, nb, nr, Nx, Ny, KzRange(iz)); endparfor每个 worker 独立组装和求解内存开销随 worker 数量线性增长所以 Nx9 时不要一次开满 8 个 worker。开 4 个 worker 的速度提升通常最明显因为eigs自身也会占用 BLAS 线程。5.3 剔除伪模式矢量 PWE 偶尔会产生掉到 0 附近的非物理高频模式判断方法是看能带图上是否有独立的孤立点突然落到带隙里而不是跟随正常色散趋势。更稳妥的办法是固定 Nx把 Kz 从 0 增加到 0.01若某个频率点跳动明显大于相邻点那大概率是伪模式应检查倒格矢截断是否对称或去掉该特征值。5.4 合并 Kz 结果与连续色散图把分 Kz 保存的临时数据合并成一张连续色散图mergedBands cat(3, bandData{:}); % bandData 是 cell 时 plotBandStruct2(mergedBands, kpath, KzRange);得到的是横轴为 k 路径、纵轴为频率、颜色或线型表示 Kz 的完整色散图。将 Kz 连续体叠加后能清楚看到带隙在哪个 Kz 处闭合从而定位光子晶体光纤或耦合器的泄漏边界。本文还有配套的精品资源点击获取 
长河君 · 2026/9/13 17:19:11 · 阅读 1.2k 平面波展开法求解光子晶体平面外能带与泄漏模
简介基于MATLAB实现的平面波扩展法PWE仿真代码可用于模拟具有平面外传播的无限二维光子晶体能带特性。主要面向光子晶体研究、光电器件设计或数值电磁计算的学生与科研人员尤其适合需要快速获取可运行PWE算法框架的中高级MATLAB使用者。压缩包共16个文件核心代码为12个m脚本涵盖参数设置、倒格矢生成、介电常数傅里叶展开、特征值求解及能带绘制等全流程另含2个mat数据文件保存特定结构参数下的工作区结果以及1个md格式说明文档和1个txt许可证文件包体仅32KB轻量易部署。已有222人学习下载。通过调整晶体周期、填充比、背景折射率与平面外波矢等参数可观察不同Kz下的能带结构与带隙变化为理解光子晶体色散关系、设计波导或滤波器件提供直观计算工具。代码结构模块化并附带多种绘制函数便于二次开发。1. 平面外传播的二维光子晶体不是把二维能带平移就能得到如果只在二维布里渊区内扫 k 点TE 和 TM 两个偏振可以各算各的本征问题矩阵规模可以压得很小。可一旦把平面外波矢 Kz 加进来电磁场的纵向分量不再是零原本解耦的偏振会混成一体若仍按二维问题处理能带图在 Kz 附近会出现连续的伪泄漏模。泄漏模恰好是光子晶体光纤、光栅耦合器设计里最关心的部分所以这个场景必须用带平面外传播项的完整 PWE。PhotonicCrystalSim 这套代码就是为这个问题准备的用平面波扩展法求解无限二维光子晶体在任意平面外传播常数下的三维本征问题支持正方、三角两种晶格附带 Kz0 到 6g、R0.30 到 0.50a、nb1.6 的预计算能带数据。对要做光源耦合、泄漏模式分析或带隙扫描的 MATLAB 使用者来说它比从零搭矩阵环境省很多事适合先跑通、再改参数、最后替换成自己的结构。2. PWE 的倒空间矩阵组装epsgg.m 与傅里叶系数2.1 麦克斯韦方程如何退化为代数本征问题对于无磁、非色散介质频率域的麦克斯韦方程可以消去磁场得到电场的旋度方程[ abla imes abla imes \mathbf{E} \frac{\omega^2}{c^2} \varepsilon(\mathbf{r}) \mathbf{E} ]周期结构下介电常数满足 (\varepsilon(\mathbf{r}\mathbf{R}) \varepsilon(\mathbf{r}))。平面波展开的做法是把 (\varepsilon^{-1}(\mathbf{r})) 在倒格矢上展开成傅里叶级数同时把电场写成 Bloch 波形式即每个平面波分量携带一个未知复振幅。代回方程后对每个倒格矢 (\mathbf{G}) 配平指数项旋度运算就变成 ((i\mathbf{k}i\mathbf{G}) \times) 的代数操作。最终所有 (\mathbf{G}) 耦合起来形成一个广义特征值问题矩阵 (M) 里除了波矢与倒格矢的叉乘项其余部分就是介电常数傅里叶系数。PhotonicCrystalSim 中的 epsgg.m 负责构造这些系数它相当于整个求解器的骨架。2.2 介电常数傅里叶系数的圆孔积分对圆孔或圆柱形散射单元介电常数傅里叶系数有解析积分圆形截面上的指数积分可以写成贝塞尔函数 (J_1) 的形式。下面是还原 epsgg.m 计算逻辑的核心片段它和我独立实现过的一段代码思路一致可直接对照function kap epsgg_circle(R, nb, nr, Gx, Gy, cell) % R : 圆截面半径 % nb : 背景折射率 % nr : 孔内/柱内折射率 % Gx : 展平的倒格矢 x 分量N*1 % Gy : 展平的倒格矢 y 分量N*1 % cell: 晶格基包含 ax, ay N numel(Gx); A_cell cell.ax * cell.ay; inv_bg 1/nb^2; inv_in 1/nr^2; kap zeros(N); for i 1:N for j 1:N dGx Gx(i) - Gx(j); dGy Gy(i) - Gy(j); G sqrt(dGx^2 dGy^2); if G 0 int_s pi * R^2; else % 圆形截面的二维傅里叶积分 int_s 2 * pi * R / G * besselj(1, G * R); end kap(i, j) inv_bg * (i j) (inv_in - inv_bg) * int_s / A_cell; end end end代码的核心是把每个倒格矢差 (\mathbf{G}-\mathbf{G}) 对应的面积分算出来。(\mathbf{G}0) 时取圆面积代表平均介电常数部分非零 (\mathbf{G}) 时用贝塞尔函数表示圆对称结构的影响。实际 PhotonicCrystalSim 中的 epsgg.m 还会对三角晶格单元做类似的形状因子修正但圆形积分的数学形式相同。参数说明R控制填充比nb、nr分别对应文件名里的背景折射率和散射体折射率Gx、Gy由倒空间网格生成矩阵维度等于平面波个数的平方。截断越大系数矩阵携带的高频傅里叶成分越完整能带计算越接近真实周期结构。2.3 倒空间截断与矩阵尺寸的伸缩PWE 的精度取决于截断后的平面波数量 (M(2N1)^2)。对于二维结构但带 Kz 的矢量本征问题每个倒格矢点有三分量电场直接组装矩阵阶数是 (3M)。下表给出典型截断下的规模每维截断 N平面波数 M三维全分量矩阵阶数 3M适用场景349147快速试探只看趋势5121363初步对比带底误差约 1%7225675折射率对比 1.6 时的收敛结果93611083发布带隙边界前精算当 Kz 不为零时矩阵计算量比 Kz0 时明显增加因为偏振不再解耦。实际运行中直接对 1000 阶稠密矩阵做eigs单个 Kz 点需要数秒到十几秒。PhotonicCrystalSim 在pwem3DIterKzR.m中采用逐 Kz 点组装、逐点求解的方式把结果写入.mat文件避免让多个 Kz 的矩阵同时驻留内存。这也是我推荐的方式不要对 Kz 连续体构造超矩阵一次性求解分点迭代是标准做法。3. pwem3DIterKzR.m平面外波矢 Kz 的迭代与能带生成3.1 主函数的三层循环结构打开pwem3DIterKzR.m后可以看到它的执行逻辑是三层循环最外层遍历 Kz中层由bz_irr_sqr.m或bz_irr_tri.m给出不可约布里渊区边界点并插值成连续 k 路径最内层通过kvect3D.m把二维 k 向量与当前 Kz 组合交给eigs3D.m求解本征频率。理解这个循环很重要因为每个 Kz 都要重新组装矩阵不能把 Kz0 的结果通过平移简单复用。3.2 运行入口和参数设置下载解压后建议先用addpath把整个目录加入搜索路径再运行主函数。README 中给出的参数定义一般可以改成下面的形式% PhotonicCrystalSim-master 根目录 addpath(fullfile(pwd, Helper Files)); R 0.40; % 圆孔半径单位取晶格常数 a nb 1.6; % 背景折射率 nr 1.0; % 孔内折射率空气孔取 1.0 Nx 7; Ny 7; % 两个方向的倒空间截断 KzRange linspace(0, 6, 7); % 平面外波矢单位 2*pi/a % 运行主函数并保存能带数据 [bandData, kpath, kzAxis] pwem3DIterKzR(R, nb, nr, Nx, Ny, KzRange);pwem3DIterKzR的函数名可以拆成三部分pwem3D表示三维平面波展开IterKz指对 Kz 迭代R则是把半径作为主扫描变量。运行后返回的kpath是归一化 k 路径坐标bandData的每一页对应一个 Kz 值下的能带集合。参数说明R直接决定带隙位置原包保存的预计算数据覆盖 0.30 到 0.50a说明这个区间是带隙容易打开的填充比范围nb1.6对应玻璃或聚合物这类低折射率背景需要较高的截断数才能看清带隙边界KzRange中的 6 本质是 (6 \cdot 2\pi/a)即波矢延伸到布里渊区边界以外用于观察泄漏模。ipanel.m还提供了一个简单滑条界面拖动 R 或 Kz 可以实时刷新能带适合直观感受参数移动对带边的影响。3.3 eigs3D.m 的矢量本征问题eigs3D.m不是把二维矩阵简单复制成三维而是重新构造包含三分量电场耦合的矩阵。Kz 加入后TE 和 TM 的解耦关系消失叉乘项 ((\mathbf{k}\mathbf{G}) \times (\mathbf{k}\mathbf{G}) \times) 会把 (E_x, E_y, E_z) 全部耦合起来。代码先形成 (M\times 3) 的完整倒空间波矢表再组装成块矩阵。若自己写替代品最容易出错的是叉乘项符号和顺序。建议先用 (N_xN_y3) 的小模型把所有矩阵元素打印出来验证 Kz0 时能严格退化为 TE/TM 两个独立本征问题再做大规模扫描。3.4 直接读取保存的 workspace预计算好的数据在SavedWorkspaces目录下文件名直接标明了参数范围Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat就是正方晶格、Kz 从 0 到 6g、半径 0.30 到 0.50a、折射率 1.6 的结果。加载方式如下S load(SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat); whos(-file, SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat)whos会列出文件内的变量名之后用S.变量名取能带值。注意文件名里有逗号和等号load必须使用完整字符串且不能省略扩展名。这些数据不需要先跑主循环就能直接用来后处理和复现论文图比重新扫描快得多。4. 正方/三角晶格的 k 路径与能带后处理bz_irr_sqr.m 到 plotBandStruct.m4.1 不可约布里渊区的高对称点差异正方晶格的二维布里渊区高对称点是 (\Gamma(0,0))、(X(\pi/a,0))、(M(\pi/a,\pi/a))三角晶格则是 (\Gamma)、(M)、(K) 三个点构成的三角形。加入 Kz 后不可约布里渊区从二维平面变成三维棱柱路径通常会变成 (\Gamma-X-M-\Gamma-Z-X) 这种带 Z 轴的组合。bz_irr_sqr.m和bz_irr_tri.m分别生成对应的 k 点连线坐标供主循环和绘图脚本共用。手动输入三角晶格 K 点坐标时最容易出错坐标偏差会让能带图出现假带隙这是移植到六角结构时最大的坑。4.2 用 plotBandStruct 系列出能带图plotBandStruct1.m和plotBandStruct2.m是两种不同风格的绘图入口。前者适合把同一 Kz 下的多条能带叠在一张图上横轴是归一化 k 路径纵轴是归一化频率 (\omega a / 2\pi c)后者适合把多个 Kz 的能带叠加到同一个坐标系观察泄漏模随 Kz 的连续变化。直接运行% 用预计算数据出图 plotBandStruct1(SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat);如果不喜欢默认配色可以自己解析数据。先加载文件再用下面的代码把所有 Kz 的能带画成灰阶线S load(SavedWorkspaces/Sqr,Kz0-6g,R0.30-0.50a,nb1.6.mat); figure; hold on; KzValues S.KzList; % 实际变量名以 whos 为准 for iz 1:numel(KzValues) plot(S.kpath, S.freq(:,:,iz), Color, [0.6 0.6 0.6]); end xlabel(Wave vector along IBS path); ylabel(Frequency \omega a / 2\pi c);这里假定能带数据按 Kz 存在第三维如果whos看到的变量名是bandData或omega把字段名换掉即可。绘图参数说明kpath是插值后的横坐标不是真实倒空间坐标KzList存储扫描的 Kz 值把线宽设成 0.5 并按 Kz 调整透明度可以更清楚地看出带隙随 Kz 的闭合过程。plotAlpha.m则是把半径 R 当作扫描量快速画出带隙宽度随 R 变化的曲线适合在参数扫描后选择最佳填充比。4.3 用 saved workspace 验证收敛性预计算数据还能用来检查收敛。取同一文件里 R0.40a 的能带对比 Nx5 和 Nx7 的结果导带和价带边频率差通常在 2% 以内。如果带隙边界随截断变化很小说明该模式已收敛如果还在抖则要增加 Nx、Ny同时留意内存占用。常见做法是先用 Nx5 扫一遍 R/Kz 参数网格把带隙最宽的候选点拎出来再用 Nx9 精算这样比一开始就上大矩阵省很多时间。5. bandGaps.m 扫描带隙时的三个经验网格、对称性与内存5.1 只扫布里渊区边界bandGaps.m不是扫描所有 k 点而是只扫描不可约布里渊区的外边界。对二维光子晶体带隙极小值只会出现在高对称点或边界中点因此沿 (\Gamma-X-M-\Gamma) 扫描即可。加入 Kz 后还要补一条沿 Kz 方向的 (\Gamma-Z) 线。实际操作中我在路径转折点附近会密集插值因为带边经常出现在高对称点附近均匀取点容易漏掉最窄带隙。5.2 用 Nx5 粗扫、Nx9 精算下表是不同截断下带隙边界误差和耗时的典型经验值对应 R0.40a、nb1.6 的低折射率对比截断 Nx带隙上边界误差单 Kz 点耗时推荐用途5约 2%0.4 秒粗扫参数域7约 0.8%2 秒带隙趋势研究9约 0.3%8 秒最终精算如果bandGaps.m里没有并行可以在外层 Kz 循环加parforparfor iz 1:numel(KzRange) [bandData(:, :, iz), kpath, kzval] pwem3DIterKzR(R, nb, nr, Nx, Ny, KzRange(iz)); endparfor每个 worker 独立组装和求解内存开销随 worker 数量线性增长所以 Nx9 时不要一次开满 8 个 worker。开 4 个 worker 的速度提升通常最明显因为eigs自身也会占用 BLAS 线程。5.3 剔除伪模式矢量 PWE 偶尔会产生掉到 0 附近的非物理高频模式判断方法是看能带图上是否有独立的孤立点突然落到带隙里而不是跟随正常色散趋势。更稳妥的办法是固定 Nx把 Kz 从 0 增加到 0.01若某个频率点跳动明显大于相邻点那大概率是伪模式应检查倒格矢截断是否对称或去掉该特征值。5.4 合并 Kz 结果与连续色散图把分 Kz 保存的临时数据合并成一张连续色散图mergedBands cat(3, bandData{:}); % bandData 是 cell 时 plotBandStruct2(mergedBands, kpath, KzRange);得到的是横轴为 k 路径、纵轴为频率、颜色或线型表示 Kz 的完整色散图。将 Kz 连续体叠加后能清楚看到带隙在哪个 Kz 处闭合从而定位光子晶体光纤或耦合器的泄漏边界。本文还有配套的精品资源点击获取