一维光子晶体传输矩阵法Matlab仿真：从原理到反射谱计算

简介：面向光学仿真研究者、高校相关专业学生及光学器件设计工程师，这套基于MATLAB的一维光子晶体传输矩阵仿真代码包含12个.m文件，覆盖能带结构、透射率与反射率计算等常见场景，适合用来理解周期性介质层中的光传播规律。压缩包整体仅9KB，脚本从折射率与厚度的周期性排布入手，包含结构参数定义、矩阵连乘与边界处理的完整范例，可作为搭建个人光子晶体仿真平台的基础模板，亦可直接修改参数开展批量分析和教学演示。目前已有818人学习下载，文件虽小但流程清晰，在光学滤波器、光开关等器件的初步设计与验证中有实际参考价值。通过阅读运行这些脚本，读者既能巩固传输矩阵方法的物理原理，又能掌握MATLAB数值仿真的实用技巧，快速获得透射光谱等关键结果，为后续深入研究与性能调优打下基础。

1. 一维光子晶体与传输矩阵方法：这个rar包能复现什么

如果你手里有一个命名为“20181022.rar_Photonic crystal_matlab_一维光子晶体_传输矩阵_光子晶体”的压缩包，打开后大概率会看到几个.m文件和一两个fig图。这个标题的核心，是用Matlab实现一维光子晶体的传输矩阵法，计算多层膜结构的反射谱和透射谱。这类代码在光学工程课设、光子晶体入门研究里出现频率极高，因为它用不到几十行矩阵运算就能复现出带隙、缺陷模这些一维光子晶体的关键现象，比直接上FDTD或COMSOL便宜得多。

它适合谁？两类人：一是正在做课程设计或毕业设计，需要快速看懂一段Matlab光子晶体代码并改参数出图；二是想验证自己推导的传输矩阵公式是否正确，拿一组已知结构当对照。接下来我会从物理模型讲到可运行的Matlab脚本，再把调参和常见翻车点拆开说，最后给你几个验证程序正确性的方法。

2. 传输矩阵为什么适合算一维光子晶体：从麦克斯韦方程到Matlab建模

2.1 传输矩阵的物理来源：每层材料只是一个2×2矩阵

一维光子晶体是折射率沿一个方向周期性变化的结构，比如常见的(高折射率层/低折射率层)N堆叠。光在其中传播时，每一层里的电场可以写成前向波和反向波的叠加。在层与层的界面上，电场和磁场的切向分量必须连续，这个边界条件就可以用一个2×2矩阵把某一侧的电场磁场状态映射到另一侧。

对第i层介质，特征矩阵的形式是：

M_i = [cos(δ_i), i·sin(δ_i)/η_i] [i·η_i·sin(δ_i), cos(δ_i) ]

其中δ_i = 2π·n_i·d_i·cos(θ_i)/λ，是光在这一层内传播的相位厚度；η_i是导纳，TE偏振下η_i = n_i·cos(θ_i)，TM偏振下η_i = n_i/cos(θ_i)。这个矩阵把层的一侧状态向量（电场、磁场）变换到另一侧，整叠多层膜的总传输矩阵就是把每一层的矩阵按光传播方向依次连乘。

为什么要强调“按光传播方向依次连乘”？因为矩阵乘法不满足交换律，先乘高折射率层还是先乘低折射率层，最终反射谱可能完全不同。后面避坑章节我会专门讲这个顺序问题。

2.2 三种常见仿真方法对照与选型：为什么代码包里几乎全是传输矩阵

我见过不少新手拿到这种代码包后第一反应是“能不能用FDTD重写一遍”，我的建议是：如果只算一维周期性结构的反射/透射谱，别换。一维光子晶体用传输矩阵法，本质上就是闭式解析解，毫秒级扫完数百个波长点；FDTD要设网格、边界条件、光源和监视器，一次仿真就得几秒甚至几分钟，而且数值色散会吃掉部分带隙细节。

对于标题里这种“一维光子晶体 + 传输矩阵”的组合，传统做法就是选传输矩阵法。它不只能算反射率，还能顺带拿到透射率、吸收率和场分布，而且把折射率改成复数就可以模拟有损耗或金属层。

2.3 从代码包视角：先确认rar里这几样东西再动手

我不确定你手里那个“20181022.rar”内部文件具体怎么命名，但按这类Matlab光子晶体代码的通用组织形式，压缩包里一般会有：一个主脚本负责定义参数和调函数，一个或多个函数文件负责单层矩阵计算、整体反射率计算，可能还有.fig图形文件保存着某个默认结构的反射谱截图。

你拿到包之后，不要急着运行，先做三件事：第一，找到主脚本，看它开头几十行是不是有n1、n2、d1、d2、N、lambda_range这些参数定义，没定义清楚的话程序跑出来也是黑匣子；第二，检查有没有独立的layerMatrix或transfer_matrix函数，没有的话主脚本里应该也有内联矩阵构建代码；第三，确认输出变量到底给的是反射率R、透射率T还是两者的对数坐标，这直接影响你后续怎么解读曲线。

3. 用Matlab跑通一维光子晶体：最小可运行的传输矩阵代码

3.1 最小脚本：正入射下的(高/低折射率)N多层膜反射谱

我直接给一份可以复制运行的最小脚本。它算的是空气 |（Si/SiO2）重复8个周期| 空气的结构，中心波长1550 nm，扫描范围800～2000 nm，正入射TE波（正入射时TE和TM没有区别）。

% demo_1d_photonic_crystal.m % 一维光子晶体传输矩阵法，正入射反射/透射谱 % 结构: Air | (H L)^N | Air, H=Si, L=SiO2 clear; clc; % ---------- 结构参数 ---------- nH = 3.48; % Si 近红外折射率 nL = 1.46; % SiO2 折射率 lam0 = 1550e-9; % 设计中心波长 1550 nm dH = lam0 / (4 * nH); % H层厚度: 四分之一波长光学厚度 dL = lam0 / (4 * nL); % L层厚度 N = 8; % 周期数 theta0 = 0; % 正入射 % ---------- 入射/出射介质 ---------- n0 = 1.0; % 空气入射 ns = 1.0; % 空气出射 % ---------- 波长扫描 ---------- lamList = linspace(800e-9, 2000e-9, 1201); R = zeros(size(lamList)); T = zeros(size(lamList)); for idx = 1:length(lamList) lam = lamList(idx); % 总传输矩阵初始化为单位阵 M_total = eye(2); % 按光入射方向逐层连乘: 先H后L, 重复N次 for k = 1:N M_H = layerMatrix(nH, dH, lam, theta0, 0); M_L = layerMatrix(nL, dL, lam, theta0, 0); M_total = M_total * M_H * M_L; end % 反射系数与透射系数 eta0 = n0 * cosd(theta0); etas = ns * cosd(theta0); F = M_total(1,1) + M_total(1,2) * etas; G = M_total(2,1) + M_total(2,2) * etas; r = (eta0 * F - G) / (eta0 * F + G); t = 2 * eta0 / (eta0 * F + G); R(idx) = abs(r)^2; T(idx) = (etas / eta0) * abs(t)^2; end % ---------- 绘图 ---------- figure('Name','1D Photonic Crystal'); plot(lamList*1e9, R, 'LineWidth', 1.5); hold on; plot(lamList*1e9, T, '--', 'LineWidth', 1.5); xlabel('Wavelength (nm)'); ylabel('Reflectance / Transmittance'); legend('R','T'); grid on;

3.2 子函数：单层特征矩阵与导纳计算

上面的脚本依赖一个子函数layerMatrix。如果你用的Matlab版本是R2016b之后，可以直接把下面这个函数写在主脚本同一个文件的末尾；如果是老版本，单独存成layerMatrix.m。

function M = layerMatrix(n, d, lam, theta0, pol) % layerMatrix 计算单层介质的特征矩阵 % 输入: % n - 该层折射率 % d - 该层几何厚度 (m) % lam - 真空波长 (m) % theta0 - 入射角 (度), 针对入射介质 % pol - 0=TE, 1=TM % 输出: % M - 2x2 特征矩阵 % 由Snell定律求层内折射角: n0*sin(theta0) = n*sin(theta) % 这里入射介质默认为空气 n0=1 sinTheta = sind(theta0) / n; theta = asind(sinTheta); % 层内角度 % 相位厚度: delta = 2*pi*n*d*cos(theta)/lambda delta = 2 * pi * n * d * cosd(theta) / lam; % 导纳: TE 用 n*cos(theta), TM 用 n/cos(theta) if pol == 0 eta = n * cosd(theta); else eta = n / cosd(theta); end % 特征矩阵 M = [cos(delta), 1i * sin(delta) / eta; 1i * eta * sin(delta), cos(delta) ]; end

3.3 代码逻辑拆解：矩阵连乘、导纳比和复数单位

这个脚本的核心就是那行M_total = M_total * M_H * M_L，它模拟的是光从空气进入第一层Si，穿过SiO2，再进入下一周期Si……直到最后从结构出射进入空气。连乘的结果是一个把出射侧电场/磁场映射到入射侧的2×2矩阵。反射系数公式r = (eta0*F - G) / (eta0*F + G)里的F和G分别是总矩阵第一行和第二行与基底导纳组合后的量，逻辑上对应入射侧与出射侧场匹配方程的解。

有几个参数值得注意。第一，eta0和etas在正入射时都等于介质折射率，所以正入射时T = abs(t)^2，但斜入射或出射介质与入射介质不同时，必须保留etas/eta0这个导纳比因子，否则透射率会违背能量守恒。第二，代码里用1i而不是i表示虚数单位，这是为了避免程序里某个循环变量恰好叫i时把虚数单位覆盖掉，这种问题在Matlab旧版本里很常见。第三，linspace(800e-9, 2000e-9, 1201)把波长范围分成了1200个间隔，对Si/SiO2这种折射率对比度很大的结构已经够密；但如果你后续要算高Q值缺陷模，这个采样密度远远不够，需要把缺陷模附近的波长范围单独加密。

这段代码跑完之后，你应该能看到一个明显的反射高平台，那就是一维光子晶体的带隙。我第一次跑通时以为代码写错了，因为带隙里反射率几乎是1，透射率几乎为0，曲线陡得不像话——后来确认这就是多层膜堆叠的正常表现。

4. 让反射谱按你的要求走：周期数、入射角与厚度参数的5个必调项

4.1 周期数N与折射率对比度：带隙深度由它们决定

周期数N直接决定带隙内的反射率能逼近1到什么程度。N=2时带隙内反射率可能只有80%，N=8基本能到99%以上，N=20后增加已经不明显，反而会让矩阵连乘次数增加、浮点误差累积。我做参数扫描时一般先把N设成1、2、4、8、16各跑一遍，看带隙中心反射率的变化趋势，再选一个曲线已经收敛的N作为正式值。

折射率对比度Δn = nH - nL决定了带隙的宽度。一维光子晶体带隙宽度有个近似估算公式：

Δω/ω0 ≈ (4/π)·arcsin(|nH - nL| / (nH + nL))

Si/SiO2的对比度约2.0，算出来相对带宽超过50%，所以800到2000 nm的谱里你会看到一段很宽的反射平台。如果换成折射率对比度很小的材料对，比如1.5和1.6，带隙会窄很多，这时候波长采样点数就要加大，否则带隙边界会走出锯齿状。

4.2 入射角与偏振：TE/TM分裂是怎么出现的

把theta0从0改成30度、45度、60度，你会发现两件事：带隙整体往短波方向移动；TE和TM的反射谱开始不一致，TM带隙会比TE窄。原因是斜入射时每个特征矩阵里的θ_i变了，TE和TM的导纳表达式不同，等效光学厚度也随之变化。

这里有个实操技巧：在Matlab里做角度扫描时，不要只跑一个theta0画一张图，而是把角度循环和波长循环嵌套起来，输出一个以波长为横轴、角度为纵轴的颜色图，带隙区域会在图上形成一条亮带。这种图在论文里很有说服力，也能直观看到TE/TM分裂。

但注意，正入射时TE和TM完全简并，没必要把pol设成1去跑两遍。只有在theta0 > 0时才需要区分偏振。

4.3 光学厚度比例：从四分之一波长堆到啁啾结构

标题里这个代码包很大概率是基于四分之一波长堆叠设计的，也就是每层的光学厚度等于λ0/4。这种结构在中心波长处反射率最高，带隙关于中心波长对称。但如果你把dH和dL改成不同的光学厚度比例，比如让dH = λ0/4nH、dL = λ0/2nL，带隙位置和宽度都会变化，甚至会出现新的透射峰。

我一般会用“光学厚度比”dHnH : dLnL作为控制变量来扫参数。固定λ0 = 1550 nm，把dL的光学厚度从0.2λ0扫到0.5λ0，可以看到带隙边界怎么移动。这个思路也是把代码从“复现别人的图”升级成“自己设计结构”的第一步。

4.4 扫描范围与采样点数：别让谱线丢细节

很多人拿到代码包后直接跑一遍出图就完事了，但如果你要看的缺陷模或带隙边界特别窄，采样点太少会把峰吃掉。经验做法是先算一次粗扫描，比如400个点，观察带隙大致位置；然后对带隙边缘或者透射峰附近做1000～2000点的密集扫描。不需要全波段都用高密度采样，那只会增加计算时间。

还有一个易忽略的细节：linspace的两个端点要覆盖你关心的特征位置。比如设计中心波长1550 nm，扫描范围如果只到1500 nm，那中心反射峰都没进去，还说“帮同学跑Matlab光子晶体课设作业”，很容易被一眼看出没理解结构含义。

4.5 Matlab版本兼容与数值精度：老代码、新环境的小问题

2018年的代码包放到现在的Matlab版本（比如2026b）里跑，最常见的问题不是语法，而是编码和函数兼容。旧脚本如果有中文注释，在老版本里用GBK编码保存，新版默认UTF-8读出来就是乱码，但不影响运行。真正影响运行的，是某些老代码用了i当循环变量，导致虚数单位被覆盖，计算结果全变成NaN——这种问题排查起来很费劲，你看到的反射谱是一堆警告和空图。

我的习惯是拿到代码后先把所有虚数单位统一改成1i，再把循环变量改成ii或kk，一行行看有没有同名冲突。数值精度上，矩阵连乘200层以内双精度足够，不需要开符号计算；但如果你用format long观察带隙边缘的反射率，会发现它不完全等于0或1，这是正常的数值现象，不是bug。

5. 一维光子晶体仿真避坑：R+T不等于1与另外5个翻车现场

5.1 透射率大于1或R+T明显失衡

现象：算出来的T在某些波长超过1，或者R+T和1差了很多。

原因：透射系数公式漏了导纳比值etas/eta0。正入射且入射、出射介质都是空气时，etas/eta0 = 1，所以问题被掩盖了；一旦你把基底换成玻璃或硅片，漏掉比值后T就会算错。

解决：统一用T = (etas / eta0) * abs(t)^2计算透射率，不要偷懒写成T = abs(t)^2。另外检查材料折射率是否加了虚部，有损耗时R+T本来就小于1，差值就是吸收。

5.2 反射谱乱跳，带隙位置完全对不上

现象：曲线形态像反射谱，但带隙出现在完全错误的波长位置，甚至和理论值差一倍。

原因：厚度d定义错了。四分之一波长堆叠的正确写法是d = lam0 / (4*n)，其中n是该层材料折射率。有人写成了d = lam0 / 4，等于把光学厚度算成了四分之一真空波长，带隙自然整体偏移。

解决：先单独检查dH、dL数值。以1550 nm为例，Si层厚度应该是111 nm左右，SiO2层约265 nm。这个量级不对，后面就不用看了。

5.3 斜入射时出现NaN或曲线断裂

现象：theta0设成60度或更大时，反射谱出现空值，或者报错说asin的参数超出[-1,1]。

原因：Snell定律n0·sin(theta0) = n·sin(theta)里，如果材料折射率小于入射介质的n0·sin(theta0)，层内折射角就不存在实数解，对应全反射和倏逝波。asind(sinTheta)里的sinTheta大于1就产生复数角度，Matlab返回NaN。

解决：在layerMatrix函数里判断abs(sinTheta) > 1时，改用复数形式的θ和η，或者限制扫描角度在临界角以下。对于空气入射的Si/SiO2结构，临界角主要受SiO2折射率限制，theta0超过约43度就会进入这个区域，和实验里观察到的现象一致。

5.4 矩阵连乘顺序反了，反射谱和透射谱完全失真

现象：代码看起来没问题，跑出来的谱也有带隙模样，但带隙位置或深浅明显不合理，尤其多层膜两边介质不对称时更明显。

原因：传输矩阵连乘顺序没按光入射方向来。有的教程从基底侧往入射侧推，公式里的M顺序和这里完全相反。M_total = M_H * M_L和M_total = M_L * M_H结果不一样，因为矩阵乘法不交换。

解决：在自己代码里固定一套约定。我推荐的是“入射侧先乘，出射侧后乘”，也就是本篇文章第3章代码里的M_total = M_total * M_H * M_L。如果你参照的教材用相反约定，那反射系数公式里的F、G定义也要跟着换，否则就会混搭出错。

5.5 矩阵乘法写成点乘，维度报错或结果诡异

现象：代码提示矩阵维度不匹配，或者不报错但反射谱上全是毛刺。

原因：特征矩阵是2×2矩阵，连乘必须用*。有些人写的时候顺手用了.*，结果变成逐元素相乘，2×2矩阵没报错但物理意义完全错了。

解决：在Matlab编辑器里高亮搜索.*，凡是用在矩阵变量之间的基本都要改回*。初始化M_total时用eye(2)而不是ones(2)，不然第一层矩阵乘下去会把两个通道混起来。

5.6 老代码里的i被循环变量覆盖，谱线全是NaN

现象：整个反射谱输出都是NaN或只有零星几个点，而且不报错。

原因：旧版代码习惯用i做循环变量，比如for i = 1:length(lamList)。循环结束后i变成数字，后面再用i当虚数单位的地方全部失效。

解决：把所有虚数单位换成1i，循环变量改成ii、kk这种。这是Matlab老代码移植到新版本最典型的隐性坑，排查时先看有没有红波浪线，如果没有再查变量覆盖。

6. 验证你的Matlab光子晶体程序：解析对照与能量守恒检查

代码跑通只是第一步，怎么确定结果可信？我给你三个验证手段，按工作量从小到大排。

第一个是单层膜解析对照。把N设成0（即只有一层膜），单层介质的反射率有解析公式，正入射时R = ((n0*ns - n^2)^2 + ((n0*ns - n^2)^2? )不用记那么复杂，你只需对比程序算出的反射率极值位置和理论值。单层介质在光学厚度为λ0/4的奇数倍时反射率达到极值，这是检验特征矩阵公式符号是否正确最灵敏的方式。

第二个是能量守恒检查。对无损耗结构，每个波长下R+T必须等于1。在程序里加一段检查代码：

energyError = abs(R + T - 1); fprintf('Max energy error: %.3e\n', max(energyError));

如果误差大于1e-10，说明导纳比或透射系数公式有问题。我第一次写这个程序时就卡在这里，透射谱在带隙边缘高出反射谱一截，最后发现就是etas/eta0写反了。

第三个是带隙边界近似公式对照。用第4章那个Δω/ω0公式估算带隙宽度，再看程序里R=50%的两个交点间距。对Si/SiO2这种高对比度结构，近似公式和传输矩阵结果会有几个百分点的偏差，但量级应该吻合。如果偏差超过30%，你就要怀疑周期数或厚度定义出了问题。

我做这类光子晶体仿真走下来的经验是：不要一上来就去调材料或者加缺陷层，先把最简单的四分之一波长堆叠跑通，验证R+T=1，再用单层膜解析式校准一次，然后才敢拿这个程序去算新结构。否则你后面加缺陷层、做啁啾结构、引入金属层时，任何一步出错都很难定位。这个习惯帮我省了很多时间，希望帮到你。

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