MATLAB激光器谐振腔模拟：ABCD矩阵、稳区与Fox-Li迭代

简介面向激光物理、光电信息专业学生及科研人员这份MATLAB资源围绕激光器谐振腔的模拟分析展开覆盖从物理建模、参数设置到数值求解的完整流程适合课程设计、科研入门或项目预研时快速上手。压缩包共5个文件包含4个MATLAB脚本和1个TXT说明文档整体仅9KB轻量而聚焦便于逐行阅读和调试。脚本分别实现谐振腔模式函数计算、高功率激光器稳定区分析、腔内多模传输求解以及模式尺寸随腔长变化等功能说明文档则对平行平面腔自再现模式的模拟要点进行了梳理可与脚本配合理解理论到代码的映射。已有354人学习下载。通过运行代码读者可直观获得光场分布、稳定区边界等关键结果并掌握基于传输矩阵法的数值迭代思路为后续设计高功率激光器、优化谐振腔结构提供实用参考。1. 激光器谐振腔模拟为什么MATLAB比光学追迹软件更适合练手拿到一套谐振腔仿真需求多数人第一反应是打开现成光学软件把腔型摆进去等结果。这套流程确实高效但整套分析里最关键的稳定性边界、束腰位置和模式损耗全都藏在操作界面背后一旦腔型超出常见范围计算结果对不上实测你根本不知道是该改参数还是改软件设置。用MATLAB从矩阵开始写谐振腔模拟等于把黑匣子拆开重装一遍g参数稳区、往返矩阵、Fox-Li迭代、q参数传输每一步都落到数组和公式上。这篇笔记面向激光器设计和光学实验方向的从业者用一套完整的平凹腔模拟案例讲清楚谐振腔仿真怎么做、参数怎么设、哪些位置最容易翻车。2. 把谐振腔写成ABCD矩阵g参数稳定判据与可执行代码2.1 谐振腔的元器件矩阵分解谐振腔仿真的起点不是渲染腔体结构而是把每一个光学元件表达成2×2的ABCD矩阵。光线在腔内往返一周依次经过自由空间传播、镜面反射、再次自由空间传播整个过程就是一个矩阵链。对平凹腔这种最常见腔型三个基本矩阵就够了。自由空间传播距离L矩阵形式为function M fresnel_prop(L) M [1, L; 0, 1]; end曲率半径为R的反射镜近轴矩阵为function M mirror_reflect(R) M [1, 0; -2/R, 1]; end注意两个细节。第一个是反射镜矩阵的写法对凹面镜和凸面镜完全不一样R为正代表凹面镜R为负代表凸面镜符号写反后续所有稳定性和光斑结果全部颠倒。第二个是平凹腔的平面镜R取Inf在MATLAB里直接传Inf计算即可-2/Inf得到0矩阵退化为单位阵。把三个矩阵按光线传播顺序相乘就得到从平面镜出发、经凹面镜反射、再回到平面镜的往返矩阵R2 5.0; % 凹面镜曲率半径单位m L 1.0; % 腔长单位m M_round fresnel_prop(L) * mirror_reflect(R2) * fresnel_prop(L);矩阵乘法顺序代表光线作用顺序第一个作用在光线上的矩阵放在最右边。很多初学者习惯从左往右读结果把传播方向搞反了算出来的光斑尺寸与实际差出几倍。我一般会在脚本开头用注释把顺序写清楚“从平面镜出射→传播L→凹面镜反射→再传播L”。这个往返矩阵是整个谐振腔数值分析的骨架。后续无论是算g参数、用q参数求束腰还是做Fox-Li迭代都要从它出发。矩阵写对后边所有计算才有意义。2.2 g参数稳定性条件与稳区图有了往返矩阵稳定性判据可以直接从矩阵元素读出也可以回到经典g参数定义。对腔长L、两镜曲率半径R1和R2g参数定义为g1 1 - L/R1 g2 1 - L/R2稳定腔条件为0 g1×g2 1等号对应临界腔——共焦腔或平行平面腔正落在边界上实际工程中极少把工作点设计在边界附近因为微小装配误差就会让g乘积越界模式损耗急剧增大。在MATLAB里画稳区图直接生成g1-g2平面的等值线g1 linspace(-2, 2, 401); g2 linspace(-2, 2, 401); [G1, G2] meshgrid(g1, g2); prod G1 .* G2; stable (prod 0) (prod 1); imagesc(g1, g2, stable); axis xy; colormap(gray); xlabel(g1 1 - L/R1); ylabel(g2 1 - L/R2); title(稳定腔区域图);linspace(-2, 2, 401)把g参数范围设得比实际腔型域更大方便观察稳区边界meshgrid生成二维网格后prod 0 prod 1是一次性做双重条件判断得到的就是稳定腔区域掩模。这里用imagesc而不是contour是因为稳定/非稳定是二值判断用图像显示边界最直观。矩形硬边界看起来不光滑把网格点加密到401×401后边界位置的视觉误差已经小于1%。画完图再把具体腔型的g工作点叠加进去R1 Inf; % 平面镜 g1_val 1 - L/R1; % g1 1 g2_val 1 - L/R2; % g2 1 - 1/5 0.8 hold on; plot(g1_val, g2_val, ro, MarkerSize, 8, LineWidth, 2);这组参数计算出来g1×g20.8落在稳定区域内部且离边界有足够余量是典型的工作点设计。平凹腔里g1恒为1所以工作点是一条竖线腔长越接近曲率半径R2g2越接近0工作点越靠近稳区边界。参数选型上需要解释一下为什么选择R25m、L1m而不是其他值因为这个腔型g1g2离边界还有20%余量镜片加工容差和热畸变导致的曲率变化不至于让腔直接进入非稳区。实际设计时我建议至少留15%的稳区余量否则激光器装调时会反复出现不出光的情况每次都要怀疑镜子装歪了其实根源是参数选得太贴边界。2.3 非稳腔为何值得留意稳定腔条件解决的是低阶模低损耗问题但很多高功率激光器用的是非稳腔。非稳腔的g1g2乘积落在0到1区间之外基模损耗大却有很好的横模鉴别能力输出耦合可以做得非常高适合大增益体积、大模体积的振荡器设计。正向非稳腔的几何放大率由往返矩阵的迹决定m (abs(M_round(1,1) M_round(2,2)) sqrt((M_round(1,1) M_round(2,2))^2 - 4)) / 2;放大率m可以粗略理解为每一往返后光斑横向扩大的倍数。这个参数对非稳腔设计很关键——输出耦合效率和镜面尺寸都要按放大率推算。不过在入门阶段先集中把稳定腔的模拟做扎实非稳腔的衍射迭代计算需要更谨慎的网格采样直接套用稳定腔的代码会出很多意想不到的数值问题。3. Fox-Li迭代求解腔模衍射积分离散化与收敛判据3.1 自再现模原理与角谱传播g参数只能回答“这个腔稳不稳定”回答不了“腔内光斑长什么样”。要得到基模横向分布和衍射损耗必须做本征模式求解。经典Fox-Li迭代的核心是自再现思想光在腔内往返足够多次后横向场分布趋于稳定每一往返只改变一个复常数增益因子分布形状不再变化。实现这一思想通常借助角谱衍射传播。把镜面M1上的场分布做二维傅里叶变换在频域乘以近轴传播传递函数再逆变换就得到传播L后在镜面M2上的场分布。单次往返过程包含两次自由空间传播和两次镜面反射。3.2 谐振腔单程传播的MATLAB实现仿真参数和网格先定义好lambda 1064e-9; % Nd:YAG激光波长单位m L 1.0; % 腔长单位m R2 5.0; % 凹面镜曲率半径单位m D 6e-3; % 凹面镜有效孔径直径单位m N 256; % 采样网格点数方形网格 dx 0.2e-3; % 空间采样步长单位m x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); r2 X.^2 Y.^2; % 径向坐标平方单位m^2 % 频域坐标 fx (-N/2 : N/2-1) / (N * dx); [FX, FY] meshgrid(fx, fx); f2 FX.^2 FY.^2; % 凹面镜反射相位因子薄透镜近似 phase_mirror exp(-1i * pi / (lambda * R2) * r2); % 凹面镜孔径遮挡 aperture zeros(N, N); aperture(r2 (D/2)^2) 1;这里有几个参数必须提前想清楚。N256是采样点数点数太少频谱分辨率低相位调制细节丢失dx0.2mm决定空间网格步长窗口总宽度是256×0.2mm51.2mm大约是孔径直径6mm的8.5倍既能完整保留边缘衍射又不会因为窗口过大浪费采样点。频域坐标fx的范围由1/(N*dx)决定最大空间频率约1/0.2mm5000m⁻¹对应衍射角约0.005弧度对5m曲率半径镜面产生的相位变化来说完全够用。初始场选择高斯分布w0 1.5e-3; % 初始光斑半径猜测值单位m u exp(-r2 / w0^2); u u / sqrt(sum(abs(u(:)).^2) * dx^2); % 能量归一化初始光斑半径的猜测值允许有误差迭代会逐渐收敛到自再现模。需要说明的是如果初始光斑给得和真实基模偏差太大比如给了平面波迭代次数会明显增加甚至需要几千次才收敛如果给得太窄高频分量太多可能激发高阶模。用高斯分布做初值是最稳妥的常见做法。接下来是迭代主循环max_iter 2000; tol 1e-6; last_field u; loss_history zeros(max_iter, 1); for iter 1:max_iter % 从平面镜传播到凹面镜角谱法 U fftshift(fft2(u)); k 2 * pi / lambda; H exp(1i * k * L) .* exp(-1i * pi * lambda * L * f2); u2 ifft2(ifftshift(U .* H)); % 凹面镜反射施加曲率相位与孔径遮挡 u2 u2 .* phase_mirror; u2 u2 .* aperture; % 从凹面镜传播回平面镜同样用角谱法 U2 fftshift(fft2(u2)); u1 ifft2(ifftshift(U2 .* H)); % 归一化记录能量损耗 energy_before sum(abs(u1(:)).^2) * dx^2; u1 u1 / sqrt(energy_before); loss_per_round 1 - energy_before; % 收敛判断归一化场分布的变化量 diff sum(abs(u1(:) - last_field(:)).^2) * dx^2; last_field u1; loss_history(iter) loss_per_round; u u1; if diff tol fprintf(迭代收敛于第%d次\n, iter); break; end end这段代码的逻辑分四步正向传播、凹面镜作用、反向传播、归一化与收敛判断。正向传播和反向传播用的是同一个传递函数H因为两次传播距离都等于腔长L腔结构对称时可以直接复用。凹面镜作用分两步实施先乘相位因子反映曲率对波前的弯折再乘孔径掩膜反映镜面有限尺寸造成的衍射损耗两者物理意义不同不能合并成一个矩阵。fftshift和ifftshift的配对容易写错——傅里叶变换前要把场分布从左下角排列挪到中心排列逆变换后再挪回来。如果只挪一次不配对结果会多出一个整体空间偏移而且这个偏移很难靠肉眼从强度图上发现通常到和解析光斑半径做对比时才暴露。3.3 收敛后的模式损耗分析与物理解读迭代结束后我们从两个维度看结果横向强度分布与衍射损耗。intensity abs(u).^2; intensity intensity / max(intensity(:)); % 沿x轴截取分布拟合光斑半径 profile abs(u(N/21, :)).^2; w_fitted dx * sqrt(2 * sum(profile .* (x - 0).^2) / sum(profile)); fprintf(Fox-Li迭代基模光斑半径%.2f mm\n, w_fitted*1e3); % 稳态损耗 steady_loss mean(loss_history(end-100:end)); fprintf(稳态单程损耗%.4f%%\n, steady_loss*100);光斑半径用二阶矩定义计算而不是直接取1/e²点因为数值解分布不严格是高斯形二阶矩能更客观地反映能量扩展。这里有个很反直觉的现象归一化后的场分布形状不再变化但每一往返能量仍然有小幅损耗这个损耗来自镜面孔径的硬边衍射只要孔径有限损耗就存在。若把孔径去掉迭代能量基本守恒场分布也会慢慢展宽且收敛不到稳定形态。硬边衍射是谐振腔模拟中少数“必须有”的耗散机制工程上用来抑制高阶模数值上则保证自再现方程有非零解。4. 高斯光束q参数传输束腰位置与光斑尺寸的快速计算4.1 q参数、光斑尺寸与波前曲率的关系Fox-Li迭代给出数值解但实际工程里更常用复光束参数q做解析计算。q参数和光斑半径w、波前曲率半径R的关系写成1/q 1/R - i·λ/(π·w²)实部对应波前曲率虚部对应光斑尺寸。这个表达方式的工程价值在于q参数经ABCD矩阵传输后用双线性变换更新q (A·q B) / (C·q D)只要知道入射光束在某个参考面的q参数任意位置的光斑尺寸、曲率半径都能用复数运算算出来不需要逐点做衍射积分。数值上它比Fox-Li快几个量级适合做参数扫描和实时反馈。4.2 平凹腔自再现q参数求解稳定腔内往返一周后q参数必须自再现即q经过往返矩阵后回到原值。这给出一元二次方程程序里直接求根R2 5.0; L 1.0; M_round fresnel_prop(L) * mirror_reflect(R2) * fresnel_prop(L); A M_round(1,1); B M_round(1,2); C M_round(2,1); D M_round(2,2); % 自再现方程 C*q^2 (D-A)*q - B 0 p [C, D-A, -B]; roots_q roots(p); % 筛选虚部为负的物理有效根 valid roots_q(imag(roots_q) 0); q0 valid(1); w0 sqrt(-lambda / (pi * imag(1/q0))); R_curv 1 / real(1/q0); fprintf(基模束腰半径 w0 %.2f mm\n, w0*1e3);这段代码的关键在roots函数两项选择一元二次方程有两个根只有虚部为负的根对应物理光束。若不加筛选直接取第一个根可能拿到虚部为正的发散解光斑半径变成复数后边所有物理量全是错的。这是一个非常典型的翻车点。自再现束腰半径计算出来后束腰位置也可以推算。平凹腔中凹面镜的曲率中心附近会形成束腰严格位置从平面镜出发找q虚部为零的传输距离z linspace(0, L, 1001); wz_list zeros(size(z)); for i 1:length(z) Mz fresnel_prop(z(i)); qz (Mz(1,1)*q0 Mz(1,2)) / (Mz(2,1)*q0 Mz(2,2)); wz_list(i) sqrt(-lambda / (pi * imag(1/qz))); end [~, idx_min] min(wz_list); z_waist z(idx_min); fprintf(束腰距平面镜距离%.3f m\n, z_waist);linspace(0, L, 1001)把平面镜到凹面镜的整个腔长范围细分成1001个截面逐点计算q参数再提取光斑半径。[~, idx_min] min(wz_list)这行返回光斑最细处的位置。需要注意这个束腰位置搜索依赖于先前解出的q0如果第2.2节往返矩阵顺序写反q0本身就是镜像解扫描出的束腰位置也会落在错误的端面。4.3 q参数解析与Fox-Li迭代结果的交叉验证解析公式和高阶数值迭代不应该彼此孤立。我在实际项目中习惯同时跑两组仿真对比结果判断参数设得对不对对比项q参数解析法Fox-Li迭代法计算开销毫秒级数毫秒到数秒物理基础近轴ABCD矩阵衍射积分自再现适用腔型任意稳定/非稳定腔含孔径、硬边效果时输出维度光斑半径、曲率半径完整横向场分布误差来源忽略衍射损耗网格离散化误差q参数法得到的光斑半径是“无限大镜面”下的理想基模Fox-Li迭代在加了6mm孔径后得到的半径通常会小一点因为硬边衍射限制了光场扩展。如果两者差距在3%以内说明网格和采样参数设置合理数值结果可信差距超过10%优先检查网格步长和窗口宽度再检查初始场是否收敛到了高阶模。这个交叉验证习惯花费的时间不到五分钟却能把Fox-Li迭代中绝大多数肉眼察觉不到的数值瑕疵暴露出来。我每次换新腔型都会强制跑一遍对比很大程度上避免了后边实验调试时把“仿真误差”和“装配误差”搞混。5. 谐振腔仿真避坑采样窗口、初始场与硬边衍射的典型翻车点5.1 采样窗口过窄导致模式能量被截断现象Fox-Li迭代收敛后光斑分布边缘出现明显的矩形条纹强度分布不光滑拟合的半径偏小且和q参数解析结果差距越来越大。原因网格总宽度N×dx小于光斑实际扩展范围场分布碰到计算窗口边缘被强制截断等效于人为加了一个不该有的矩形光阑。解决先按q参数解析公式估算光斑半径再设定窗口宽度为预估光斑半径的6到8倍。具体操作是用当前代码跑一次在收敛后统计光斑半径w_fitted如果w_fitted大于窗口宽度的1/6把dx调大或增加N重算重新对比解析值。5.2 网格步长过大导致频域传递函数采样不足现象迭代不收敛损耗曲线振荡且幅度不下降。频率谱在高频端出现折叠特征强度图上有精细的莫尔条纹。原因角谱传播中频域相位φπλLf²随f²增长当dx过大最高空间频率分量对应的相位变化超过π频域欠采样造成相位混叠。解决检查频域坐标最大值的相位是否满足πλLf_max² π/2。以λ1064nm、L1m为例f_max必须小于约3.86×10⁴ m⁻¹对应dx必须大于1/(2·f_max·N)约0.65μm——实际工程中取dx0.2mm已经远大于这个下限问题更多出现在dx设置过小导致f_max过大。这里容易搞反dx过小同样有害。我踩过的坑是优先把dx往小调结果N固定时窗口变小反而触发5.1的问题。5.3 初始场给随机噪声导致收敛缓慢甚至落入高阶模现象相同腔型参数换成随机相位初始场后迭代两千次还没收敛或者收敛后光斑分布带有一个明显的相位旋转损耗值明显偏高。原因随机相位场包含大量高频分量和多个横模成分Fox-Li迭代收敛的是损耗最低的横模但随机初值不容易把能量集中到该模式上迭代过程在模式竞争的“中间态”里停留很久。解决初始场用高斯分布宽度取预估光斑半径的0.8到1.5倍。就算给宽一点迭代前几十轮会快速收敛到基模邻域。想验证高阶模存在性时可故意用厄米-高斯组合做初值但要清楚此时得到的可能是多模叠加而不是纯基模不能直接当作单模结果用。5.4 忘记施加镜面孔径迭代永远不收敛现象去掉aperture后每轮损耗接近零光斑半径不断展宽且没有稳定趋势迭代差值diff始终在10⁻³量级不下来。原因无孔径谐振腔在近轴近似下无衍射损耗自再现方程没有有限尺寸的非零解。能量不断向外扩展等价于场分布永不闭合。这在物理上是“理想腔没有稳态解”在数值上就是迭代发散。解决给两个镜面中至少一个设定实际通光孔径。如果真实系统里镜架、泵浦模块本身有明显限孔可以把等效孔径直接设为机械限制中最窄的一处直径。这里给一个经验值孔径直径取预估基模光斑直径的3到5倍既保留明显的基模损耗优势又不至于把损耗压得太大导致能量利用率下降。5.5 矩阵乘法顺序写反腔内往返方向颠倒现象q参数法解出的w0和Fox-Li迭代结果差得离谱甚至算出的束腰负值或虚值。原因公式上往返矩阵写成M_round fresnel_prop(L) * mirror_reflect(R2) * fresnel_prop(L)但如果习惯性地把矩阵沿着“平面镜→凹面镜→平面镜”的物理顺序从左往右排就变成了T(L) * R(R2) * T(L)的自右向左作用方向反了等价于交换了R1和R2对非对称腔结果完全不同。解决给矩阵链加注释标明靠近向量的矩阵是第一个作用在光线上的元件。验证方法很简单把R1和R2交换位置看两个镜面分别是无穷大平镜和5m凹镜时的结果是否交换。如果交换后结果一样说明本来就结构对称检验失效对平凹腔这种不对称腔交换后w0应当显著变化若没变化说明矩阵顺序写错。6. 从稳区图到参数扫描验证仿真结果的三个进阶习惯习惯一做参数扫描前先把单点自洽检查跑通。我见过不少人在参数扫描里得到一条光滑曲线却忘了验证其中一个点是否物理正确。正确流程是先固定一组工作点用q参数法算光斑半径和束腰位置再让Fox-Li迭代收敛比较两个结果差异在3%以内然后才开始扫描。习惯二腔长扫描时注意稳区边界突变。以R25m、R1Inf的平凹腔为例腔长L从0.5m扫到4.9mg2从0.9单调降到0.02g1g2乘积从0.9降到0.02束腰半径在L接近5m时发散。用一个循环把这段趋势画出来L_list linspace(0.5, 4.9, 200); w0_list zeros(size(L_list)); for i 1:length(L_list) M fresnel_prop(L_list(i)) * mirror_reflect(R2) * fresnel_prop(L_list(i)); A M(1,1); B M(1,2); C M(2,1); D M(2,2); rts roots([C, D-A, -B]); rts rts(imag(rts) 0); if isempty(rts) w0_list(i) NaN; else w0_list(i) sqrt(-lambda / (pi * imag(1/rts(1)))); end end plot(L_list, w0_list*1e3); xlabel(腔长 L/m); ylabel(束腰半径 w0/mm);这段扫出来的曲线会在L接近R2时急剧抬升曲线尽头w0发散。看到这种现象不要慌这不是程序错误而是谐振腔物理上在临界点失效的体现。用这个图选工作点我会把L设在1m到4m之间的中段区域这个范围既远离稳区边界又有相对小的光斑工程调试容差好。习惯三双参数稳区扫描确认设计点不是孤点。把g1和g2同时扫描画g1g2乘积的等高线看目标工作点周围是否存在连续的稳定区域。只做单参数扫描容易漏掉二维空间里的窄带通道。双参数扫描用contourf实现稳定区内填充绿色不稳定区留白目标工作点用星号标记一眼就能看出余量方向。从那以后我每次搭建新的谐振腔模型都强制走一遍这套流程g参数验稳区、q参数定束腰、Fox-Li验证横向分布、扫描确认设计余量单点合格再做扫描。这套组合拳帮我躲过了很多数不清的翻车现场尤其是那些仿真图上看起来挺漂亮、实际装调却根本出不了光的腔型。希望帮到你。本文还有配套的精品资源点击获取