简介资源为基于FDTD时域有限差分方法的MATLAB仿真项目专注于多层膜结构的反射与透射特性分析可计算不同波长下的反射谱与透射谱适用于光学涂层、滤光片、太阳能电池等器件的设计与优化包括高反射膜、增透膜等常见设计场景也适合光电子、电磁场专业的课程设计和科研入门。压缩包内有4个MATLAB脚本整体仅5KB轻量易读其中包含一维FDTD核心算法实现以及网格参数、边界条件、光源和结果处理等辅助代码通过修改脚本中的膜层厚度、折射率、入射角等参数即可快速获得新的反射透射结果。目前已有517人学习下载。借助该项目不仅可以直观理解光在多层膜中的干涉机制还可将一维核心代码作为基础进一步扩展至二维/三维FDTD仿真或膜系自动优化兼具教学和科研参考价值。1. 为什么用 1D FDTD 算多层膜的反射谱做光学涂层设计的人都有这种体会传输矩阵法算一维多层膜的反射谱只要几毫秒但当膜系里出现吸收介质、非线性材料或者想观察光在膜层内部的瞬态传播过程时它就不够用了。这时 FDTD时域有限差分反而更直观它直接把麦克斯韦方程组在空间和时间上离散化光怎么进去、怎么在界面上来回反射、最后怎么透出去每一步场分布都看得见。这个fdtd-matlab.rar包里恰好就是一套用 MATLAB 写的 1D FDTD 实现四个脚本分别覆盖主求解器、参数设置和数据后处理。选择 1D FDTD 而不是 2D/3D是因为多层膜问题的场只在膜层法线方向变化横向无限延伸降维到一维后计算量骤减普通笔记本几秒钟就能跑完一个宽谱脉冲的反射谱。这套代码的价值不在于算法本身有多新而在于它把「光源 → 膜系 → 反射/透射监测」这条链路完整打通了。新手可以直接改折射率和厚度看反射谱变化老手则可以把它当成骨架往上加色散模型、斜入射和 PML 优化。下面从网格、时间步、边界条件开始逐个拆解这四个脚本的实际用途。2. 网格、时间步与边界条件1D FDTD 在多层膜中的参数骨架2.1 空间步长与 CFL 条件FDTD 的精度受限于空间步长 $\Delta x$。在多层膜仿真中膜层材料折射率通常为 1.03.5如 SiO₂ 约 1.46、TiO₂ 约 2.35、Si 约 3.5介质中的最小波长是 $\lambda_0 / n_{max}$。要获得可接受的反射谱精度通常要求 $\Delta x \le \lambda_{min} / 20$更保守的做法是取 $\lambda_{min} / 30$。时间步长 $\Delta t$ 受 CFL 条件约束在一维情况下为c · Δt / Δx ≤ 1为了保证数值稳定性工程上一般取 $S c \cdot \Delta t / \Delta x 0.5$这个值在大多数 1D FDTD 实现中足够安全还能兼顾计算效率。% 基础参数定义常见做法 lambda_ref 0.55e-6; % 参考波长单位米 n_max 3.5; % 材料最大折射率 dx lambda_ref / (n_max * 30); % 空间步长保证最密材料也有30个点 dt dx / (3e8 * 0.5); % CFL0.5安全余量 total_T 4000; % 总时间步数需覆盖脉冲在膜系中多次反射完成这段代码建立了两组核心参数空间步长决定入射光在每层膜内被离散成多少个网格点时间步长决定场推进的节奏。n_max要取整个膜系里折射率最高的材料否则高频成分在该材料内会出现明显的数值色散。total_T的取值取决于膜系总光程——膜层越厚光在里面往返所需步数越多这个值取小了反射谱会出现截断伪影。2.2 界面处的介电常数处理多层膜的每个界面上折射率发生跳变。FDTD 电场更新公式中会出现介电常数 $\varepsilon$如果直接在每个网格点上用该点所属材料的 $\varepsilon$界面处会产生非物理的微小反射。常见的做法是取相邻两个网格点介电常数的平均值也就是在界面位置做平滑处理。% 折射率分布构建界面处做算术平均 Nx 2000; % 总网格数 n ones(1, Nx) * 1.0; % 默认空气 for k 1:numel(layer_thickness) idx_start round(sum(layer_thickness(1:k-1)) / dx) 1; idx_end round(sum(layer_thickness(1:k)) / dx); n(idx_start:idx_end) layer_n(k); end n_avg (n(1:end-1) n(2:end)) / 2; % 界面平均折射率 eps n_avg.^2; % 相对介电常数这里把离散的折射率数组转换成相邻网格的平均值再用平方得到介电常数。这样处理之后电场更新方程中每个网格位置上的 $\varepsilon$ 都对应实际材料的局部平均值避免了阶梯近似造成的虚假反射。layer_thickness数组需要用户在脚本里按实际膜系填写单位与dx保持一致。2.3 PML 吸收边界与光源位置计算域两端必须设置吸收边界否则反射波碰到边界会反弹回来污染结果。一维 FDTD 中常用 PML完美匹配层在边界处放置若干层吸收介质。层数太少、吸收系数太小边界反射会叠加到反射谱上形成振荡纹波。PML 层数适用场景边界反射水平8快速验证、粗看谱形约 1% 量级16常规反射谱计算约 0.1% 量级32高精度反射率标定可低于 0.01%光源的位置也应该规划清楚。把光源放在 PML 区域外、第一个反射监测点和膜系之间入射电场沿着正向传播先到达膜系反射波再回到监测点。通常设置两个监测点一个在光源与膜系之间记录入射场加反射场另一个在膜系后端记录透射场。监测点距离膜系至少 10 个网格避免近场效应对结果的干扰。3. fdtd_1d_3.m 求解器拆解高斯脉冲激励与 FDTD 反射谱提取3.1 主循环中的电磁场交替推进fdtd_1d_3.m是这个包里编号为 3 的 1D FDTD 主求解器文件名里的_3说明它经过了至少两轮迭代修改。核心循环在每一时间步内依次更新磁场、电场再处理 PML 吸收层。以下是一个极简但完整的 1D FDTD 核心循环骨架% 场数组初始化 Ez zeros(1, Nx); % 电场 Hy zeros(1, Nx); % 磁场 % PML 导电率分布多项式渐变 sigma_max 0.75; % 最大电导率按经验取值 sigma zeros(1, Nx); pml_width 32; sigma(1:pml_width) sigma_max * ((pml_width:-1:1)/pml_width).^2; sigma(end-pml_width1:end) sigma_max * ((1:pml_width)/pml_width).^2; % 时间步进主循环 for t 1:total_T % 磁场更新Hy(n) Hy(n) - (dt/(mu*dx)) * (Ez(n1) - Ez(n)) Hy(1:end-1) Hy(1:end-1) - (dt / (mu0 * dx)) * (Ez(2:end) - Ez(1:end-1)); % 电场更新含 PML 损耗项 Ez(2:end-1) Ez(2:end-1) - (dt ./ eps(2:end-1) / dx) .* (Hy(2:end-1) - Hy(1:end-2)); % PML 内电场衰减 Ez(1:pml_width) Ez(1:pml_width) ./ (1 sigma(1:pml_width) * dt / eps(1:pml_width)); Ez(end-pml_width1:end) Ez(end-pml_width1:end) ./ (1 sigma(end-pml_width1:end) * dt / eps(end-pml_width1:end)); % 高斯脉冲源注入硬源或柔软源 t0 60; tau 15; source exp(-((t - t0)/tau)^2); Ez(src_pos) source; end磁场与电场在时间上差半步更新这是 FDTD 的 Yee 网格格式的基本特点保证了二阶精度。PML 层的电导率从边界向内渐变吸收系数过大会产生数值反射过小则吸收不充分sigma_max 0.75是一维问题中比较稳的起点值。高斯脉冲的tau直接决定频谱宽度tau越小脉冲越窄、频谱越宽能覆盖的波长范围越大但需要更小的空间步长来保证高频成分不失真。3.2 反射谱的频域提取反射谱不能在时域直接看需要把监测点记录的时域电场做傅里叶变换再用反射场频谱除以入射场频谱得到反射系数。这里有三种做法先跑一次无膜系纯空气的仿真记录参考波 $E_{ref}(t)$再跑有膜系的仿真记录 $E_{total}(t)$相减得到反射波 $E_{refl}(t)$然后分别做 FFT 求比值。在膜系两侧各设一个监测点分别得到反射场和透射场直接用透射频谱除以参考频谱得到透射率。用连续波扫描一个波长一个波长地算稳态反射率精度最高但耗时。方案 1 是实际最常用的因为一次性宽谱激励就能覆盖整个目标波段代码如下% 监测点时域记录假设已经同时记录参考场 ref_Ez 和总场 total_Ez refl_Ez total_Ez - ref_Ez; % 反射场 总场 - 入射场 % FFT 提取频域响应 Fs 1 / dt; % 采样率 N_FFT 2^nextpow2(length(refl_Ez)); freq (0:N_FFT-1) / N_FFT * Fs; lambda 3e8 ./ freq(2:end); % 频率转波长 R_spectrum abs(fft(refl_Ez, N_FFT)) ./ abs(fft(ref_Ez, N_FFT)); R_spectrum R_spectrum(2:end); % 去掉直流分量这段代码的核心逻辑是「反射场等于加载膜系后的总场减去未加载膜系时的纯入射场」。这样处理自动消除了光源本身的频谱形状得到的R_spectrum就是膜系对入射光的反射响应的频域表示。FFT 点数取 2 的幂次是为了调用快速算法波长数组从频率数组转换而来。实际使用中需要注意total_Ez和ref_Ez的记录必须在同一个监测位置、同一时间窗口否则相位不一致会导致反射谱严重失真。3.3 透射谱计算与能量守恒校验透射谱的提取方式与反射谱类似在膜系后端设置监测点记录透射场 $E_{trans}(t)$做 FFT 后对参考频谱归一化。归一化后每个波长点的透射率 $T(\lambda)$、反射率 $R(\lambda)$ 满足T(λ) R(λ) A(λ) ≈ 1其中 $A(\lambda)$ 是吸收率。这个关系式是结果校验的黄金标准。如果膜系材料无损耗把反射谱和透射谱加起来发现明显偏离 1那说明计算域长度不够透射波还没完全离开监测点或者 PML 吸收边界反射干扰了结果。% 能量守恒校验脚本 sum_RT R_spectrum T_spectrum; % 无吸收时应接近 1 deviation abs(sum_RT - 1); bad_idx find(deviation 0.05); % 超过 5% 的波长点标记为异常 if isempty(bad_idx) fprintf(能量守恒校验通过\n); else fprintf(波长索引 %d 处 RT 偏离 1 超过 5%%\n, bad_idx); end能量守恒校验通常放在参数扫描之前。如果校验不通过先不要急着改膜系参数优先检查时间步数是否跑够、PML 层数和吸收系数是否合理。这是排查 FDTD 仿真问题最有效的起点。4. 把 25.m 与 Untitled 系列脚本串联起来参数配置和数据后处理4.1 脚本角色的推断与重组压缩包内的脚本命名并不规范25.m、Untitled.m、Untitled2.m显然不是按功能命名的。从 FDTD 仿真工作流的经验判断这四个脚本的角色大概率如下文件推断角色工作流位置fdtd_1d_3.mFDTD 主求解器含时间循环核心计算25.m参数设置脚本折射率、膜厚、波长范围前处理Untitled.m反射/透射谱后处理与绘图后处理Untitled2.m参考场计算或扫描辅助脚本辅助计算建议把它们整合成两个文件一个config.m存放所有物理参数另一个run_fdtd.m包含求解器和后处理消除散落的未命名脚本带来的路径依赖问题。整合方式是在每个脚本头部加入段注释标明输入输出变量。4.2 膜系参数化与快速验证25.m如果是一个参数配置文件通常包含膜系层数、每层厚度和折射率三个数组。多层膜设计中常见的验证对象是四分之一波堆QWOT即每层光学厚度等于设计波长的四分之一。以下是参数配置示例% 膜系参数配置高反膜设计λ0 632.8nm lambda_design 632.8e-9; % 设计波长 n_H 2.35; % TiO2 高折射率层 n_L 1.46; % SiO2 低折射率层 layer_n [n_H, n_L, n_H, n_L, n_H, n_L, n_H, n_L, n_H]; % 9层高低交替 layer_thickness (lambda_design / 4) ./ layer_n; % 每层光学厚度均为 λ0/4 % 入射与出射介质 n_inc 1.0; n_sub 1.52; % 空气入射、玻璃基底这里每层厚度取 $d \lambda_0 / (4n)$使得各界面反射波在空气中同相位叠加从而获得高反射率。此配置可在设计波长附近形成反射带带宽由折射率比值 $n_H / n_L$ 决定带宽范围 ≈ (2 / π) · arcsin[(n_H - n_L) / (n_H n_L)]当折射率对比度越大反射带越宽。这个参数配置可以直接作为fdtd_1d_3.m的输入跑完后反射谱带中心位置应与 $632.8nm$ 吻合偏差超过 $5nm$ 就需要回头检查空间步长的取值。4.3 反射透射曲线绘图与保存Untitled.m和Untitled2.m可能是两个版本的绘图后处理脚本。反射谱和透射谱的横坐标可以是波长也可以是频率光学工程师习惯用波长。绘图时最常见的失误是 FFT 频率轴顺序反了——MATLAB 的fft输出从 0 到采样率画图前要用fftshift或切片处理。% 后处理绘图反射谱与透射谱 figure(Color, w, Position, [100 100 800 500]); plot(lambda * 1e9, R_spectrum, r-, LineWidth, 1.6); hold on; plot(lambda * 1e9, T_spectrum, b--, LineWidth, 1.6); xlabel(波长 (nm)); ylabel(强度); legend(反射率, 透射率, Location, Best); xlim([0.4 1.0] * 1e3); ylim([0 1]); grid on; set(gca, FontSize, 12);绘图脚本的核心参数是xlim的范围它决定了你关心的波段窗口。把波长换算成纳米单位是为了从视觉上对齐材料手册中习惯的横轴刻度。透射谱在膜系吸收为零时与反射谱应基本互补——反射峰对应透射谷如果图像中两者趋势不一致说明监测点位置或归一化基准出了问题。5. 进阶变厚度扫描与反射谱调试的三个实用技巧5.1 用传输矩阵法做交叉验证1D FDTD 结果跑出来之后最稳妥的验证是用传输矩阵法TMM对拍。传输矩阵法把每层膜表示成一个 $2 \times 2$ 矩阵级联后直接求反射系数。对于无吸收多层膜FDTD 与 TMM 的结果在波段中心应高度一致% 传输矩阵法快速计算反射率用于对拍验证 function R_tmm tmm_reflection(layer_n, layer_d, n_inc, n_sub, lambda) M eye(2); k0 2*pi / lambda; for i 1:numel(layer_n) delta k0 * layer_n(i) * layer_d(i); % 相位厚度 M M * [cos(delta), 1j*sin(delta)/layer_n(i); ... 1j*layer_n(i)*sin(delta), cos(delta)]; end r (M(1,1)*n_inc M(1,2)*n_inc*n_sub - M(2,1) - M(2,2)*n_sub) / ... (M(1,1)*n_inc M(1,2)*n_inc*n_sub M(2,1) M(2,2)*n_sub); R_tmm abs(r)^2; endTMM 计算速度快可以在同一个波长网格上跑完整个反射谱然后与 FDTD 结果画在同一张图里对比。两者偏差在 1% 以内说明 FDTD 设置正确偏差明显则需要检查 PML 反射或时间截断。5.2 反射谱高频振荡的排查思路反射谱上出现密集的周期性振荡通常有三个原因。第一是 PML 吸收不完全反射波从边界弹回叠加形成驻波表现为谱上的高频波纹排查方法是将 PML 层数从 16 提高到 32看振荡幅度是否下降。第二是时间窗口不够时域反射波被矩形窗截断FFT 后产生泄漏旁瓣补救办法是增加total_T或对时域信号加 Hanning 窗% 对时域信号加汉宁窗抑制泄漏 win hanning(length(refl_Ez)); refl_Ez_win refl_Ez .* win;第三是空间步长偏大高频成分在介质中色散畸变后在膜系内部产生虚假谐振表现为特定波长带的非物理尖峰。此时把dx减小一半重新计算若尖峰消失即为色散问题。5.3 变厚度扫描脚本的收敛技巧在多层膜优化设计里常常需要对某一层厚度做扫描观察反射谱峰位移动规律。扫描操作可以直接在 MATLAB 循环中动态修改layer_thickness每次重建折射率分布数组后再次调用 FDTD 求解器。如果每次扫描都从零开始跑脉冲到 4000 步总时长会上升常见做法是先跑一次参考场固定不变的部分在循环外算好循环内只更新膜系区域的场。厚度扫描得到的数据用二维矩阵存储横轴是扫描厚度纵轴是波长颜色代表反射率得到反射谱随厚度演化的热力图。这种图上可以直观看出反射峰位随厚度近似线性的移动规律对比理论预期峰位正比于光学厚度变化能快速验证 FDTD 结果的物理合理性。热力图绘制可以用pcolor或imagesc注意横纵坐标都要换算成直观单位厚度用纳米、波长用纳米。本文还有配套的精品资源点击获取 用1D FDTD计算多层膜反射谱：MATLAB代码与参数解析
简介资源为基于FDTD时域有限差分方法的MATLAB仿真项目专注于多层膜结构的反射与透射特性分析可计算不同波长下的反射谱与透射谱适用于光学涂层、滤光片、太阳能电池等器件的设计与优化包括高反射膜、增透膜等常见设计场景也适合光电子、电磁场专业的课程设计和科研入门。压缩包内有4个MATLAB脚本整体仅5KB轻量易读其中包含一维FDTD核心算法实现以及网格参数、边界条件、光源和结果处理等辅助代码通过修改脚本中的膜层厚度、折射率、入射角等参数即可快速获得新的反射透射结果。目前已有517人学习下载。借助该项目不仅可以直观理解光在多层膜中的干涉机制还可将一维核心代码作为基础进一步扩展至二维/三维FDTD仿真或膜系自动优化兼具教学和科研参考价值。1. 为什么用 1D FDTD 算多层膜的反射谱做光学涂层设计的人都有这种体会传输矩阵法算一维多层膜的反射谱只要几毫秒但当膜系里出现吸收介质、非线性材料或者想观察光在膜层内部的瞬态传播过程时它就不够用了。这时 FDTD时域有限差分反而更直观它直接把麦克斯韦方程组在空间和时间上离散化光怎么进去、怎么在界面上来回反射、最后怎么透出去每一步场分布都看得见。这个fdtd-matlab.rar包里恰好就是一套用 MATLAB 写的 1D FDTD 实现四个脚本分别覆盖主求解器、参数设置和数据后处理。选择 1D FDTD 而不是 2D/3D是因为多层膜问题的场只在膜层法线方向变化横向无限延伸降维到一维后计算量骤减普通笔记本几秒钟就能跑完一个宽谱脉冲的反射谱。这套代码的价值不在于算法本身有多新而在于它把「光源 → 膜系 → 反射/透射监测」这条链路完整打通了。新手可以直接改折射率和厚度看反射谱变化老手则可以把它当成骨架往上加色散模型、斜入射和 PML 优化。下面从网格、时间步、边界条件开始逐个拆解这四个脚本的实际用途。2. 网格、时间步与边界条件1D FDTD 在多层膜中的参数骨架2.1 空间步长与 CFL 条件FDTD 的精度受限于空间步长 $\Delta x$。在多层膜仿真中膜层材料折射率通常为 1.03.5如 SiO₂ 约 1.46、TiO₂ 约 2.35、Si 约 3.5介质中的最小波长是 $\lambda_0 / n_{max}$。要获得可接受的反射谱精度通常要求 $\Delta x \le \lambda_{min} / 20$更保守的做法是取 $\lambda_{min} / 30$。时间步长 $\Delta t$ 受 CFL 条件约束在一维情况下为c · Δt / Δx ≤ 1为了保证数值稳定性工程上一般取 $S c \cdot \Delta t / \Delta x 0.5$这个值在大多数 1D FDTD 实现中足够安全还能兼顾计算效率。% 基础参数定义常见做法 lambda_ref 0.55e-6; % 参考波长单位米 n_max 3.5; % 材料最大折射率 dx lambda_ref / (n_max * 30); % 空间步长保证最密材料也有30个点 dt dx / (3e8 * 0.5); % CFL0.5安全余量 total_T 4000; % 总时间步数需覆盖脉冲在膜系中多次反射完成这段代码建立了两组核心参数空间步长决定入射光在每层膜内被离散成多少个网格点时间步长决定场推进的节奏。n_max要取整个膜系里折射率最高的材料否则高频成分在该材料内会出现明显的数值色散。total_T的取值取决于膜系总光程——膜层越厚光在里面往返所需步数越多这个值取小了反射谱会出现截断伪影。2.2 界面处的介电常数处理多层膜的每个界面上折射率发生跳变。FDTD 电场更新公式中会出现介电常数 $\varepsilon$如果直接在每个网格点上用该点所属材料的 $\varepsilon$界面处会产生非物理的微小反射。常见的做法是取相邻两个网格点介电常数的平均值也就是在界面位置做平滑处理。% 折射率分布构建界面处做算术平均 Nx 2000; % 总网格数 n ones(1, Nx) * 1.0; % 默认空气 for k 1:numel(layer_thickness) idx_start round(sum(layer_thickness(1:k-1)) / dx) 1; idx_end round(sum(layer_thickness(1:k)) / dx); n(idx_start:idx_end) layer_n(k); end n_avg (n(1:end-1) n(2:end)) / 2; % 界面平均折射率 eps n_avg.^2; % 相对介电常数这里把离散的折射率数组转换成相邻网格的平均值再用平方得到介电常数。这样处理之后电场更新方程中每个网格位置上的 $\varepsilon$ 都对应实际材料的局部平均值避免了阶梯近似造成的虚假反射。layer_thickness数组需要用户在脚本里按实际膜系填写单位与dx保持一致。2.3 PML 吸收边界与光源位置计算域两端必须设置吸收边界否则反射波碰到边界会反弹回来污染结果。一维 FDTD 中常用 PML完美匹配层在边界处放置若干层吸收介质。层数太少、吸收系数太小边界反射会叠加到反射谱上形成振荡纹波。PML 层数适用场景边界反射水平8快速验证、粗看谱形约 1% 量级16常规反射谱计算约 0.1% 量级32高精度反射率标定可低于 0.01%光源的位置也应该规划清楚。把光源放在 PML 区域外、第一个反射监测点和膜系之间入射电场沿着正向传播先到达膜系反射波再回到监测点。通常设置两个监测点一个在光源与膜系之间记录入射场加反射场另一个在膜系后端记录透射场。监测点距离膜系至少 10 个网格避免近场效应对结果的干扰。3. fdtd_1d_3.m 求解器拆解高斯脉冲激励与 FDTD 反射谱提取3.1 主循环中的电磁场交替推进fdtd_1d_3.m是这个包里编号为 3 的 1D FDTD 主求解器文件名里的_3说明它经过了至少两轮迭代修改。核心循环在每一时间步内依次更新磁场、电场再处理 PML 吸收层。以下是一个极简但完整的 1D FDTD 核心循环骨架% 场数组初始化 Ez zeros(1, Nx); % 电场 Hy zeros(1, Nx); % 磁场 % PML 导电率分布多项式渐变 sigma_max 0.75; % 最大电导率按经验取值 sigma zeros(1, Nx); pml_width 32; sigma(1:pml_width) sigma_max * ((pml_width:-1:1)/pml_width).^2; sigma(end-pml_width1:end) sigma_max * ((1:pml_width)/pml_width).^2; % 时间步进主循环 for t 1:total_T % 磁场更新Hy(n) Hy(n) - (dt/(mu*dx)) * (Ez(n1) - Ez(n)) Hy(1:end-1) Hy(1:end-1) - (dt / (mu0 * dx)) * (Ez(2:end) - Ez(1:end-1)); % 电场更新含 PML 损耗项 Ez(2:end-1) Ez(2:end-1) - (dt ./ eps(2:end-1) / dx) .* (Hy(2:end-1) - Hy(1:end-2)); % PML 内电场衰减 Ez(1:pml_width) Ez(1:pml_width) ./ (1 sigma(1:pml_width) * dt / eps(1:pml_width)); Ez(end-pml_width1:end) Ez(end-pml_width1:end) ./ (1 sigma(end-pml_width1:end) * dt / eps(end-pml_width1:end)); % 高斯脉冲源注入硬源或柔软源 t0 60; tau 15; source exp(-((t - t0)/tau)^2); Ez(src_pos) source; end磁场与电场在时间上差半步更新这是 FDTD 的 Yee 网格格式的基本特点保证了二阶精度。PML 层的电导率从边界向内渐变吸收系数过大会产生数值反射过小则吸收不充分sigma_max 0.75是一维问题中比较稳的起点值。高斯脉冲的tau直接决定频谱宽度tau越小脉冲越窄、频谱越宽能覆盖的波长范围越大但需要更小的空间步长来保证高频成分不失真。3.2 反射谱的频域提取反射谱不能在时域直接看需要把监测点记录的时域电场做傅里叶变换再用反射场频谱除以入射场频谱得到反射系数。这里有三种做法先跑一次无膜系纯空气的仿真记录参考波 $E_{ref}(t)$再跑有膜系的仿真记录 $E_{total}(t)$相减得到反射波 $E_{refl}(t)$然后分别做 FFT 求比值。在膜系两侧各设一个监测点分别得到反射场和透射场直接用透射频谱除以参考频谱得到透射率。用连续波扫描一个波长一个波长地算稳态反射率精度最高但耗时。方案 1 是实际最常用的因为一次性宽谱激励就能覆盖整个目标波段代码如下% 监测点时域记录假设已经同时记录参考场 ref_Ez 和总场 total_Ez refl_Ez total_Ez - ref_Ez; % 反射场 总场 - 入射场 % FFT 提取频域响应 Fs 1 / dt; % 采样率 N_FFT 2^nextpow2(length(refl_Ez)); freq (0:N_FFT-1) / N_FFT * Fs; lambda 3e8 ./ freq(2:end); % 频率转波长 R_spectrum abs(fft(refl_Ez, N_FFT)) ./ abs(fft(ref_Ez, N_FFT)); R_spectrum R_spectrum(2:end); % 去掉直流分量这段代码的核心逻辑是「反射场等于加载膜系后的总场减去未加载膜系时的纯入射场」。这样处理自动消除了光源本身的频谱形状得到的R_spectrum就是膜系对入射光的反射响应的频域表示。FFT 点数取 2 的幂次是为了调用快速算法波长数组从频率数组转换而来。实际使用中需要注意total_Ez和ref_Ez的记录必须在同一个监测位置、同一时间窗口否则相位不一致会导致反射谱严重失真。3.3 透射谱计算与能量守恒校验透射谱的提取方式与反射谱类似在膜系后端设置监测点记录透射场 $E_{trans}(t)$做 FFT 后对参考频谱归一化。归一化后每个波长点的透射率 $T(\lambda)$、反射率 $R(\lambda)$ 满足T(λ) R(λ) A(λ) ≈ 1其中 $A(\lambda)$ 是吸收率。这个关系式是结果校验的黄金标准。如果膜系材料无损耗把反射谱和透射谱加起来发现明显偏离 1那说明计算域长度不够透射波还没完全离开监测点或者 PML 吸收边界反射干扰了结果。% 能量守恒校验脚本 sum_RT R_spectrum T_spectrum; % 无吸收时应接近 1 deviation abs(sum_RT - 1); bad_idx find(deviation 0.05); % 超过 5% 的波长点标记为异常 if isempty(bad_idx) fprintf(能量守恒校验通过\n); else fprintf(波长索引 %d 处 RT 偏离 1 超过 5%%\n, bad_idx); end能量守恒校验通常放在参数扫描之前。如果校验不通过先不要急着改膜系参数优先检查时间步数是否跑够、PML 层数和吸收系数是否合理。这是排查 FDTD 仿真问题最有效的起点。4. 把 25.m 与 Untitled 系列脚本串联起来参数配置和数据后处理4.1 脚本角色的推断与重组压缩包内的脚本命名并不规范25.m、Untitled.m、Untitled2.m显然不是按功能命名的。从 FDTD 仿真工作流的经验判断这四个脚本的角色大概率如下文件推断角色工作流位置fdtd_1d_3.mFDTD 主求解器含时间循环核心计算25.m参数设置脚本折射率、膜厚、波长范围前处理Untitled.m反射/透射谱后处理与绘图后处理Untitled2.m参考场计算或扫描辅助脚本辅助计算建议把它们整合成两个文件一个config.m存放所有物理参数另一个run_fdtd.m包含求解器和后处理消除散落的未命名脚本带来的路径依赖问题。整合方式是在每个脚本头部加入段注释标明输入输出变量。4.2 膜系参数化与快速验证25.m如果是一个参数配置文件通常包含膜系层数、每层厚度和折射率三个数组。多层膜设计中常见的验证对象是四分之一波堆QWOT即每层光学厚度等于设计波长的四分之一。以下是参数配置示例% 膜系参数配置高反膜设计λ0 632.8nm lambda_design 632.8e-9; % 设计波长 n_H 2.35; % TiO2 高折射率层 n_L 1.46; % SiO2 低折射率层 layer_n [n_H, n_L, n_H, n_L, n_H, n_L, n_H, n_L, n_H]; % 9层高低交替 layer_thickness (lambda_design / 4) ./ layer_n; % 每层光学厚度均为 λ0/4 % 入射与出射介质 n_inc 1.0; n_sub 1.52; % 空气入射、玻璃基底这里每层厚度取 $d \lambda_0 / (4n)$使得各界面反射波在空气中同相位叠加从而获得高反射率。此配置可在设计波长附近形成反射带带宽由折射率比值 $n_H / n_L$ 决定带宽范围 ≈ (2 / π) · arcsin[(n_H - n_L) / (n_H n_L)]当折射率对比度越大反射带越宽。这个参数配置可以直接作为fdtd_1d_3.m的输入跑完后反射谱带中心位置应与 $632.8nm$ 吻合偏差超过 $5nm$ 就需要回头检查空间步长的取值。4.3 反射透射曲线绘图与保存Untitled.m和Untitled2.m可能是两个版本的绘图后处理脚本。反射谱和透射谱的横坐标可以是波长也可以是频率光学工程师习惯用波长。绘图时最常见的失误是 FFT 频率轴顺序反了——MATLAB 的fft输出从 0 到采样率画图前要用fftshift或切片处理。% 后处理绘图反射谱与透射谱 figure(Color, w, Position, [100 100 800 500]); plot(lambda * 1e9, R_spectrum, r-, LineWidth, 1.6); hold on; plot(lambda * 1e9, T_spectrum, b--, LineWidth, 1.6); xlabel(波长 (nm)); ylabel(强度); legend(反射率, 透射率, Location, Best); xlim([0.4 1.0] * 1e3); ylim([0 1]); grid on; set(gca, FontSize, 12);绘图脚本的核心参数是xlim的范围它决定了你关心的波段窗口。把波长换算成纳米单位是为了从视觉上对齐材料手册中习惯的横轴刻度。透射谱在膜系吸收为零时与反射谱应基本互补——反射峰对应透射谷如果图像中两者趋势不一致说明监测点位置或归一化基准出了问题。5. 进阶变厚度扫描与反射谱调试的三个实用技巧5.1 用传输矩阵法做交叉验证1D FDTD 结果跑出来之后最稳妥的验证是用传输矩阵法TMM对拍。传输矩阵法把每层膜表示成一个 $2 \times 2$ 矩阵级联后直接求反射系数。对于无吸收多层膜FDTD 与 TMM 的结果在波段中心应高度一致% 传输矩阵法快速计算反射率用于对拍验证 function R_tmm tmm_reflection(layer_n, layer_d, n_inc, n_sub, lambda) M eye(2); k0 2*pi / lambda; for i 1:numel(layer_n) delta k0 * layer_n(i) * layer_d(i); % 相位厚度 M M * [cos(delta), 1j*sin(delta)/layer_n(i); ... 1j*layer_n(i)*sin(delta), cos(delta)]; end r (M(1,1)*n_inc M(1,2)*n_inc*n_sub - M(2,1) - M(2,2)*n_sub) / ... (M(1,1)*n_inc M(1,2)*n_inc*n_sub M(2,1) M(2,2)*n_sub); R_tmm abs(r)^2; endTMM 计算速度快可以在同一个波长网格上跑完整个反射谱然后与 FDTD 结果画在同一张图里对比。两者偏差在 1% 以内说明 FDTD 设置正确偏差明显则需要检查 PML 反射或时间截断。5.2 反射谱高频振荡的排查思路反射谱上出现密集的周期性振荡通常有三个原因。第一是 PML 吸收不完全反射波从边界弹回叠加形成驻波表现为谱上的高频波纹排查方法是将 PML 层数从 16 提高到 32看振荡幅度是否下降。第二是时间窗口不够时域反射波被矩形窗截断FFT 后产生泄漏旁瓣补救办法是增加total_T或对时域信号加 Hanning 窗% 对时域信号加汉宁窗抑制泄漏 win hanning(length(refl_Ez)); refl_Ez_win refl_Ez .* win;第三是空间步长偏大高频成分在介质中色散畸变后在膜系内部产生虚假谐振表现为特定波长带的非物理尖峰。此时把dx减小一半重新计算若尖峰消失即为色散问题。5.3 变厚度扫描脚本的收敛技巧在多层膜优化设计里常常需要对某一层厚度做扫描观察反射谱峰位移动规律。扫描操作可以直接在 MATLAB 循环中动态修改layer_thickness每次重建折射率分布数组后再次调用 FDTD 求解器。如果每次扫描都从零开始跑脉冲到 4000 步总时长会上升常见做法是先跑一次参考场固定不变的部分在循环外算好循环内只更新膜系区域的场。厚度扫描得到的数据用二维矩阵存储横轴是扫描厚度纵轴是波长颜色代表反射率得到反射谱随厚度演化的热力图。这种图上可以直观看出反射峰位随厚度近似线性的移动规律对比理论预期峰位正比于光学厚度变化能快速验证 FDTD 结果的物理合理性。热力图绘制可以用pcolor或imagesc注意横纵坐标都要换算成直观单位厚度用纳米、波长用纳米。本文还有配套的精品资源点击获取