简介本资源是一套基于MATLAB实现的传输矩阵法TMM计算工具面向光学、电磁学方向的研究生、科研人员及工程技术人员用于快速建模与分析波导结构的反射率、透射率等关键光学特性。压缩包为RAR格式共2个文件均为MATLAB源代码.m文件体积仅1KB轻量高效涵盖单层传输矩阵构建、全局矩阵组装及有效折射率搜索等核心功能模块如logF.m负责场传播计算Search_neff.m用于求解波导模式有效折射率。目前已有311人学习下载体现了其在教学演示与基础科研中的实用价值。用户可直接运行脚本输入波导几何参数与材料属性获得定量光学响应结果并结合MATLAB绘图功能直观呈现频谱特性是理解TMM物理内涵与掌握数值实现路径的精简型实践范例。1. 波导光学仿真不靠商业软件用这套MATLAB TMM代码3分钟跑出TE/TM模式透射谱你手头有一段多层介质波导结构——比如SiO₂/Si/SiO₂脊形波导或者AlGaAs/GaAs量子阱堆叠想快速知道它在1550 nm附近有没有导模、有效折射率是多少、透射率随波长怎么变别急着打开COMSOL或Lumerical——这套名为TMM_WG.rar的MATLAB代码包就是专为这类问题设计的轻量级求解器。它不依赖PDE工具箱不调用任何外部编译器纯矩阵运算单文件启动输入厚度、折射率、入射角三组参数就能输出neff、R、T、相位响应全量结果。适合光学器件初筛、课程设计建模、论文附录验证也适合作为自研光子仿真框架的底层TMM模块。如果你熟悉复数运算和特征方程求根但不想重写400行边界条件拼接逻辑这套代码就是现成的“可读、可改、可嵌入”的TMM最小可行实现。2. TMM核心原理与MATLAB实现逻辑为什么用2×2传输矩阵而不是直接解麦克斯韦方程2.1 为什么TMM是波导分析的“黄金折中”对均匀各向同性介质中的平面波电磁场满足亥姆霍兹方程。若结构沿传播方向z分段均匀每段内场可表示为正向反向行波叠加$$ E(z) A e^{i\beta z} B e^{-i\beta z} $$其中$\beta k_0 n_{\text{eff}}$为传播常数。TMM的核心洞察在于相邻界面处的电场与磁场连续构成线性约束而每层内部的场演化是确定性相位旋转。因此无需在整个空间离散化求解偏微分方程只需将每层抽象为一个2×2复数矩阵描述该层对入射/反射波振幅的线性变换关系$$ \begin{bmatrix} E^ \ E^- \end{bmatrix}_{\text{out}}\mathbf{M}j \begin{bmatrix} E^ \ E^- \end{bmatrix}{\text{in}} $$提示这里的$E^$和$E^-$不是电场分量而是沿z和-z方向传播的复振幅系数。TMM本质是把物理问题转化为线性代数问题——这正是MATLAB最擅长的领域。2.2 TE与TM模式的矩阵形式差异及MATLAB编码要点TEs偏振和TMp偏振模式因边界条件不同其单层传输矩阵结构不同。关键区别在于TM模式需引入阻抗修正因子即$\eta_j \frac{\cos\theta_j}{n_j}$而TE模式直接使用$n_j$。MATLAB中必须严格区分% TE模式电场平行于界面s-polarization M_TE [exp(1i*beta_j*d_j), 0; ... 0, exp(-1i*beta_j*d_j)]; % TM模式磁场平行于界面p-polarization eta_j cos(theta_j) / n_j; % 阻抗归一化因子 M_TM [exp(1i*beta_j*d_j), 0; ... 0, exp(-1i*beta_j*d_j)] * diag([1, eta_j/eta_{j1}]);2.2.1logF.m中的特征方程构建逻辑logF.m并非直接求解本征值而是构造波导横向谐振条件的标量函数$$ F(n_{\text{eff}}) \log\left| \det\left( \mathbf{M}{\text{total}} - \mathbf{I} \right) \right| $$当$F(n{\text{eff}})$出现极小值接近负无穷即对应满足闭合路径相位条件的导模有效折射率。MATLAB中采用fminbnd配合自适应采样% 在合理区间内搜索neff例如1.4~3.5 neff_range [1.4, 3.5]; neff_guess fminbnd((neff) logF(neff, layers, lambda, theta), ... neff_range(1), neff_range(2), ... optimset(TolX, 1e-6, MaxIter, 200));注意logF函数返回的是对数模值避免数值下溢fminbnd比fzero更鲁棒因特征方程在neff虚部非零时可能无实根而log|det|始终有定义。2.3 全局矩阵组装的MATLAB向量化实现传统TMM代码常采用for循环逐层相乘易出错且无法利用MATLAB的BLAS加速。本包采用预分配累积乘法% 初始化全局矩阵为单位阵 M_total eye(2); % 按物理顺序从入射侧到出射侧遍历各层 for j 1:length(layers) M_j compute_layer_matrix(layers(j), lambda, theta, mode); % 返回2x2矩阵 M_total M_j * M_total; % 注意乘法顺序先入射层后透射层 end2.3.1 层参数结构体的设计与校验layers必须是结构体数组每个元素含n复折射率、d厚度、typedielectric/metal字段layers(1).n 1.45; layers(1).d 1e-6; layers(1).type dielectric; % SiO2 layers(2).n 3.48 1e-3i; layers(2).d 220e-9; layers(2).type dielectric; % Si layers(3).n 1.45; layers(3).d Inf; layers(3).type substrate;提示Inf厚度表示半无限衬底其传输矩阵简化为单位阵复折射率n nr 1i*ni自动支持吸收损耗计算ni0时透射率自然衰减。3. 实战从解压到绘制TE₀模透射谱的完整MATLAB工作流3.1 环境准备与代码加载解压TMM_WG.rar后得到TMM_WG/目录内含logF.m特征方程目标函数Search_neff.mneff搜索主函数TMM_WG.m或.mat主调用脚本或工作区变量确保MATLAB路径包含该目录addpath(D:\TMM_WG); % 替换为你的实际路径 which Search_neff % 应返回 D:\TMM_WG\Search_neff.m3.2 定义三层SiO₂/Si/SiO₂波导结构创建waveguide_config.m% 波长扫描范围单位米 lambda_vec linspace(1.5e-6, 1.6e-6, 201); % 结构参数SI单位 layers struct(); layers(1).n 1.444; layers(1).d 1e-6; layers(1).type cladding; layers(2).n 3.477 1e-5i; layers(2).d 220e-9; layers(2).type core; layers(3).n 1.444; layers(3).d Inf; layers(3).type substrate; % 入射条件 theta 0; % 正入射 mode TE; % 或 TM % 预分配结果数组 neff_vec zeros(size(lambda_vec)); T_vec zeros(size(lambda_vec)); R_vec zeros(size(lambda_vec));3.3 批量计算neff与透射率for idx 1:length(lambda_vec) lambda lambda_vec(idx); % 搜索基模有效折射率TE₀ [neff, ~] Search_neff(layers, lambda, theta, mode, max_mode, 1); neff_vec(idx) neff; % 计算该波长下的透射率T和反射率R [T, R, ~] TMM_calculate(layers, lambda, theta, mode, neff); T_vec(idx) T; R_vec(idx) R; end3.3.1TMM_calculate.m的关键参数说明该函数需传入已知neff由Search_neff返回内部执行计算每层传播常数beta_j 2*pi/lambda * sqrt(n_j^2 - neff^2)构造每层2×2矩阵TE/TM分支累积得全局矩阵M_total由M_total(1,1)和M_total(2,1)推导入射端反射系数$r M_{21}/M_{11}$透射系数$t 1/M_{11}$最终$R |r|^2$$T \frac{\text{Re}(n_{\text{out}} \cos\theta_{\text{out}})}{\text{Re}(n_{\text{in}} \cos\theta_{\text{in}})} |t|^2$能量守恒修正注意TMM_calculate中n_out取最后一层折射率n_in取第一层当theta0时$\cos\theta1$公式简化为$T |t|^2$但代码仍保留通用形式以支持斜入射。3.4 可视化与物理验证figure(Position, [100, 100, 1200, 400]); subplot(1,3,1); plot(lambda_vec*1e6, neff_vec, b-o, MarkerSize, 4); xlabel(\lambda (\mum)); ylabel(n_{eff}); title(Effective Index); grid on; subplot(1,3,2); plot(lambda_vec*1e6, T_vec, r-s, MarkerSize, 4); xlabel(\lambda (\mum)); ylabel(T); title(Transmittance); ylim([0, 1.05]); grid on; subplot(1,3,3); plot(lambda_vec*1e6, R_vec, k-d, MarkerSize, 4); xlabel(\lambda (\mum)); ylabel(R); title(Reflectance); ylim([0, 1.05]); grid on;3.4.1 验证能量守恒R T ≈ 1添加断言检查数值精度energy_error max(abs(R_vec T_vec - 1)); fprintf(Max energy conservation error: %.2e\n, energy_error); % 合理阈值1e-12双精度极限至1e-8浮点累积误差 if energy_error 1e-7 warning(Energy not conserved! Check layer impedance matching.); end4. 进阶技巧处理金属包层波导与多模竞争的稳定收敛策略4.1 金属层带来的数值病态性及MATLAB应对方案当结构含Au、Ag等金属层n nr i*ni,ni 1时beta_j变为强衰减复数导致exp(i*beta_j*d_j)在指数项实部很大时产生数值溢出。标准做法是分离相位与衰减% 原危险写法可能导致Inf或NaN phase_term exp(1i * real(beta_j) * d_j - imag(beta_j) * d_j); % 安全写法用real/exp单独处理 real_part cos(real(beta_j)*d_j) * exp(-imag(beta_j)*d_j); imag_part sin(real(beta_j)*d_j) * exp(-imag(beta_j)*d_j); phase_term real_part 1i*imag_part;TMM_WG包中compute_layer_matrix函数已内置此保护但需确认其是否启用——检查源码中是否有if imag(n_j) 1e-2分支。4.2 多模搜索的初始化陷阱与Search_neff.m参数调优Search_neff默认只找一个根但波导常支持多个导模TE₀, TE₁, TM₀...。要获取全部需分段扫描去重% 将neff区间划分为5段每段独立搜索 neff_segments linspace(1.4, 3.5, 6); all_neff []; for k 1:length(neff_segments)-1 range_k [neff_segments(k), neff_segments(k1)]; neff_k fminbnd((neff) logF(neff, layers, lambda, theta), ... range_k(1), range_k(2), opts); % 去重与已有结果差值0.01则跳过 if isempty(all_neff) || all(abs(neff_k - all_neff) 0.01) all_neff [all_neff, neff_k]; end end4.2.1Search_neff.m关键可调参数表参数名默认值作用调整建议max_mode1最大搜索模式数设为3可强制找前3个nefftol_neff1e-4neff收敛容差高精度需求设为1e-6init_step0.1初始步长金属波导建议降至0.01use_logFtrue是否用logdet4.3 快速验证用已知解析解交叉检验对于单层介质波导空气/玻璃/空气TE₀模有效折射率近似为$$ n_{\text{eff}} \approx n_{\text{film}} - \frac{0.63\lambda}{4\pi d} $$取n_film1.5,d1um,lambda1.55um理论值≈1.499。运行代码layers struct(n,{1,1.5,1}, d,{Inf,1e-6,Inf}, type,{cladding,core,cladding}); [neff, ~] Search_neff(layers, 1.55e-6, 0, TE); fprintf(Computed neff %.6f, Theory ≈ 1.499\n, neff); % 输出应为 1.4989xx误差0.001若偏差超过0.01检查layers.d单位是否为米非微米、lambda是否为米、mode是否匹配TE/TM。本文还有配套的精品资源点击获取