简介这是一份用于连续抽运被动调Q激光脉冲数值仿真的代码面向激光物理、光电子及相关专业学生与科研人员针对可饱和吸收体激发态吸收的影响建立了完整的被动调Q速率方程组并通过常微分方程初值问题求解函数对耦合方程进行数值求解从而模拟调Q脉冲序列的建立与演化过程。压缩包共3个文件均为脚本整体仅2千字节包含速率方程定义、主仿真程序以及结果分析脚本文件划分清晰代码结构紧凑便于学习者按需修改和快速移植。该资源已有447人学习下载。运行代码后能够直观获得被动调Q激光脉冲的数值结果进一步提取脉冲重复频率和脉冲宽度等关键参数为理解调Q机理、完成课程设计或开展相关仿真研究提供了可直接使用的工具与参考实现同时便于进一步探索不同抽运功率与可饱和吸收体参数下的脉冲特性适合具备激光原理基础并希望动手进行数值仿真的读者。1. 连续抽运被动调Q可饱和吸收体的重复频率为什么不能靠解析公式连续抽运下的被动调Q激光器输出脉冲序列的重复频率和脉冲宽度不是由某个单一公式决定的而是由增益介质的储能、可饱和吸收体的初始透过率、腔内光子寿命以及抽运速率共同耦合出来的瞬态过程。很多人上来就想套用谐振腔的稳态条件结果算出的是一个平均功率而不是脉冲行为。beidongtiaoQ.zip 这组 MATLAB 代码给出了更实际的路径用rate_eq.m建立速率方程组用Qswitch.m调用 ode23() 对脉冲微分方程做初值求解再用Qswitch_Analyze.m从光子密度时间序列里提取重复频率和脉宽。对做激光器仿真的光学工程师或者需要标定 Nd:YAG/Cr4:YAG 这类被动调Q系统的开发人员来说这套思路比反复调参数试硬件要快得多。2. 可饱和吸收体被动调Q的速率方程激发态吸收项怎么写2.1 三变量与两变量的取舍被动调Q的基本过程并不复杂抽运持续为上能级积累粒子数可饱和吸收体的基态吸收让腔内损耗保持在高位光子密度被压制当反转粒子数超过阈值一个涨落光子触发了雪崩式的受激辐射光子密度在几个腔往返时间内上升到峰值同时把吸收体漂白。问题在于这个过程中吸收体并不是一个简单的“有损/漂白”两态开关而是存在基态和激发态两个布居数且激发态对腔内光子也有吸收截面。所以在速率方程里不能只写一个可饱和吸收体的总粒子数守恒而要把吸收体的基态粒子数密度作为独立变量。加上光子密度和增益介质的反转粒子数密度至少需要三个一阶常微分方程。如果忽略激发态吸收则在脉冲峰值附近损耗会被低估计算出来的脉冲宽度会比实测偏窄重复频率也会因为损耗恢复过快而偏高。rate_eq.m里最常见的结构就是把光子密度、增益粒子数密度、吸收体基态粒子数密度放在列向量y中返回它们的导数。三方程的好处是可以直接体现基态吸收和激发态吸收分别对腔损耗的贡献也方便后续把吸收体换成半导体可饱和吸收镜时保留同一个框架。2.2 rate_eq.m 的右端函数与参数约定在beidongtiaoQ.zip中rate_eq.m作为 ode23() 的右端函数输入是当前时刻、状态向量以及结构体p。下面是一个在工程上拿过来就能改的写法% rate_eq.m — 连续抽运被动调Q速率方程组右端函数 % y(1): 腔内光子密度 % y(2): 增益介质反转粒子数密度 % y(3): 可饱和吸收体基态粒子数密度 function dydt rate_eq(t, y, p) phi y(1); n_gain y(2); n_gsa y(3); n_esa p.n_sat - n_gsa; % 激发态吸收体粒子数密度 tr 2 * p.L_cav / p.c0; % 腔内光子往返时间 gain 2 * p.sigma_g * n_gain * p.l_g; loss_sa 2 * p.sigma_gs * n_gsa * p.l_s ... 2 * p.sigma_es * n_esa * p.l_s; loss_mir -log(p.R_out) p.L_int; dphi phi / tr * (gain - loss_sa - loss_mir); % 抽运项 自发辐射上能级衰减 受激辐射消耗 dn_gain p.Rp - n_gain / p.tau_g - p.sigma_g * p.c0 * phi * n_gain; % 基态被漂白激发态通过寿命恢复回基态 dn_gsa -p.sigma_gs * p.c0 * phi * n_gsa n_esa / p.tau_s; dydt [dphi; dn_gain; dn_gsa];这里的p是参数结构体而不是解向量ode23() 会把它透传给右端函数。dphi使用了单个光子密度增长率公式括号里的gain - loss_sa - loss_mir是净增益除以往返时间loss_mir已经含了输出耦合的自然对数项。dn_gsa中受激吸收带走基态粒子数激发态通过寿命项恢复这个恢复速度直接决定了同一抽运周期内还能不能产生第二个脉冲。下面这张表是 Nd:YAG / Cr4:YAG 这类典型系统的初始参数范围也是我在调rate_eq.m时最常改的几个量参数物理含义典型量级sigma_g增益介质受激发射截面2.8e-23 m^2sigma_gs饱和吸收体基态吸收截面4.0e-22 m^2sigma_es饱和吸收体激发态吸收截面9.0e-23 m^2tau_g增益介质上能级寿命230e-6 stau_s饱和吸收体基态恢复时间3.5e-6 sn_sat饱和吸收体总掺杂浓度4e23 m^-3R_out输出镜反射率0.85L_int腔内线性散射损耗0.02Rp连续抽运速率1e22 m^-3 s^-1参数结构体p.c0最好写成增益介质内的等效光速而不是真空光速否则在调整腔长时会把折射率效应重复计入。Rp不写成功率而写成抽运速率是因为方程里每一项的单位必须是同一个数量级如果把抽运功率直接除到有效模体积里也可以但要保证n_gain和phi都定义在同一个模体积上。3. 用 ode23() 解脉冲微分方程Qswitch.m 与事件函数的配合3.1 为什么是 ode23 而不是直接循环被动调Q的脉冲微分方程既不是严格刚性的也不是全程平滑的。脉冲上升沿只有几十纳秒而抽运阶段的弛豫时间在微秒到百微秒量级。ode23() 在误差控制上比定步长 RK4 灵活又比 ode15s 在非刚性段开销小。对Qswitch.m来说真正要注意的不是求解器选型而是“如何让一个长时间仿真覆盖几十个脉冲”。常见做法是使用事件函数把连续时间区间切成一段段脉冲。每次 ode23() 解到一个脉冲结束就停止一次记录该脉冲的峰值位置和峰后状态然后把光子密度重置到一个极低值再继续积分下一个抽运周期。这么做是因为被动调Q的脉冲波形宽度与抽运周期相差 4 个数量级如果一次性用一个大时间跨度的tspanode23() 的步长控制会在脉冲快速上升段反复缩小步长最终把大部分计算时间消耗在寻找脉冲波形上而且容易漏掉峰值处的极窄尖峰。Qswitch.m的框架通常是这样% Qswitch.m — 分段调用ode23求解被动调Q脉冲序列 options odeset(RelTol, 1e-6, AbsTol, 1e-8, ... Events, pulse_detect, ... MaxStep, p.tau_g / 20); t_current 0; t_end 5e-3; y [1e-10; 1e20; p.n_sat]; t_collect []; y_collect []; pulse_times []; while t_current t_end [t, y, te, ye, ie] ode23(rate_eq, ... [t_current, t_end], ... y, options, p); t_collect [t_collect; t(2:end)]; y_collect [y_collect; y(2:end, :)]; if isempty(te) break; end % 事件返回的te是本次脉冲结束时刻 pulse_times(end1, 1) te(1); y [1e-10; ye(1,2); ye(1,3)]; t_current te(1); endte是事件发生时刻ye是事件发生时的状态向量。这里把光子密度重置成1e-10不是为了表示真实物理而是避免刚过完脉冲的小幅度弛豫振荡再次触发事件函数从而把下一个脉冲的涨落提前。实际光子密度从零涨到阈值需要一个抽运积累期这正好被后续 ode23() 的正常积分模拟出来。事件函数本身不复杂关键是direction的设置% 事件检测光子密度回落穿过阈值时停住 function [value, isterminal, direction] pulse_detect(t, y, p) value y(1) - p.peak_threshold; isterminal 1; direction -1; % 只在下降沿触发direction -1让事件只在光子密度由高于阈值变为低于阈值时触发这样稳态噪声不会被当成脉冲。p.peak_threshold应该至少比最大光子密度低 3 个数量级否则一个正常脉冲结束后光子密度还没降到阈值以下就被迫续算下一次事件会在同一个脉冲的回摆中反复触发。Qswitch.m里这几个选项的取值直接决定整段仿真是否收敛odeset 选项推荐值作用RelTol1e-6相对误差容限太小会让积分变慢AbsTol1e-8绝对误差容限防止光子密度接近0时步长过小MaxSteptau_g/20限制抽运阶段最大步长避免跨过脉冲建立点Eventspulse_detect每个脉冲结束后停止本次积分3.2 从时间序列中提取脉冲序列Qswitch.m输出的是一个连续的时间序列而不是像增益开关那种单脉冲解。Qswitch_Analyze.m拿到t_collect和y_collect后先找每个脉冲的边界再计算脉冲峰值位置和半高宽。这部分我用的是阈值标记法比直接调用findpeaks更可控也不依赖额外工具箱% Qswitch_Analyze.m — 提取重复频率和脉冲宽度 phi y_collect(:, 1); ts t_collect(:, 1); base max(phi) * 1e-4; mask (phi base); idx_start find(diff([0; mask]) 1); idx_stop find(diff([0; mask]) -1); for k 1:length(idx_start) seg phi(idx_start(k):idx_stop(k)); tseg ts(idx_start(k):idx_stop(k)); [pk, idx] max(seg); peak_time(k) tseg(idx); half pk / 2; ileft find(seg(1:idx) half, 1, first); iright find(seg(idx:end) half, 1, first) idx - 1; fwhm_s(k) tseg(iright) - tseg(ileft); end rep_time diff(peak_time); rep_freq 1 ./ rep_time;这段代码先按阈值把光子密度序列切成每个脉冲的区间再在单脉冲区间内找最大值。ileft和iright分别是半高位置索引fwhm_s的单位与时间向量一致。要注意的是启动阶段的前几个脉冲峰值时间间隔会与稳定阶段差很多直接mean(rep_freq)会低估重频后面我会说怎么处理这个瞬态偏差。4. 重复频率与脉冲宽度的参数标定从速率方程到可复现数值4.1 参数标定的先后顺序把rate_eq.m里的参数直接填进仿真大部分时候结果不收敛问题不在算法而在量纲。比如 Cr4:YAG 的小信号透过率是由sigma_gs、n_sat、l_s三者乘积决定的而脉冲宽度主要由增益介质的储能和腔内光子寿命决定重复频率则由抽运速率和上能级寿命共同决定。这三个维度在数值上相互耦合但实际调参时应按“先定损耗、再定存储、最后定抽运”的顺序调整。先把sigma_es设成sigma_gs的 1/4 左右因为大多数实用可饱和吸收体的激发态吸收截面小于基态吸收截面但又不至于小到可以忽略。然后调loss_mir使光子密度在小信号状态下保持不增长把Rp设成零时phi必须单调衰减且衰减时间常数等于tr / loss_mir。如果光子密度在无抽运时反而增长说明R_out对应的损耗与增益介质残留吸收的叠加算错了。在beidongtiaoQ.zip的代码里Qswitch.m和rate_eq.m的参数互相引用时容易出现变量未初始化问题。我一般会在引用p结构体前用一句p struct(c0, 3e8, L_cav, 0.2);兜底防止 MATLAB 第一次执行时把某个字段误判成函数名。4.2 时间单位统一与刚性切换速率方程里时间基准很混乱脉冲上升沿用纳秒描述抽运阶段用微秒描述。如果全部使用秒作为单位tau_g 230e-6和tau_s 3.5e-6在 ode23() 的误差控制里并没有问题有问题的是AbsTol取太小会让积分器在光子密度接近零时强制缩小步长。更工程化的做法是把时间刻度统一换算到微秒纳秒量级比如外部用秒做输出单位但在微分方程内部把tau_g除以 1e6 等比例缩放。这样做能让 ode23() 的相对误差分布更均匀。如果换成 ode23s 或 ode15s则要额外注意MaxStep不能太大否则刚性求解器会用大步长直接略过脉冲尖峰。这里给出一行可比的切换方式options odeset(options, MaxOrder, 2); [t, y] ode23s(rate_eq, tspan, y0, options, p);ode23s 对刚性很强的速率方程更稳但多数被动调Q场景并不是真刚性而是“局部刚性”也就是只有在脉冲上升沿附近的 Jacobian 特征值差距很大。直接换 ode23s 会让非刚性段计算时间增加 30% 到 50%因此我建议先用 ode23() 跑通只有遇到脉冲宽度抖动、步长频繁报警时才切换。4.3 仿真失败的常见信号第一个典型问题是“脉冲宽度不随抽运速率变化”。这通常不是求解器错误而是dphi里缺少光子密度在达到峰值后的快速下降驱动力。被动调Q的脉冲宽度由可饱和吸收体恢复时间和光子寿命共同决定如果tau_s设置得过长吸收体来不及在脉冲尾部重新吸收光子脉冲会被拉得很扁。调低tau_s至少 5 倍看fwhm_s是否收敛。第二个典型问题是“重复频率跳变”表现为相邻两个峰值间隔忽大忽小。我遇到过的情况大多是事件阈值peak_threshold太低导致同一个脉冲的下降沿噪声被识别成新事件。解决办法是把阈值提高到最大光子密度的 1e-5 到 1e-4 量级并在Qswitch_Analyze.m里加一个最小峰值间隔约束比如要求相邻peak_time的间隔不得小于tau_g / 5小于这个值就合并为同一个脉冲。出现第三类问题ode23()在 while 循环里卡死多半是重置后的光子密度1e-10还没到达上升沿就不断被事件检测判定为“下降沿”。此时可以把事件函数里的方向条件改成direction 0并配合isterminal 1这样不管是上升还是下降穿过阈值都会截断但需要在截断后检查ie的值只在下降沿记录脉冲。5. 用 Qswitch_Analyze.m 把脉冲宽度和重复频率算准的三种取巧方式5.1 重频取中位数而不是均值连续抽运被动调Q在刚开始的几个脉冲之间重复频率会比稳态值高。这是因为初始反转粒子数接近零第一个脉冲实际上消耗了更多存储能量恢复周期更长这一瞬态过程会拉低均值频率。分析Qswitch_Analyze.m里得到的rep_freq时只用后 50% 的数据并且用median(rep_freq(end-20:end))来估值。这样得到的重复频率对离散毛刺不敏感也更接近实验里用频率计测得的结果。5.2 脉冲宽度用插值代替最近点索引上面给出的fwhm_s直接用时间索引求差误差会受 ode23() 自适应步长影响。脉冲宽度通常只有几十纳秒而抽运阶段步长可能到几百纳秒索引法算出来的半宽经常是步长整数倍伪影。解决方式是用interp1对单个脉冲波形做三次插值后再找半高交点t_fine linspace(tseg(1), tseg(end), 1000); phi_fine interp1(tseg, seg, t_fine, pchip); half max(phi_fine) / 2; cross_down find(phi_fine half, 1, first); fwhm_interp (t_fine(cross_down) - t_fine(1)) * 2;这段代码假设脉冲波形左右对称且从波形上升起点到落下到半高的时间近似等于半宽的一半。对于对称度高的时候这个近似误差可以忽略。对不太对称的脉冲改用左右分别找半高点的做法先找到max(phi_fine)的时间索引再向左找最后一个大于等于半高的点向右找第一个小于等于半高的点两者做差结果更稳。5.3 把关键量落盘再对照实验同一组仿真参数多跑几遍会发现重复频率对Rp呈近似线性脉冲宽度对Rp则是一条缓慢下降的曲线。与其在 MATLAB 里反复打断看图不如让Qswitch_Analyze.m直接把这些量写出来writematrix([peak_time(:), rep_freq(:), fwhm_interp(:)], ... passive_q_switch_output.csv);CSV 导出到 Excel 或 Origin 后再画成曲线既能直观看出脉冲序列的建立过程又能把启动段和稳态段分别选中做线性拟合。对照示波器读数时把横轴单位从秒换算成纳秒峰值光强单位换算成采样电压值就能以极小的计算成本验证速率方程里的tau_s、sigma_es这些参数是否和实际晶体一致。本文还有配套的精品资源点击获取