企业官网建设流程全解析

简介：这份MATLAB仿真资料面向学习固体激光器原理与数值建模的高校学生、科研入门者及工程技术人员，围绕主动调Q固体激光器的高功率短脉冲输出过程，提供可运行的数值实验工具，用于理解四能级系统、粒子数反转与声光调Q机制。压缩包为rar格式，共2个文件，均为m脚本，体积约919B，轻量便于在MATLAB中直接打开、修改和复现。已有2102人学习，说明其在激光器仿真入门场景中具有一定参考价值。内容以四能级速率方程为主线，可调整增益、损耗、泵浦功率及声光调制器开关时序，观察脉冲宽度、峰值功率和能量效率的变化；同时配有被动调Q例程作对照，便于比较主动与被动调Q在脉冲形成和参数敏感性上的差异，适合课程实验、自学推演与设计优化前的方案验证。

1. 从示波器上那条“对不上”的波形说起

调 Q 波形到底长什么样，很多时候不是算不出来，是算出来的那一版不敢信。前阵子帮一个做固体激光器的组看主动调 Q 的 MATLAB 仿真，他们卡在一个特别典型的点上：泵浦时间、重复频率、Q 开关开门时刻这三个参数怎么配，示波器上的脉宽和峰值功率就跟仿真对不上，反复调也调不拢。这套 matlab 固体激光器仿真文件解决的正是这件事：用四能级速率方程，把 Nd:YAG 这类增益介质、声光或电光调 Q 的开关时序、腔内损耗和输出耦合串成一条完整链路，跑出脉宽、峰值功率、单脉冲能量随参数变化的曲线。适合两类人：手里有激光器但在调参上反复试错的实验同学，以及想先把物理图像搞明白、再决定上不上实验台的研究生。它不能替代 Zemax 或 Fox-Li 那类横向模式计算，但它能把“开门早一点还是晚一点”“输出耦合从 10% 提到 30% 能量掉多少”这种问题，用几秒钟回答一次。

2. 项目定位与技术选型：速率方程求解器比全波仿真更适合调 Q 场景

2.1 这套仿真到底解的是什么物理问题

主动调 Q 的物理过程可以切成三段。第一段是泵浦积累：Q 开关处于关断状态，腔内附加损耗把阈值抬得很高，反转粒子数只能一路往上存，存到远高于正常阈值也不敢起振。第二段是开门：损耗在几十纳秒内塌下来，腔内那几个自发辐射光子开始被放大，增益远大于损耗。第三段是雪崩与泄放：光子数暴涨，反转粒子数被迅速抽干，增益跌回阈值以下，脉冲结束，剩下的反转粒子数进入下一轮。

速率方程模型就是给这三段写两组常微分方程，一组描述反转粒子数密度 N₂ 的收支，一组描述腔内光子数密度 φ 的收支。它做了一个很强的假设：增益介质内部的反转粒子数和光子密度是空间均匀的，只保留时间维度。这个假设换来的是计算量下降好几个数量级。上一套 Fox-Li 加迭代的办法，算一次驻波场分布要几分钟，而速率方程跑一次完整的调 Q 过程，在普通笔记本上是秒级。做参数扫描、调重复频率、试不同输出耦合，靠的就是这个速度。

代价也要说清楚：空间平均意味着它算不出横模结构、算不出热透镜导致的模式体积变化、也算不出增益介质里的纵向反转分布不均。它给你的是一条时间波形和几个积分量，不给光斑。想清楚这一点，后面参数怎么设就有数了。

2.2 文件与目录结构：五个 .m 文件各管一段

这份仿真文件的组织方式很朴素，没有 App Designer，也没有 Simulink 模型，就是一组函数加脚本，按职责切开：

q_switch_laser/ ├── main.m % 入口：选单段还是两段式，调求解器，出图 ├── params_default.m % 全部物理参数集中在这里，只返回一个结构体 ├── rate_eq.m % 速率方程右端项，被求解器反复调用 ├── q_loss.m % Q 开关损耗随时间的变化（矩形/高斯/自定义） ├── post_process.m % 光子密度 -> 输出功率 -> 能量/峰值/脉宽 └── scan/ ├── scan_repeat_freq.m % 扫重复频率 └── scan_coupling.m % 扫输出耦合透过率

这种拆法的好处是改一处不影响别处。比如你想把理想矩形开关换成带有限上升沿的高斯型，只需要动 q_loss.m，rate_eq.m 里那行调用不用碰。想把四能级换成三能级，改 params_default.m 里的 gamma 和方程里的一项就行，后处理完全不受影响。

我一般会额外加一个results/目录存 .mat 和 .png，脚本跑完自动写进去。参数扫描动辄几十组，不落盘的话中途崩了就得重来，这个亏吃过一次就够了。

2.3 技术栈边界：只用 base MATLAB 能走到哪一步

整个工程不依赖任何工具箱，用的都是ode15s、odeset、trapz、interp1这些基础函数。这一点很重要，因为很多实验室的 MATLAB 授权并不全，跑一个仿真还要去申请 Optimization Toolbox 是很烦的事。R2016b 之后的版本可以直接用，因为代码里用了匿名函数把参数结构体塞进求解器，老版本需要把@(t,y) rate_eq(t,y,p)拆成一个带全局变量的独立函数文件。

数值方法上，速率方程在开门瞬间是典型的刚性方程：反转粒子数的时间尺度是上能级寿命（几百微秒），光子数的时间尺度是腔内往返时间（纳秒量级），两者差了四五个数量级。用 ode45 硬跑，要么步长被压到极小、跑半小时出不来，要么步长跟不上、光子数直接跑成负数。所以默认用 ode15s，这是变阶的刚性求解器，专门对付这种快慢分离的系统。

需要注意的一点：这份文件里没有做横向模式，也没有做热效应。如果你的激光器是高重频、高平均功率，热透镜焦距可能从几百毫米缩到几十毫米，模式面积变化会显著改变阈值和脉冲能量。这种情况要么在外部先用 ABCD 矩阵算出模式半径，再作为参数喂进来，要么就得换含横向维度的模型，那是另一个量级的工程。

3. 核心模块拆解：从泵浦到脉冲输出的一条完整计算链路

3.1 参数集中管理：一个结构体装下所有物理量

参数散落在各个脚本里是仿真代码最常见的烂尾原因。这份文件用一个函数返回结构体，好处是扫描参数时可以整体复制、只改一个字段，不会漏改。默认值按 Nd:YAG 的量级给，全部用 SI 单位：

function p = params_default() % 主动调Q固体激光器仿真默认参数（四能级、Nd:YAG 量级） % 统一 SI 单位：长度 m，时间 s，面积 m^2，粒子数密度 m^-3 % ---- 增益介质 ---- p.sigma = 2.8e-23; % 受激发射截面 m^2（即 2.8e-19 cm^2） p.tau = 230e-6; % 上能级寿命 s p.gamma = 1.0; % 反转因子：四能级取 1，三能级取 2 p.n_ref = 1.82; % 折射率 p.l = 0.01; % 增益介质长度 m p.A = pi*(0.003)^2; % 模式面积 m^2，约对应 Φ3 mm 光斑 p.V = p.A * p.l; % 模式体积 m^3 % ---- 谐振腔 ---- p.Lcav = 0.15; % 几何腔长 m p.Lopt = p.Lcav + (p.n_ref-1)*p.l; % 光学腔长 m p.t_r = 2*p.Lopt/3e8; % 往返时间 s p.T = 0.10; % 输出耦合透过率 p.loss_int = 0.02; % 腔内往返寄生损耗 % ---- Q 开关 ---- p.loss_hold = 0.60; % 关断期附加损耗（声光开关关断衍射损耗） p.t_open = 200e-6; % 开门时刻 s % ---- 泵浦 ---- p.Rp = 5e27; % 泵浦速率 m^-3 s^-1 p.beta = 1e-6; % 自发辐射耦合因子，用来"点火" % ---- 仿真控制 ---- p.t_end = 260e-6; % 单段模式的总时长 s p.two_stage= 1; % 1 = 泵浦段解析、开门段数值，推荐 p.win = 200e-9; % 开门段积分窗口 s end

gamma这个参数在四能级系统里取 1，原因是激光下能级基本是空的，粒子数反转就等于上能级粒子数；三能级系统下能级有热布居，要取 2 来补偿。写成参数而不是写死在方程里，是因为有人会拿去试 Nd:YVO₄ 或者掺铒体系，换个值就能跑。

t_r是往返时间，光学腔长已经把介质折射率算进去了，这一点容易漏。如果只拿几何腔长算，往返时间会偏小几个百分点，光子寿命跟着偏，脉冲宽度就对不上。

3.2 速率方程组：两组常微分方程和它们的量纲

右端项函数是整个工程的心脏。四能级速率方程写成下面这个样子，每一项的单位都必须是 m⁻³s⁻¹：

function dydt = rate_eq(t, y, p) % y(1) = N2 反转粒子数密度 m^-3 % y(2) = phi 腔内光子数密度 m^-3 N2 = y(1); phi = y(2); % 光子寿命随 Q 开关状态切换，这是"主动"二字的全部含义 loss = q_loss(t, p) + p.T + p.loss_int; t_c = p.t_r / loss; % 反转粒子数收支：泵浦 - 自发辐射 - 受激辐射 dN2 = p.Rp ... - N2 / p.tau ... - p.gamma * p.sigma * 3e8 * N2 * phi; % 光子数收支：受激放大 - 腔内损耗 + 自发辐射种子 dphi = (2 * p.sigma * N2 * p.l / p.t_r) * phi ... - phi / t_c ... + p.beta * N2 / p.tau; dydt = [dN2; dphi]; end

逐项对一遍量纲。p.Rp单位是 m⁻³s⁻¹，N2/p.tau是 m⁻³ 除以 s，同样 m⁻³s⁻¹。第三项sigma * c * N2 * phi：m²·(m/s)·m⁻³·m⁻³ = m⁻³s⁻¹，对得上。光子方程里2*sigma*N2*l/t_r是无量纲的增益系数（每次往返的净增益），乘上 φ 得到 m⁻³s⁻¹；phi/t_c是损耗项，beta*N2/tau是自发辐射往激光模式里漏的那一点点，量级极小但必不可少——没有它，初值给 0 光子，方程永远解出 0。

2*sigma*N2*l里的 2 是驻波场的空间烧孔补偿因子。如果你的腔是环形腔行波运转，这个 2 要去掉，去掉之后阈值和脉冲能量都会变，这是改模型时最容易忘掉的一处。

q_loss.m把开关状态折成损耗：

function loss = q_loss(t, p) % 理想矩形开关：开门前高损耗，开门后归零 % 想模拟有限上升沿，把下面换成 erfc 或分段线性即可 loss = p.loss_hold * (t < p.t_open); end

一行就够。真实的声光开关上升沿在几十到一百纳秒，如果你的脉宽本身只有十几纳秒，用理想阶跃会低估脉宽、高估峰值功率，这时候把(t < p.t_open)换成0.5*erfc((t - p.t_open)/p.tau_sw)更靠谱。

3.3 求解器调用：为什么默认走两段式

单段仿真是把 0 到t_end整个 260 微秒丢给求解器，让它在里面同时处理泵浦积累和纳秒脉冲。这样做物理上没问题，工程上很痛苦：要分辨纳秒级的脉冲，最大步长必须压到 10⁻¹¹ 量级，260 微秒就是两千六百万步，跑一次十几分钟，参数扫描基本没戏。

两段式把这个坑绕开了。Q 开关关断期间，腔内光子数可以认为约等于 0，受激项全为零，反转粒子数方程退化成一条纯泵浦的解析解：

N2_open = p.Rp * p.tau * (1 - exp(-p.t_open / p.tau)); % 开门瞬间的反转粒子数

只有开门瞬间，才对 [t_open, t_open + win] 这个两百纳秒的窗口做数值积分：

y0 = [N2_open; 1]; % 初值：解析反转 + 1 个光子/m^3 的种子 opts = odeset('RelTol', 1e-8, ... 'AbsTol', [1e18, 1e-6], ... % 两个状态量量级差很多，分开设 'MaxStep', 2e-11); % 约 1/50 个往返时间，够分辨波形 [t, y] = ode15s(@(t,y) rate_eq(t,y,p), [p.t_open, p.t_open+p.win], y0, opts);

三个参数值得单独说。RelTol给 1e-8 是经验值，再松脉冲峰值会抖，再紧收益不大。AbsTol必须写成两元素向量：N₂ 的量级在 10²³ 上下，φ 的峰值可能到 10²⁰ 也可能只有 10¹⁵，用同一个绝对容差会出事——给大了光子数早早被判成零，脉冲起不来；给小了求解器在脉冲尾部死磕，步长越走越小。MaxStep卡的是分辨率，一个往返时间至少要采 20 个点，脉宽才能算准。

这里有个反直觉的地方：两段式的精度并不比单段差多少，因为泵浦段的解析解是精确的，唯一被切掉的是开门前后的那一小段过渡。真实的 Q 开关也不是瞬间打开，这个切割反而比强行数值积分更接近实际。

3.4 后处理：光子密度怎么变成示波器上的功率

求解器吐出来的是光子数密度，单位 m⁻³，没人看得懂。要转成输出功率，得算清楚每秒有多少光子从输出镜跑掉：

function [P, E, Ppk, tpk, dt] = post_process(t, y, p) phi = max(y(:,2), 0); % 数值误差可能给一点点负值，先夹住 h_nu = 6.626e-34 * 3e8 / 1064e-9; % 1064 nm 光子能量，约 1.87e-19 J % 输出镜每秒放走的光子份额 = T / t_r，乘上门限内的总光子数就是功率 P = h_nu * phi * p.V * (p.T / p.t_r); E = trapz(t, P); % 单脉冲能量，对功率积分 [Ppk, i_pk] = max(P); tpk = t(i_pk); % 半高全宽：先找跨过半高的两个样本点，再线性插值 half = Ppk / 2; i1 = find(P >= half, 1, 'first'); i2 = find(P >= half, 1, 'last'); t1 = t(i1-1) + (half - P(i1-1)) * (t(i1) - t(i1-1)) / (P(i1) - P(i1-1)); t2 = t(i2) + (half - P(i2)) * (t(i2+1) - t(i2)) / (P(i2+1) - P(i2)); dt = t2 - t1; end

p.V是模式体积，T/t_r是输出镜的等效泄漏率，两者相乘再加光子能量，就是功率。这个换算里隐含一个假设：腔内的光子密度在模式体积内是均匀的，跟 2.1 节的假设一致，不要在这里又用另一个体积，否则能量对不上。

脉宽的插值容易翻车。用interp1直接插会要求自变量单调，而 P 在半高附近可能有几个点的抖动，导致插值点选错。老老实实取相邻两点做线性插值，代码长一点但是稳。

后处理里还能顺手做一个自检：把仿真得到的单脉冲能量，跟 Degnan 给的那条解析式对一下：

E ≈ (A·hν / (2·σ·γ)) · ln(1/R) · ln(N_i / N_f)

N_i 是开门瞬间的反转粒子数密度，N_f 是脉冲结束后的残余。两者对得上（误差 15% 以内）说明方程、单位、体积三项都没写错。这个交叉验证后面还会用到。

4. 跑通流程：从零到第一张调 Q 波形

4.1 环境准备与首跑

需要的东西很干净：MATLAB R2016b 或更新版本，不需要 Simulink，也不需要任何附加工具箱。把整个q_switch_laser/目录连同scan/子目录一起放到工作路径下，cd进去，命令行敲main。第一次跑建议保持所有默认参数不动，先确认能出图。

main.m的骨架控制两件事：走哪条路径，以及出什么图。

clear; clc; close all; p = params_default(); if p.two_stage % 路径 A：泵浦段解析 + 开门段数值，速度快，推荐 N2_open = p.Rp * p.tau * (1 - exp(-p.t_open / p.tau)); y0 = [N2_open; 1]; opts = odeset('RelTol',1e-8, 'AbsTol',[1e18 1e-6], 'MaxStep',2e-11); [t, y] = ode15s(@(t,y) rate_eq(t,y,p), ... [p.t_open, p.t_open + p.win], y0, opts); else % 路径 B：单段全时域积分，慢，但能看泵浦段的完整过程 y0 = [0; 1]; opts = odeset('RelTol',1e-6, 'AbsTol',[1e18 1e-6], 'MaxStep',1e-10); [t, y] = ode15s(@(t,y) rate_eq(t,y,p), [0, p.t_end], y0, opts); end [P, E, Ppk, tpk, dt] = post_process(t, y, p); fprintf('开门时刻反转粒子数 %.3e m^-3\n', y(1,1)); fprintf('峰值功率 %.2f kW\n', Ppk/1e3); fprintf('单脉冲能量 %.3f mJ\n', E*1e3); fprintf('脉宽 FWHM %.2f ns\n', dt*1e9); fprintf('峰值时刻延迟 %.2f ns\n', (tpk - p.t_open)*1e9); figure('Position',[100 100 900 420]); subplot(1,2,1); plot((t - p.t_open)*1e9, P/1e3, 'LineWidth', 1.4); xlabel('相对开门时刻 (ns)'); ylabel('输出功率 (kW)'); title('调Q脉冲波形'); grid on; subplot(1,2,2); plot((t - p.t_open)*1e9, y(:,1)/1e23, 'LineWidth', 1.4); hold on; plot((t - p.t_open)*1e9, y(:,2)/max(y(:,2)), 'LineWidth', 1.2); xlabel('相对开门时刻 (ns)'); ylabel('归一化幅度'); legend('N_2 (10^{23} m^{-3})','\phi (归一化)'); title('反转粒子数与光子数'); grid on;

跑通之后控制台应该打出类似这样的东西：开门时刻反转粒子数 6.7e23 m⁻³，峰值功率几百 kW 量级，单脉冲能量十几 mJ，脉宽十几纳秒。数值不用跟这里完全一致，因为参数是"某台中小能量脉冲激光器"的量级，量级对就行。

一张波形图和一个双曲线图能看出很多信息：光子数曲线应该比反转粒子数曲线晚起、陡降；反转粒子数在脉冲结束时应跌到阈值以下，之后基本平掉。如果两条曲线都是平的，说明没起振，去 5.3 节。

4.2 参数逐项拆解：哪些必须改、哪些可以先放

默认参数是为了"跑得动"而不是"对得上你手上那台机器"，真要用起来，下面这几个必须换成实测值：

sigma和tau一般查得到，Nd:YAG 在 1064 nm 上是 2.8×10⁻¹⁹ cm² 和 230 μs 量级。loss_int是唯一一个经常说不清的：镜片吸收、介质散射、衍射损耗全混在一起，实验上测的往返损耗常常是理论值的两三倍。我一般先按 0.02 给，跑出阈值泵浦功率，再跟实测的阈值泵浦功率反推一个值回来，这一步比什么都管用。

T值得单独扫一遍。调 Q 激光器的输出耦合有一个最优值，太小则腔内功率密度高、损耗占比大，能量出不来；太大则阈值抬高、开门前存的反转不够。这个最优点通常在 10% 到 30% 之间，取决于泵浦强度，扫描脚本几秒钟就能画出来。

4.3 结果判读：三条线确认仿真没跑歪

波形出来之后，别急着信。三个检查挨个过一遍。

第一，脉冲峰值时刻的延迟。从开门到峰值，物理上需要经过若干次往返把光子数从种子放大到饱和，这个延迟通常是脉宽的 1 到 3 倍，也就是几十纳秒量级。如果算出来是几百皮秒，多半是增益系数写大了或者初值给太高；如果算出来是微秒，增益不够。

第二，反转粒子数的跌落深度。开门前的 N₂ 和脉冲结束后的 N₂ 之比，应该跟初始反转比阈值高多少直接相关。经验上，泵浦超过阈值 2 到 4 倍时，脉冲结束后的残余反转应该接近阈值附近，跌得不深说明能量没抽干净，跌过头（跌到零甚至负数）说明数值出问题了。

第三，能量守恒的粗账。存进去的储能是 N₂·V·hν，出来的脉冲能量除以储能，效率通常在 30% 到 70% 之间（剩下的被自发辐射和寄生损耗吃掉）。如果算出来效率 200%，一定是体积或者单位搞错了。

三条都过，再去看参数扫描的曲线。这时候曲线才有意义。

5. 避坑与排查：调 Q 仿真最容易翻车的几处

5.1 求解器报错或光子数跑成负数

现象：ode15s返回Unable to meet integration tolerances，或者跑完了但 φ 里有负值，画出来波形在零轴下面挂一条。

原因：单段模式下，从 0 到 260 微秒全程积分，开门瞬间的刚性太强，求解器为了压制误差把步长压到 10⁻¹⁵ 以下，超过了最小步长限制。负值则是绝对容差给得太松，脉冲尾部 φ 已经掉到 10⁵ 量级，还被允许有同样的绝对误差，数值噪声就把符号翻过去了。

解决：第一步，改用两段式，让泵浦段走解析解，只对开门后的两百纳秒做数值积分，这个改动通常能解决八成问题。第二步，AbsTol拆成两元素向量，给 φ 单独设一个小值。第三步，实在要用单段，把MaxStep显式卡在往返时间的 1/50，并且把NonNegative选项打开（R2018b 之后支持）。最后那句max(y(:,2),0)只是给画图用的补丁，不要拿它掩盖求解器的问题。

5.2 单位混用：cm 和 m 打架

现象：算出来的脉宽是几百微秒，或者能量只有几纳焦，跟预期差三四个数量级，但方程看上去没写错。

原因：受激发射截面教科书上普遍给 cm²（2.8×10⁻¹⁹ cm²），粒子数密度也常按 cm⁻³ 写，而长度、时间又是 SI。混着用的时候，增益系数会差 10⁴ 到 10⁶ 倍，阈值跟着全错。

解决：进params_default.m的第一行就写死注释"全部 SI"，所有查来的 cm 单位数字当场换算再填进去：截面乘 10⁻⁴，密度乘 10⁶，长度乘 10⁻²。换算完做一次量纲自检——把sigma*N2*l单独算出来，这个无量纲增益系数应该在 0.01 到 0.5 之间，超出这个范围一定是单位错了。我见过最隐蔽的一次是面积用 cm² 填、体积又用 m³ 算，单看两行都没问题。

5.3 参数扫不出脉冲：阈值和开关时刻对不上

现象：改完参数，波形是一条贴着零轴的直线，或者只有一个极小的尖峰，能量比预期低两个数量级。

原因：三种可能。一是开门时刻反转粒子数没到阈值，N2_open < N_th = loss/(2·sigma·l)，开门了也不起振。二是关断损耗loss_hold给得太小，关断期阈值不够高，泵浦还没结束就起振了，能量提前泄掉。三是泵浦速率和重复频率不匹配——重频越高，留给泵浦的时间越短，N₂ 存不满，高重频下能量下降是物理规律而不是 bug。

解决：先做一次阈值体检，把两个阈值都算出来打印：

N_th_open = (p.T + p.loss_int) / (2*p.sigma*p.l); N_th_hold = (p.T + p.loss_int + p.loss_hold) / (2*p.sigma*p.l); N2_open = p.Rp * p.tau * (1 - exp(-p.t_open/p.tau)); fprintf('开门阈值 %.3e，关断阈值 %.3e，实际反转 %.3e\n', ... N_th_open, N_th_hold, N2_open);

三条判据：N2_open要明显大于N_th_open（至少 2 倍），否则能量上不去；N2_open要小于N_th_hold，否则关断期就漏了；loss_hold通常是T + loss_int的 3 到 5 倍，按这个比例反推。做重复频率扫描时，t_open按0.8/f设，留 20% 时间给脉冲泄放和冷却，否则上一发还没冷下来就开始下一发。

5.4 扫描跑得慢、内存爆

现象：扫 50 组重复频率，跑了四十分钟还没完；或者循环跑到一半 MATLAB 提示内存不足。

原因：默认两段式的窗口给了 200 纳秒，但高重频、低能量时脉冲只有几纳秒，剩下 190 纳秒全是无效积分，白白浪费。另外每组结果都存了完整的t和y，50 组下来就是几十万个点，加上扫描维度很容易过 GB。

解决：窗口设成自适应的，先跑一发看脉宽，再按8~10 倍脉宽重设p.win，能省一个数量级的计算量。存储只留标量统计量，波形按需单独跑：

results = zeros(numel(f_list), 5); for k = 1:numel(f_list) pk = p; pk.t_open = 0.8 / f_list(k); pk.win = 10 * p.win; % 先用宽窗口探一次 % ... 求解 ... results(k,:) = [f_list(k), E, Ppk, dt, y(1,1)]; end

想要更快，就把for换成parfor，每个 worker 各自持有一份params结构体，没有共享状态，并行起来很干净。这一步在 6 核机器上大概能拿 4 倍加速，代价是启动并行池的十几秒。

6. 二次开发：把单点仿真改成能回答问题的工具

跑通一次波形只是起点，真正有价值的是把它变成一个能回答"如果……会怎样"的工具。下面几条改造路线，按投入产出比排。

最值得先做的是可饱和吸收体，因为它把模型从"主动"扩到"被动"，适用面一下宽了很多。改动很小，加一个基态粒子数密度 Ng 的方程：

% y(3) = Ng 吸收体基态粒子数密度 m^-3 Ng = y(3); % 吸收体漂白：光子越强，基态被抽空得越快 dNg = (p.Ng0 - Ng)/p.tau_g - p.sigma_g * 3e8 * Ng * phi; % 腔内总损耗 = 固定损耗 + 吸收体的可饱和损耗 loss = p.T + p.loss_int + 2 * p.sigma_g * Ng * p.lg;

p.Ng0是初始掺杂浓度，p.sigma_g是基态吸收截面，p.lg是吸收体厚度。这三个参数一起决定初始透过率T0 = exp(-2·sigma_g·Ng0·lg)，实验上直接能测，所以标定起来比主动调 Q 的关断损耗容易得多。要注意的是被动调 Q 的脉冲重复频率不是外部给的，而是由泵浦和漂白动力学自己决定的，所以两段式那套"开门时刻解析"就不适用了，得回到单段积分，步长和容差都要重新调一遍。

验证方法我一般用两条腿走路。一条是数值自洽：把时间步长和容差各收紧一倍，跑出来的峰值功率变化应该在 1% 以内，变化超过 5% 说明还没收敛。另一条是跟解析式对：用 3.4 节那条 Degnan 公式算单脉冲能量，跟数值积分结果比，误差 15% 以内算过。这两条都过不了的时候，问题基本都在单位或者模式体积上，回去查 5.2 节。

最后一个习惯。每次改完参数准备拿去跟实验对比之前，我都会强制把N_th_open、N_th_hold、N2_open这三个数打出来看一眼——它们能在三十秒内告诉你这组参数是不是物理上根本不可能起振，省下的是一整晚的无效仿真。这套 matlab 固体激光器仿真文件本身不大，真正花时间的是把参数调成你手上那台机器的样子，而这部分只能靠一次次对波形。希望帮到你。

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