企业官网建设流程全解析

简介：这份资源面向光子晶体光学特性研究与数值计算的学习者，核心是二维平面波扩展法（PWEM）的MATLAB实现，用于分析三类典型晶体的完全禁带。包内共27个m文件，压缩包约78KB，全部为MATLAB脚本，涵盖方形、三角、蜂窝、椭圆、矩形、柱状等多种晶格结构的计算程序，以及若干测试脚本，便于对比不同结构下的禁带结果。程序围绕Bloch定理构建哈密顿量矩阵，通过求解特征值获得频谱分布，进而识别禁带宽度与中心频率，可帮助读者理解平面波基函数选取、周期边界条件处理与结果可视化等关键环节。已有260人学习，适合熟悉MATLAB、希望快速上手光子晶体禁带计算并在此基础上做参数扫描与结构优化的研究者参考。

1. 从一份 .rar 说起：PWEM 到底在算什么

如果你手头正好有一份名为「二维平面波法讨论及光子晶体完全禁带分析程序.rar」的压缩包，里面大概率是一套用 MATLAB 写的平面波展开法（Plane Wave Expansion Method，PWEM）代码，用来算二维光子晶体（Photonic Crystal）的能带结构，进而判断完全禁带（PBG，Photonic Band Gap）落在哪个频率区间。这件事听起来像纯理论，但落到工程上非常具体：你要设计一个光子晶体波导、微腔或者滤波器，第一步就是知道这个周期结构到底在哪个归一化频率范围内对所有偏振、所有传播方向都不让光通过——这就是完全禁带。PWEM 就是干这个的：把麦克斯韦方程在周期介质里展开成傅里叶级数，把介电常数分布写成倒格矢上的平面波叠加，最后把问题变成一个矩阵本征值求解。适合谁？做微纳光学、光子集成、超材料的研究生和工程师，以及需要快速验证结构参数、不想从零推公式的人。它不神秘，但坑不少，下面按我实际跑代码的顺序拆开讲。

2. 平面波法把麦克斯韦方程变成矩阵：原理与最小可跑框架

2.1 为什么是平面波展开，而不是 FDTD 或 FEM

光子晶体能带计算常见三条路：时域有限差分（FDTD）、有限元（FEM）、平面波展开（PWEM）。FDTD 适合宽带响应和透射谱，但算禁带要扫很多频点，还得处理周期边界和 Bloch 相位；FEM 精度高但网格剖分和本征值求解对二维周期结构来说偏重。PWEM 的优势在于：它天然利用周期性，把问题直接写成倒格矢空间里的本征值问题，一次求解就能得到整条能带，判断完全禁带非常直接。代价是介电常数必须是分段常数或能用少量傅里叶分量描述，遇到强色散或金属色散要额外处理。我一般做二维介质柱/空气孔结构时首选 PWEM，快速扫参数；确认候选结构后再用 FDTD 验证透射。

2.2 从波动方程到本征矩阵：E 偏振和 H 偏振分开写

二维光子晶体通常分两种偏振：E 偏振（电场沿柱轴方向，TM 模）和 H 偏振（磁场沿柱轴方向，TE 模）。对 E 偏振，波动方程可以写成标量形式：

$$ -\nabla \cdot \left( \frac{1}{\varepsilon(\mathbf{r})} \nabla E_z \right) = \frac{\omega^2}{c^2} E_z $$

把 $1/\varepsilon(\mathbf{r})$ 和 $E_z$ 都按倒格矢 $\mathbf{G}$ 展开，得到矩阵本征方程：

$$ \sum_{\mathbf{G}'} \left| \mathbf{k}+\mathbf{G} \right| \left| \mathbf{k}+\mathbf{G}' \right| \eta(\mathbf{G}-\mathbf{G}') E_{\mathbf{G}'} = \frac{\omega^2}{c^2} E_{\mathbf{G}} $$

其中 $\eta(\mathbf{G})$ 是 $1/\varepsilon$ 的傅里叶系数。H 偏振的方程形式不同，矩阵元素涉及 $\varepsilon$ 的傅里叶系数和 $(\mathbf{k}+\mathbf{G})\cdot(\mathbf{k}+\mathbf{G}')$ 的叉乘结构。实际写代码时，很多人只算 E 偏振就下结论，这是完全禁带分析里最常见的翻车点——完全禁带要求两种偏振在所有方向上都没有交叠的带隙，只算一个偏振等于没算。

2.3 最小可跑框架：倒格矢生成、介电常数傅里叶系数、本征求解

下面这段 MATLAB 代码是我从零搭 PWEM 时的最小骨架，不依赖任何工具箱，纯矩阵运算。它生成二维正方晶格的倒格矢，计算圆柱形介质柱的介电常数傅里叶系数，然后对给定波矢 k 求解本征频率。

% 二维正方晶格 PWEM 最小框架（E 偏振） clear; clc; a = 1; % 晶格常数，归一化 N = 15; % 每个方向平面波数，总平面波数 (2N+1)^2 r = 0.2*a; % 介质柱半径 eps_bg = 1; % 背景介电常数（空气） eps_rod = 12; % 介质柱介电常数 % 生成倒格矢 G = m*b1 + n*b2，正方晶格 b1=(2pi/a,0), b2=(0,2pi/a) m = -N:N; n = -N:N; [M,Nn] = meshgrid(m,n); Gx = 2*pi/a * M(:); Gy = 2*pi/a * Nn(:); G = [Gx, Gy]; numPW = size(G,1); % 计算 1/eps 的傅里叶系数 eta(G) % 对圆柱结构，eta(G) 有解析式；这里用数值积分做通用写法 eta = zeros(numPW,1); for i = 1:numPW if norm(G(i,:)) < 1e-12 % G=0 时，eta = 面积平均 1/eps eta(i) = (1 - pi*r^2/a^2)/eps_bg + (pi*r^2/a^2)/eps_rod; else % 圆柱的 1/eps 傅里叶系数解析式 f = (1/eps_rod - 1/eps_bg) * 2*pi*r^2/(a^2); eta(i) = f * besselj(1, norm(G(i,:))*r) / (norm(G(i,:))*r); end end % 构造傅里叶系数矩阵 eta(G-G') etaMat = zeros(numPW, numPW); for i = 1:numPW for j = 1:numPW dG = G(i,:) - G(j,:); % 找到 dG 对应的 eta 值，这里用最近邻匹配 [~, idx] = min(sum((G - dG).^2, 2)); etaMat(i,j) = eta(idx); end end % 对给定 k 点求解本征值 k = [0, 0]; % Gamma 点 K = repmat(k, numPW, 1) + G; Knorm = sqrt(sum(K.^2, 2)); A = zeros(numPW, numPW); for i = 1:numPW for j = 1:numPW A(i,j) = Knorm(i) * Knorm(j) * etaMat(i,j); end end eigvals = eig(A); freqs = sqrt(abs(eigvals)) * a / (2*pi); % 归一化频率 a/lambda freqs = sort(freqs(freqs > 1e-6)); disp('前 10 个归一化本征频率：'); disp(freqs(1:min(10, end)));

这段代码的逻辑说明：先按倒格矢定义生成 $(2N+1)^2$ 个平面波，$N$ 控制截断精度；然后算 $1/\varepsilon$ 的傅里叶系数，圆柱结构有解析式，用贝塞尔函数表示，避免数值积分误差；接着构造矩阵 $\eta(\mathbf{G}-\mathbf{G}')$，这里用最近邻匹配是因为倒格矢是离散集合，差值必然落在集合里；最后对给定 k 构造本征矩阵并求本征值，开方得到频率。参数方面，$N=15$ 对折射率对比 3 左右的结构通常够用，$N$ 加到 25 结果变化小于 1% 就可以认为收敛；$r/a$ 和介电常数对比是决定禁带位置和宽度的核心参数。

提示：这段代码只算了 Gamma 点。要画能带图，需要对不可约布里渊区边界上的多个 k 点循环求解，再把同一带序的频率连起来。

3. 完全禁带分析：从能带图到禁带判据的完整流程

3.1 不可约布里渊区路径怎么取，k 点怎么扫

二维正方晶格的不可约布里渊区是 Gamma-X-M-Gamma 三角形。实际扫 k 时，我一般在这三条边上各取 20 到 40 个点，总点数 60 到 120 个。点太少会漏掉带隙交叠的细节，点太多对判断禁带没有额外收益，反而拖慢速度。对六角晶格，路径是 Gamma-M-K-Gamma，形状不同但逻辑一样。关键是要把路径上的 k 点按顺序排列，画图时横轴用累计路径长度，这样能带图才连续。

% 生成正方晶格不可约布里渊区 k 路径 Nk = 30; % 每条边取点数 Gamma = [0, 0]; X = [pi/a, 0]; M = [pi/a, pi/a]; path1 = [linspace(Gamma(1), X(1), Nk)', linspace(Gamma(2), X(2), Nk)']; path2 = [linspace(X(1), M(1), Nk)', linspace(X(2), M(2), Nk)']; path3 = [linspace(M(1), Gamma(1), Nk)', linspace(M(2), Gamma(2), Nk)']; kpath = [path1; path2(2:end,:); path3(2:end,:)];

这段代码生成 k 路径，注意去掉重复端点，避免画图时出现竖直跳变。累计路径长度用cumsum(sqrt(sum(diff(kpath).^2,2)))算，横轴标签在 Gamma、X、M、Gamma 处打点。

3.2 两种偏振分别算，再取交集判断完全禁带

完全禁带的定义是：在同一个归一化频率范围内，E 偏振和 H 偏振的带隙在所有 k 方向上都不存在本征态。操作上就是分别算两条能带图，找出各自的带隙区间，然后取交集。如果交集为空，就没有完全禁带。下面这段代码演示如何从本征频率矩阵中提取带隙。

% 假设 freqsE 和 freqsH 分别是 E 偏振和 H 偏振的本征频率矩阵 % 每行对应一个 k 点，每列对应一个带序 numBands = 8; % 只看前 8 条带 gapE = findGap(freqsE(:,1:numBands)); gapH = findGap(freqsH(:,1:numBands)); gapCommon = [max(gapE(1), gapH(1)), min(gapE(2), gapH(2))]; if gapCommon(1) < gapCommon(2) fprintf('完全禁带范围：%.4f ~ %.4f (a/lambda)\n', gapCommon(1), gapCommon(2)); else fprintf('没有完全禁带\n'); end function gap = findGap(freqs) % 找最大带隙：相邻带之间频率不交叠的区间 numBands = size(freqs, 2); gap = [0, 0]; maxWidth = 0; for b = 1:numBands-1 lower = max(freqs(:,b)); upper = min(freqs(:,b+1)); if lower < upper && (upper - lower) > maxWidth maxWidth = upper - lower; gap = [lower, upper]; end end end

逻辑说明：findGap对每对相邻带，取第 b 带在所有 k 点的最大值和第 b+1 带在所有 k 点的最小值，如果前者小于后者，说明这两条带之间在所有 k 方向上都空着，就是一个带隙。取宽度最大的那个作为该偏振的主带隙。然后两个偏振的带隙取交集，交集非空才是完全禁带。参数方面，numBands取 8 到 12 通常能覆盖前几个禁带；如果禁带出现在更高带序，需要增加平面波数 $N$ 和带序数。

3.3 收敛性验证：平面波数 N 和带序数怎么定

PWEM 的结果对平面波截断数 $N$ 敏感。我一般做三步验证：先用 $N=11$ 粗算，看禁带大概位置；再用 $N=15$ 和 $N=19$ 各算一次，如果禁带上下边界变化小于 1%，认为收敛；如果变化大，继续加到 $N=25$。带序数也要够，至少覆盖禁带上下各两条带。常见错误是 $N$ 太小导致虚假带隙——截断误差把高频带压低，看起来像禁带，实际不存在。验证方法是把 $N$ 翻倍，如果禁带消失或大幅移动，说明之前的结果不可信。

4. 避坑与排查：PWEM 算完全禁带时最容易翻车的 5 个地方

4.1 只算一个偏振就宣布有完全禁带

现象：E 偏振能带图显示明显带隙，直接认为结构有完全禁带。原因：完全禁带要求两种偏振同时存在带隙且频率交叠，E 偏振有带隙不代表 H 偏振也有。解决：必须分别算 E 和 H 偏振，取交集后再下结论。我见过不少论文只给一个偏振的能带图就说完全禁带，审稿人一眼就能看出来。

4.2 傅里叶系数截断导致虚假带隙

现象：$N=11$ 时看到很宽的禁带，$N=15$ 时禁带变窄甚至消失。原因：平面波数不足时，高频分量被截断，矩阵本征值出现系统性偏差，可能把两条带错误地分开。解决：做收敛性测试，$N$ 从 11 加到 25，观察禁带边界变化。如果变化超过 2%，继续加 $N$。对高介电常数对比（比如 12 以上），$N$ 要取到 25 甚至 31。

4.3 倒格矢匹配用错，矩阵元素全乱

现象：本征频率出现负值或复数，能带图完全不可读。原因：构造 $\eta(\mathbf{G}-\mathbf{G}')$ 时，差值 $\mathbf{G}-\mathbf{G}'$ 没有正确映射回倒格矢集合，或者用了错误的索引。解决：倒格矢是离散集合，差值必然等于某个 $\mathbf{G}$，用最近邻匹配时检查距离是否接近零；更稳妥的做法是用整数索引直接算，$\mathbf{G}{m,n} - \mathbf{G}{m',n'} = \mathbf{G}_{m-m',n-n'}$，避免浮点匹配误差。

4.4 归一化频率单位搞混，禁带位置对不上实验

现象：代码算出的禁带在 0.3~0.4，但实验测到的是 0.5 附近。原因：归一化频率 $a/\lambda$ 和实际频率 $\omega$ 的换算搞错，或者晶格常数 $a$ 的单位不一致。解决：统一用归一化频率 $a/\lambda$，其中 $\lambda$ 是真空波长。如果要用实际频率，$f = c \cdot (a/\lambda) / a$，注意 $a$ 用米。我一般全程用归一化频率，最后再换算。

4.5 介质柱形状傅里叶系数用错公式

现象：圆柱结构算出的禁带和文献对不上，偏差超过 10%。原因：圆柱的 $1/\varepsilon$ 傅里叶系数解析式用错，比如贝塞尔函数阶数或参数写错。解决：对圆柱，$\eta(G) = (1/\varepsilon_{rod} - 1/\varepsilon_{bg}) \cdot 2\pi r^2 / a^2 \cdot J_1(Gr)/(Gr)$，其中 $J_1$ 是一阶贝塞尔函数。检查 $G=0$ 时的面积平均公式是否正确。对非圆柱形状，用数值积分算傅里叶系数，别硬套解析式。

5. 进阶技巧：用对称性加速和验证完全禁带

5.1 利用点群对称性缩减矩阵规模

正方晶格有 $C_{4v}$ 对称性，六角晶格有 $C_{6v}$。利用对称性可以把平面波基组按不可约表示分块，矩阵规模缩小到原来的 1/4 到 1/8。对 $N=25$ 的情况，总平面波数 2601，直接求本征值很慢；分块后每块几百阶，速度快很多。实现上，把倒格矢按对称操作分组，构造对称化基函数，再投影到各不可约表示。这一步对快速扫参数很有用，但代码复杂度上升，建议先用全矩阵跑通再优化。

5.2 用带隙图验证：扫参数看禁带宽度变化

完全禁带分析不只是算一个结构，更常见的是扫参数：固定晶格常数，扫介质柱半径 $r/a$ 从 0.1 到 0.4，看禁带宽度怎么变。我一般写一个循环，对每个 $r/a$ 算 E 和 H 偏振的带隙，取交集宽度，画成热力图。这样能快速找到最优 $r/a$。下面是一个扫参数的骨架。

% 扫 r/a 看完全禁带宽度 rOverA = 0.1:0.02:0.4; gapWidth = zeros(size(rOverA)); for idx = 1:length(rOverA) r = rOverA(idx) * a; % 调用 PWEM 计算函数，返回 E 和 H 偏振的带隙 [gapE, gapH] = computeGaps(a, r, eps_bg, eps_rod, N, kpath); gapCommon = [max(gapE(1), gapH(1)), min(gapE(2), gapH(2))]; if gapCommon(1) < gapCommon(2) gapWidth(idx) = gapCommon(2) - gapCommon(1); else gapWidth(idx) = 0; end end plot(rOverA, gapWidth, 'o-'); xlabel('r/a'); ylabel('完全禁带宽度 (a/\lambda)');

这段代码的逻辑：对每个 $r/a$ 调用一次完整的 PWEM 计算，返回两个偏振的带隙，取交集宽度。参数说明：rOverA范围根据结构定，正方晶格介质柱通常 0.1 到 0.4；步长 0.02 够粗扫，精扫可以到 0.005。注意每次调用都要重新算傅里叶系数，因为 $r$ 变了。如果速度慢，可以把傅里叶系数写成 $r$ 的函数预计算。

5.3 和 FDTD 交叉验证：透射谱对不上怎么办

PWEM 算出完全禁带后，我习惯用 FDTD 算一段有限结构的透射谱做交叉验证。如果透射谱在禁带范围内确实掉到 -30 dB 以下，说明 PWEM 结果可信；如果透射谱还有明显透射峰，可能是 PWEM 的禁带判断有误，或者有限结构的尺寸不够、边界反射没处理好。常见对不上的原因：PWEM 用的是无限周期结构，FDTD 用的是有限周期数，禁带边缘的态在有限结构里可能泄漏。我一般把周期数加到 10 以上再对比。

注意：PWEM 和 FDTD 的归一化频率定义要一致，都用电场或磁场分量对应同一个偏振，别一个算 E 偏振一个算 H 偏振。

5.4 我踩过的坑和现在的习惯

最早做光子晶体禁带分析时，我图快只算 E 偏振，结果设计出来的波导在实验里根本不通光，血泪教训。后来养成习惯：任何完全禁带的结论，必须两种偏振都算，必须做 $N$ 收敛，必须用至少一种其他方法交叉验证。MATLAB 代码我会把参数写在脚本开头，跑完自动保存能带图和禁带范围到文本文件，方便回溯。现在看到别人只给一张能带图就说完全禁带，我第一反应就是问：另一个偏振呢？$N$ 取多少？收敛了吗？希望帮到你。

本文还有配套的精品资源，点击获取