简介本资源基于MATLAB实现一维光子晶体传输矩阵法的仿真程序集适合光学工程、凝聚态物理方向的科研人员及相关领域高年级学生使用。程序通过定义各介质层的折射率与厚度计算光子晶体的透射率、反射率与吸收率等关键光学参数可帮助使用者理解光子禁带与局域模态的形成机理并为光学滤波器、光开关等器件设计提供数值参考。压缩包内共12个文件全部为.m脚本涵盖多个版本的主程序、测试函数与功能子程序便于对照不同建模思路包体大小仅9KB清晰简洁适合快速上手和二次修改。目前已有818人学习下载。通过分析这些MATLAB源码读者不仅能掌握传输矩阵方法的核心推导与矩阵运算技巧还能获得一套可复用的光子晶体仿真框架用于参数扫描、结构优化及结果可视化从而在周期性介质结构设计与光学特性调控方面开展进一步探索。1. 一维光子晶体MATLAB仿真传输矩阵代码包到底帮你省掉哪些事一维光子晶体最反直觉的地方在于只要把两种折射率不同的介质按四分之一波长光学厚度交替叠起来光就会在某个波长范围被整片反射哪怕每一层都比头发丝细得多。更麻烦的是带隙边缘对层厚误差极其敏感厚度差 1%透射峰可能移动几十纳米手工估算完全靠不住。这套资源就是一组用 MATLAB 实现传输矩阵法TMM的脚本核心是 new4.m、v20181022.m 等 .m 文件通过矩阵连乘计算一维光子晶体的透射率、反射率和吸收率。它能直接回答“层数加到多少带隙才够深”“中心波长怎么移动”这类器件设计问题。适合做光子晶体课程设计、光学滤波器和光开关预研的人也适合刚接触 TMM 想找一份能复现参考代码的 MATLAB 使用者。2. 传输矩阵方法把多层膜拆成矩阵连乘MATLAB 里的实现顺序2.1 为什么一维光子晶体用 TMM 而不是 FDTD光子晶体仿真常见的方法有时域有限差分FDTD、严格耦合波分析RCWA和传输矩阵法。FDTD 能处理复杂几何结构但在一维周期性膜堆上反而显得笨重需要在空间和时间上离散网格边界条件处理不好还会产生伪反射每扫一个波长都要重新跑完整时域演化。一维光子晶体的每一层都是均匀介质光在其中传播可以写成解析形式的平面波用 TMM 直接得到频域响应计算量小两个数量级。TMM 的核心思想是把整堆薄膜看成一系列 2×2 矩阵的连乘。每一层对应一个传输矩阵描述光波在该层内的相位累积和界面处的电场连续性把所有层矩阵按光传播方向依次相乘得到整堆膜的总体矩阵。最终透射系数、反射系数都从这个总体矩阵的对应元素里读出来。这个过程和 MATLAB 的矩阵运算天然匹配这也是这套资源选择 TMM 的根本原因。用 TMM 还有一个工程上的好处它天然适合参数扫描。改变某一层厚度、某一层折射率只需要重算对应层的矩阵不必重建整个电磁模型。对于光学滤波器设计里“改一个参数看一次曲线”的迭代过程TMM 几乎是零成本。RCWA 虽然精度更高但对一维垂直入射场景属于杀鸡用牛刀。所以这套脚本的选型是合理的不是图省事。2.2 单层传输矩阵的推导与 MATLAB 表达先看单层介质。假设光从折射率 n1 的区域垂直进入厚度为 d、折射率为 n2 的介质层再进入 n3 区域。在 TMM 中每一层内部只发生两件事界面处的折射/反射以及层内传播的相位变化。垂直入射时TE 和 TM 没有区别光学导纳 η 就等于该层折射率 n。更一般地斜入射时 TE 模式取 η n·cosθTM 模式取 η n/cosθ其中 θ 是该层内的折射角。每一层可以用一个 2×2 矩阵表示最常用的形式是function M layer_matrix(n, d, lambda, theta) % 单层介质传输矩阵 % n: 介质折射率 % d: 介质物理厚度单位与 lambda 保持一致 % lambda: 真空中波长 % theta: 层内折射角默认 0 表示垂直入射 delta 2 * pi * n * d * cos(theta) / lambda; M [cos(delta), -1i * sin(delta) / n; -1i * n * sin(delta), cos(delta)]; end这段代码里的 delta 就是光波在该层内积累的相位厚度。如果 d 和 lambda 都用纳米delta 就是无量纲的弧度。cos(theta) 在垂直入射时等于 1矩阵退回最简形式。矩阵的第二行第一列带虚数单位 i这是为了同时表达界面反射的相位跳变。注意这里把导纳 n 直接写进矩阵所以斜入射时要根据偏振替换成对应的导纳值否则 TE/TM 会混在一起。我一般会把这个函数单独存成 layer_matrix.m然后在主脚本里反复调用。这样改材料参数时只改主脚本不碰矩阵逻辑。这份资源里的 new4.m、new5.m 大概率就是这种结构主脚本负责定义结构参数和波长扫描TMM 函数独立存在。如果你的版本里函数和主脚本混在一个文件里复制出来单独保存即可MATLAB 会自动搜索同目录下的 .m 文件。2.3 全套代码里最值得读的循环结构TMM 主循环的核心不是矩阵求逆而是“连续右乘”。假设一维光子晶体由 N 对高低折射率层组成入射光从左边进来那么最先接触光的层必须放在连乘顺序的最左边也就是第一个右乘到单位矩阵上。% 光子晶体结构N 对高低折射率层 N 8; M eye(2); % 2x2 单位矩阵 for k 1:N M M * layer_matrix(nH, dH, lambda, 0); % 先高折射率层 M M * layer_matrix(nL, dL, lambda, 0); % 再低折射率层 end % 从入射介质到出射介质的透射系数 t 1 / M(1,1); T abs(t)^2;这里的 M eye(2) 是连乘的起点表示还没有经历任何层时电场保持不变。每经过一层就把该层的传输矩阵右乘到当前总矩阵 M 上。为什么要右乘而不是左乘因为光传播方向的顺序和矩阵乘法顺序相反先经过的层放在更右边。如果写反带隙位置和深度都会错后面第 4 章会专门讲这个坑。T abs(t)^2 成立的前提是入射介质和出射介质折射率相同。如果出射介质是衬底透射率还要乘上折射率比修正项代码会变成T abs(t)^2 * real(n_out) / real(n_in);n_in 是入射介质折射率n_out 是出射介质折射率。这一点在阅读 new4.m 时一定要先确认如果脚本只写 abs(t)^2说明作者默认两侧是空气如果写了折射率比说明考虑了衬底。不要强行套用。我在复现这种老脚本时第一件事就是搜索代码里有没有出现 n_sub、n_out 这类变量以此判断边界条件。3. 跑通 new4.m参数定义、循环计算与透射率曲线3.1 文件归属与运行入口判断这套资源里包含 new.m、new1.m、new2.m、new3.m、new4.m、new5.m、new6.m、v20181007.m、v20181009.m、v20181011.m、v20181022.m、test1.m。从命名规律看v 开头的是带日期的版本new 系列是功能迭代版本test1.m 很可能是调试脚本。运行入口我建议优先看 new4.m或者 v20181022.m因为标题里带 20181022且 new4.m 正好对应这个编号。如果你直接运行某个脚本报错“函数或变量未定义”不要急着改代码先检查 MATLAB 当前路径是否包含所有 .m 文件。老脚本经常拆成多个函数文件漏一个就会断链。常见做法是把整个解压目录添加到 MATLAB 路径或者用 cd 切到该目录下。我一般会先跑 test1.m因为它通常是最简单的验证脚本能快速检查环境是否通。test1.m 如果画出一条曲线说明基础函数没问题再跑 new4.m 看完整结果。3.2 介质层参数怎么设折射率、厚度、波长范围跑 TMM 之前要确定四个参数高低折射率层的折射率、每层物理厚度、周期数、波长扫描范围。最常见的材料组合是 TiO2 和 SiO2折射率分别取 2.35 和 1.46。也有用 Si 和 SiO2 的折射率差更大带隙更宽。厚度通常按四分之一波长波片设计也就是 d λ0 / (4n)其中 λ0 是目标中心波长。% 结构参数设定单位统一用纳米 nH 2.35; % 高折射率层 nL 1.46; % 低折射率层 lambda0 600; % 目标中心波长 600 nm dH lambda0 / (4 * nH); % 高折射率层厚度 dL lambda0 / (4 * nL); % 低折射率层厚度 N 8; % 周期数8 对高低折射率层 lambda linspace(400, 800, 801); % 波长扫描范围这里的 dH 约等于 63.8 nmdL 约等于 102.7 nm。注意 dH 和 dL 不一样因为四分之一波长光学厚度相同物理厚度要除以各自的折射率。很多新手把两层设成一样厚结果带隙位置偏到完全找不到。如果脚本里能看到类似 dH lambda0 / (4 * nH) 的赋值说明作者遵守了四分之一波长条件如果直接写 d 100那就要小心结果是否可靠。波长扫描范围可以取目标中心波长附近比如 400 nm 到 800 nm。步长用 801 个点分辨率约 0.5 nm对大多数滤波特性分析够用。如果你只需要看带隙深度可以加密到 2001 个点但计算量大一倍而曲线形状几乎不变。我在工程中习惯先用粗步长找趋势再在带隙边缘加密这比一开始就用超高分辨率效率高得多。3.3 从矩阵连乘到透射率输出的代码骨架为了让你快速跑起来我把这套资源的核心逻辑整理成一个可运行的简化骨架假设出射介质和入射介质都是空气也就是悬浮膜模型% 简化版一维光子晶体 TMM 透射率计算 clear; clc; % 材料参数 nH 2.35; nL 1.46; lambda0 600; dH lambda0 / (4 * nH); dL lambda0 / (4 * nL); N 8; % 波长扫描范围 lambda linspace(400, 800, 801); T zeros(size(lambda)); for idx 1:length(lambda) M eye(2); for k 1:N M M * layer_matrix(nH, dH, lambda(idx), 0); M M * layer_matrix(nL, dL, lambda(idx), 0); end t 1 / M(1,1); T(idx) abs(t)^2; end % 画图 plot(lambda, T, b-, LineWidth, 1.5); xlabel(Wavelength (nm)); ylabel(Transmittance); title(1D Photonic Crystal TMM Result); grid on;这段代码的运行逻辑很直接外层按波长循环内层把 8 对高低折射率层依次右乘到总矩阵上。每个波长点都重新算一遍矩阵连乘所以第 9 行到第 14 行是计算最密集的部分。theta 传 0 表示垂直入射没有必要为每个层单独计算角度。如果要在衬底上计算例如入射侧是空气、出射侧是玻璃需要把最后一行改成 T(idx) abs(t)^2 * 1.5 / 1并且把层数为奇数还是偶数考虑进去。很多实际器件都长在衬底上所以请务必分清脚本里的 n_out 变量。v20181022.m 这类文件可能就包含多个版本有的带衬底修正有的不带运行时注意区分。3.4 如何判断这条曲线算对了跑出曲线之后第一件事不是调参数而是验证结果是否符合物理直觉。你可以在 600 nm 附近看到明显的反射带也就是透射率凹陷区域。带隙中心应该大致落在 lambda0 附近带隙宽度取决于折射率对比度。如果曲线整体接近直线大概率是矩阵顺序写反或者折射率单位错了。如果透射率超过 1则可能是边界修正项没有除干净。如果带隙边缘出现尖锐的震荡那往往是正常现象这是有限层数带来的法布里-珀罗干涉条纹。我判断脚本是否跑通的步骤通常是先看曲线在目标波长附近有没有凹陷再看背景透射率是否接近 1最后看带隙内的最小透射率是否低于 1%对于 8 对层来说带隙深度一般能做到小于 1%。如果这些都满足说明代码逻辑和参数设置基本正确。如果曲线完全向左或向右平移再检查厚度公式是不是 d lambda0 / (4n)而不是 d lambda0 / 4。4. 一维光子晶体仿真避坑传输矩阵最容易翻车的五个细节4.1 现象透射率超过 1 或出现负值我第一次跑通 TMM 时透射率曲线在高透区直接冲到 1.3一度怀疑是光子晶体自带增益。后来才意识到这是振幅透射率和功率透射率没分清。abs(t)^2 给出的是电场振幅平方只有当入射介质和出射介质折射率相同时它才等于功率透射率。两侧折射率不同光通过界面时波阻抗会变化功率流密度不等于振幅平方。原因就是漏掉了折射率比修正项。常规做法是T abs(t)^2 * real(n_out) / real(n_in)。如果出射介质是衬底n_out 取衬底折射率如果出射介质是空气n_out 取 1。对应到这份资源里打开脚本先看有没有乘这个修正项。没有的话就手动加上或者把两侧介质设成一致。4.2 现象带隙中心波长整体偏移几十纳米最近复现一个滤波器时我要的中心波长是 532 nm算出来带隙中心却落在 498 nm。查了很久才发现厚度公式里把 n 当成了波长系数写成了 d lambda0 / 4而不是 d lambda0 / (4n)。高折射率层如果按错误公式算厚度会偏大太多光学厚度直接翻倍带隙自然往长波方向跑。解决方法是统一公式和单位。先明确 lambda0 的单位然后 dH lambda0 / (4nH)dL lambda0 / (4nL)。如果忘了哪一层该除以哪个折射率记住一个原则光学厚度 n*d 必须等于 lambda0/4。这个条件满足时带隙中心才会落在设计波长上。4.3 现象矩阵连乘顺序写反带隙消失有的脚本算出来透射率曲线平平无奇完全没有带隙但参数看起来都对。问题往往出在矩阵连乘顺序上。TMM 中光传播方向极其重要先接触光的层必须先右乘。写成 M M * layer_matrix(nH, ...) 是标准做法如果写成 M layer_matrix(nH, ...) * M就相当于从反方向入射整个结构的光学响应完全变了。我在 new3.m 这类早期版本里就见过这种问题。解决方法是加一行注释明确指出“先进入光线的层放在右边”。验证方法很简单如果高低折射率层互换位置透射率曲线应当发生变化如果曲线完全不变说明矩阵可能是对角矩阵这不可能一定是哪里写错了。4.4 现象斜入射时 TE/TM 结果混为一谈当你把代码从垂直入射扩展到斜入射时层矩阵里的导纳项会出问题。垂直入射时 TE 和 TM 没有区别直接取 n 就行。斜入射时TE 模式的等效导纳是 n*cosθTM 模式是 n/cosθ。如果还用原来的 n算出来的透射率随角度变化完全不对。常见做法是在 layer_matrix 函数里加一个偏振参数s 表示 TEp 表示 TM。theta 是各层内的折射角需要通过 Snell 定律从入射角求出。如果不打算做斜入射那正入射的脚本不用动但如果要做光学开关或角度调谐必须把偏振项拆开。这份资源大概率只覆盖正入射所以你扩展时要注意。4.5 现象波长扫描步长太粗带隙边缘细节被吞我见过有人在 400 nm 到 800 nm 范围内只取 100 个点步长 4 nm算出来的带隙边缘看起来像一条平滑的直线完全看不出振铃结构。带隙边缘的陡峭程度和干涉振荡是特征指标步长太粗会把最大透射率点漏掉导致后续分析得出错误结论。解决方法是分阶段扫描先 linspace(400, 800, 401) 粗扫找到带隙大致范围再对该范围加密比如 linspace(550, 650, 2001)。这样既保证精度又不浪费算力。对于普通 8 层结构2001 个波长点一次扫描也就几秒钟完全值得。5. 用这套代码做滤波器和参数扫描从单点计算到带隙优化5.1 透射率随周期数变化的扫描脚本拿到 new4.m 之后你很快会想回答一个问题周期数 N 到底取多少合适。N 太小带隙不深N 太大制备困难且计算变慢。我一般会把主脚本包成一个函数输入波长数组和周期数返回透射率数组然后循环扫描 N。% 固定波长范围扫描不同周期数 lambda linspace(400, 800, 801); Ns [4 6 8 12 16]; figure; for i 1:length(Ns) T tmm_sweep(lambda, Ns(i), nH, dH, nL, dL); plot(lambda, T, LineWidth, 1.2); hold on; end legend(N4,N6,N8,N12,N16); xlabel(Wavelength (nm)); ylabel(Transmittance); grid on;这里 tmm_sweep 是把第 3 章骨架封装成的函数内部实现就是双层循环。从结果会看到N4 时带隙底部还有不少透射N8 时带隙深度已经低于 1%N16 以后边缘变陡但中心深度改善有限。对于带通滤波器设计8 到 12 对是比较经济的区间。再多层数不仅计算慢实际镀膜时的厚度误差也会让理论曲线失真。这个扫描脚本最大的用途是帮你建立直觉带隙深度的增长不是线性的头几层效果最明显后面边际收益递减。如果你设计的光开关需要 20 dB 以上的抑制比那至少需要 10 对以上如果只是观察带隙现象4 对就足够。5.2 中心波长位置怎么调厚度缩放与折射率对比度调整中心波长有两种办法一是整体缩放所有层的物理厚度二是改变目标波长公式。工程上最常用的是厚度缩放因为折射率一般是材料属性不容易连续变化。% 把中心波长从 600 nm 调整到 850 nm scale_factor 850 / 600; dH2 dH * scale_factor; dL2 dL * scale_factor;这里的 scale_factor 是目标波长和当前波长的比值。只要保持光学厚度为 lambda0/4带隙位置就成比例移动。注意折射率对比度不变时带隙的相对宽度大致不变但绝对宽度会随中心波长增加而变宽。比如 600 nm 处带隙宽度如果是 80 nm搬到 850 nm 后可能变成 113 nm 左右。如果你想增大带隙宽度只能提高折射率对比度。比如把高折射率层从 TiO2 换成 SinH 从 2.35 提到 3.5带隙会明显变宽。但 Si 在可见波段吸收较大只适合近红外。这就是材料和结构参数的权衡目标波长、带隙宽度、吸收损耗互相制约。用 TMM 可以快速扫出曲线再结合材料实际吸收选择折中方案。5.3 把结果导出成设计图表仿真参数调好后下一步往往是要把曲线放进报告或论文里。MATLAB 的 figure 窗口直接截图分辨率不够我习惯先用 writematrix 把数据存成 CSV再用 Origin 或 Python 重新绘图。也可以直接在 MATLAB 里用 print 导出高分辨率图片。% 保存透射率数据 T_out [lambda(:), T(:)]; writematrix(T_out, photonic_crystal_transmittance.csv); % 导出 300 dpi 的 PNG 图片 set(gcf, Color, white); print(gcf, -dpng, -r300, bandgap_result.png);writematrix 是 R2019a 之后的函数如果你的 MATLAB 版本较老可以用 csvwrite 或 xlswrite。保存 CSV 的另一个好处是方便对比不同周期数、不同厚度下的曲线不必每次都重新打开脚本。我一般会把关键参数也写进文件名比如 bandgap_N8_dH63_experiment1.csv这样后期整理数据时不会混乱。导出之后还要注意一张图不要塞太多曲线。5.1 节里的 N 扫描图最多四条曲线再多就分不清了。工程报告里更清晰的做法是画成 2×2 子图分别展示不同周期数的透射率、反射率和相位变化。相位数据可以从 M(1,1) 和 M(1,2) 的复角提取不过在滤波分析中透射率通常已经够用。6. 验证你的结果能量守恒与网格独立性检查6.1 能量守恒T R A 1拿到一组透射率曲线后我会立刻算一下反射率。TMM 算出总矩阵后反射系数 r -M(2,1)/M(1,1)前提是入射介质和出射介质相同。如果不考虑吸收透射率和反射率之和应当等于 1。实际材料总会有点吸收但纯介质膜在可见波段的吸收很小TR 低于 0.99 就说明代码里有问题。r -M(2,1) / M(1,1); R abs(r)^2; disp([TR , num2str(sum(T(1:10:end) R(1:10:end)) / length(1:10:end))]);这个检查能一次性暴露两个问题一是透射率修正项有没有漏二是矩阵元素取错了。我在跑 v20181022.m 时就用这个方法发现某个版本在长波端 TR 只有 0.94原因是波长数组里有几个点算出来 M(1,1) 特别小数值稳定性出问题。解决方法是把波长范围稍微缩小或者改用双精度格式。6.2 网格独立性加密波长点曲线应稳定第二个验证是网格独立性检查。如果你的透射率曲线在波长步长从 801 点加密到 2001 点后带隙边缘位置变化小于 0.5 nm说明结果已经收敛。如果加密后曲线明显变化说明之前的步长太粗带隙边缘细节不可信。我一般会跑两次第一次用 linspace(400, 800, 501)第二次用 linspace(400, 800, 2001)然后把两条曲线画在同一张图里对比。如果两条线完全重合就用粗的那版做后续优化速度快如果不重合说明 501 点的结果不可信必须用更细的网格。从那以后我每次跑完 TMM 都会强制走一遍能量守恒和网格加密验证再用这套流程检查新写的代码。虽然多花两分钟但能避免把错曲线写进报告里的尴尬。希望帮到你。本文还有配套的精品资源点击获取