尧图网站建设 
从研究生一年级第一次在实验室里看到空间光调制器把一束普通高斯光“拧”成一个甜甜圈形状开始我就对Hermite-GaussianHG和Laguerre-GaussianLG这两种光束模式产生了兴趣。后来做光束整形和光镊实验几乎所有需要模式结构的地方都离不开这两类解。说实话纯解析推导容易把人劝退但Matlab模拟可以很直观地把公式变成眼前跳动的光斑图案。这篇文章就把我自己常用的模拟思路、代码细节和踩坑记录整理出来写给正在做光学模拟、激光模式分析或相关课程设计的同学参考。先说清楚这是个什么东西Hermite-Gaussian光束是傍轴波动方程在直角坐标系下的本征解光斑呈矩形网格状排列Laguerre-Gaussian光束是同一方程在柱坐标系下的本征解带有螺旋相位、中心强度为零就是我们常说的涡旋光。两类模式在激光谐振腔设计、光镊、超分辨显微、光通信模式复用等领域都是核心角色。Matlab模拟的核心价值在于不需要搭真实光路就能研究参数变化对光束行为的影响也能直接验证实验方案里某个相位板或者空间光调制器加载的图案到底该长什么样。尤其适合光学工程、物理专业的学生和刚接触光场调控的工程师上手。下面我从理论公式讲起把代码实现、结果分析和问题排查都过一遍。整个思路是一边解释“为什么这么写”一边给出可以直接抄作业的Matlab代码。1. 理论公式与参数设计思路1.1 HG光束的数学表达与物理直觉Hermite-Gaussian光束在z平面上的复振幅可以写成[ E_{nm}(x,y,z)A\frac{w_0}{w(z)}H_n\left(\frac{\sqrt{2}x}{w(z)}\right)H_m\left(\frac{\sqrt{2}y}{w(z)}\right)\exp\left(-\frac{x^2y^2}{w^2(z)}\right)\exp\left(ik\frac{x^2y^2}{2R(z)}\right)\exp\left(-i(nm1)\psi(z)\right) ]初学者看到这个公式往往直接头皮发麻我拆开给你看。(H_n)和(H_m)是n阶和m阶厄米多项式决定了光场在x方向和y方向上的振荡结构。n0时没有横向振荡n1时有一个过零节点n2时有两个节点以此类推。所以HG光束的强度图案就像一张网格纸暗线对应于厄米多项式的零点。(w(z))是光束在传播距离z处的束半径(w_0)是束腰处的半径。它描述了光束横向尺寸的演化远离束腰时光斑会逐渐发散。(R(z))是波前曲率半径表示等相位面的弯曲程度。束腰处波前是平面R为无穷大远离束腰波前逐渐变成球面。(\psi(z)\arctan(z/z_R))是Gouy相位其中(z_R\pi w_0^2/\lambda)是瑞利距离。这个相位累加项在单纯观察横向光强分布时可以忽略但你如果要模拟干涉或者计算谐振腔模式就必须保留。直观理解HG光束可以这样想拿一块矩形的弹性膜不同振动模式就对应不同的HG阶数。00阶就是普通高斯光斑一个圆点10阶是左右两个瓣中间一条暗线01阶是上下两个瓣11阶是四瓣。阶数越高光斑里的“格子”越多能量分布的范围也越大。模拟这类光束时最关键的参数有三个波长(\lambda)、束腰半径(w_0)、模式阶数((n,m))。另外模拟的横向计算区域大小和采样点数直接决定了图案分辨率和混叠程度这个后面展开说。1.2 LG光束的数学表达与相位奇点Laguerre-Gaussian光束在柱坐标下表示为[ E_{pl}(r,\phi,z)A\frac{w_0}{w(z)}\left(\frac{\sqrt{2}r}{w(z)}\right)^{|l|}L_p^{|l|}\left(\frac{2r^2}{w^2(z)}\right)\exp\left(-\frac{r^2}{w^2(z)}\right)\exp\left(-il\phi\right)\exp\left(ik\frac{r^2}{2R(z)}\right)\exp\left(-i(2p|l|1)\psi(z)\right) ]这里出现了两个新的角色。(L_p^{|l|})是广义拉盖尔多项式p是径向指数决定了光强在径向上有几个亮环。p0时径向没有零点和暗环p1时会出现一个径向暗环光斑呈现“单环套双瓣”的结构。(l)是拓扑荷数也叫涡旋阶数它对应的相位因子(\exp(-il\phi))是LG光束的灵魂。旋转方位角(\phi)一圈相位改变(2\pi l)等相位面是一个螺旋面。中心处因为相位不确定振幅必须为零所以强度图案中心是暗的形成“甜甜圈”。如果你用过空间光调制器产生涡旋光你会发现在中心总有一个黑点那就是相位奇点。拓扑荷(l)的绝对值越大暗斑半径越大甜甜圈越“胖”。p的数值则影响径向次级结构比如p1时中心暗斑外围会再出现一个暗环。HG和LG之间不是孤立的它们都是同一个光学系统在不同坐标系下的完备基。用适当系数将一个方向的HG模式叠加可以合成LG模式。这个关系在复杂光场调控中很有用比如你想生成任意空间结构光束本质上就是在一个基底下做模式展开。Matlab模拟时理解这一点能帮你检查程序结果是否合理——比如00阶HG和00阶LG完全一样都是普通高斯光。1.3 为什么模拟而不是直接实验有人会问公式都有了直接拿激光器和空间光调制器做实验不就行了模拟的意义在哪里一个很现实的原因是成本。搭建一套完整的实验光路需要激光器、透镜组、空间光调制器、CCD相机设备动辄几十万而且每一次调节都要花大量时间。相比之下Matlab模拟只需要一台普通电脑5分钟就能把任意阶数的HG和LG光斑画出来。在项目预研阶段模拟可以快速验证想法这个光束经过透镜聚焦后光斑长什么样如果我在光路上加一个与LG相位共轭的波片能不能把涡旋光还原成高斯光这些问题在代码里修改几个参数就能得到初步答案再决定值不值得上实验。另一个原因是精度可控。实验里总存在像差、噪声、对准误差很难单独观察某一个参数变化的影响。模拟则可以精确控制(w_0)、(\lambda)、模式阶数单独剥离某个变量这是理论推导和实验之间很好的桥梁。我自己习惯把所有模拟结果保存成标准格式的图片和数据方便在写论文或者做汇报时直接引用。2. Matlab数值实现方案2.1 模拟参数的选取依据写Matlab代码之前先把模拟区域和采样参数确定下来。这里的参数选择直接决定结果是否可信我常用的配置如下波长(\lambda)取632.8 nm也就是常用氦氖激光器的波长。这个值在模拟中只影响瑞利距离和传播相位不改变横向光斑形状但后续计算衍射要用。束腰(w_0)取1 mm。这个尺寸和实验室常见的激光器光斑大小接近方便直观对照。横向计算区域L取20 mm × 20 mm边长是束腰的20倍。这个范围足够容纳到高阶模式的整个光斑也不会出现能量“溢出”边界的问题。采样点数N取512。这是一个比较好的平衡点既能看清精细结构计算速度也快。选择合适的计算区域很关键。如果L相对于束腰太小高阶模式的外围旁瓣会被截断强度图案出现明显的矩形边界伪影。如果L太大而N不够中心区域的采样稀疏光斑细节就没了。经验法则是L至少是最高阶模式整体宽度外包络的2到3倍采样间隔要小于最小结构的四分之一周期。坐标网格用meshgrid生成。注意这里我用的是物理坐标单位是米后续显示图像时再转换成毫米方便阅读。lambda 632.8e-9; % 波长单位m w0 1e-3; % 束腰半径单位m z 0; % 观察位置腰斑处 L 20e-3; % 计算区域边长单位m N 512; % 每个方向采样点数 x linspace(-L/2, L/2, N); y x; [X, Y] meshgrid(x, y); [R, Phi] cart2pol(X, Y); % 柱坐标下的径向距离和方位角这里使用cart2pol把直角坐标网格转换成柱坐标网格后面计算LG光束会大量用到径向距离(R)和方位角(\Phi)。很多人会直接写R sqrt(X.^2Y.^2)和Phi atan2(Y,X)效果一样但cart2pol可读性更好一点我习惯用它。束腰处z0时(w(z)w_0)(R(z)\to\infty)Gouy相位(\psi(z)0)公式会大幅简化。实际模拟中我最常在z0处画横向光斑因为这是最典型的观察面也最容易和实验中的CCD采集图像对照。如果要模拟传播过程再引入z循环即可。2.2 厄米多项式的高效递推实现Matlab的符号工具箱里确实有hermiteH函数可以求厄米多项式比如hermiteH(3, x)。但符号计算在512×512的网格上会非常慢而且还要求安装Symbolic Math Toolbox。我自己的项目里从来不用它做二维扫描而是用递推关系实现数值版本。厄米多项式满足递推公式[ H_0(x)1,\quad H_1(x)2x,\quad H_{n1}(x)2xH_n(x)-2nH_{n-1}(x) ]这个递推只需要几次数组乘法和加法比符号计算快几个数量级。我封装成一个函数function H hermite_poly(n, x) % Hermite多项式数值递推 % n: 多项式阶数x: 自变量矩阵 if n 0 H ones(size(x)); elseif n 1 H 2 * x; else H0 ones(size(x)); H1 2 * x; for k 1 : n-1 H2 2 * x .* H1 - 2 * k .* H0; H0 H1; H1 H2; end H H1; end end注意这段代码里2kH0中的k是循环变量代表当前递推步数。厄米多项式递推公式里的系数是2nn就是正在构造的多项式序号所以循环变量k从1跑到n-1时正好对应生成H2之前的那一步。实际用的时候这个函数返回的是和输入x同尺寸的矩阵因为x是meshgrid产生的二维数组ones(size(x))保证每一步都是矩阵运算。这里有一个容易犯的错误递推的初始值H0和H1必须是和x同样大小的矩阵不能用标量1和2*x代替否则矩阵维度对不上。有些初学同学在这里栽跟头报错说“矩阵维度必须一致”其实就是因为初始值没写成矩阵形式。2.3 广义拉盖尔多项式的递推实现广义拉盖尔多项式(L_p^\alpha(x))也有现成的递推公式这里(\alpha|l|)[ L_0^\alpha(x)1 ] [ L_1^\alpha(x)1\alpha-x ] [ L_{k1}^\alpha(x)\frac{(2k1\alpha-x)L_k^\alpha(x)-(k\alpha)L_{k-1}^\alpha(x)}{k1} ]实现的时候需要注意一个细节公式里的x是归一化变量在LG公式里取(x2r^2/w^2(z))。调用递推函数时传入的应当是整个归一化变量矩阵而不是径向距离本身。很多教程直接传r结果画出来的光斑径向结构完全不对就是这个原因。function L laguerre_poly(p, alpha, x) % 广义拉盖尔多项式递推 % p: 多项式阶数, alpha: 参数, x: 自变量矩阵 if p 0 L ones(size(x)); elseif p 1 L 1 alpha - x; else L0 ones(size(x)); L1 1 alpha - x; for k 1 : p-1 L2 ((2*k1alpha-x) .* L1 - (kalpha) .* L0) / (k1); L0 L1; L1 L2; end L L1; end end这段代码里除以(k1)是唯一一个除法操作用的右除/对矩阵来说相当于每个元素都除以同一个标量结果还是同尺寸矩阵。递推过程涉及x与L1的元素级相乘所以不要用而用.这是Matlab新手最常见的坑。如果你把.写成了Matlab会尝试做矩阵乘法矩阵尺寸不对就报错或者更奇怪的是在方阵情况下能“乘”出来但结果完全错误。我见过有人在代码里用*把整个结果算得面目全非还不自知的情况所以这里特意强调。2.4 HG光束复振幅计算与可视化有了递推函数HG光束的复振幅就是一行数学式的直接翻译。在束腰位置z0处(w(z)w0)(R(z))不考虑Gouy相位为0公式简化为[ E_{nm}(x,y)H_n(\sqrt{2}x/w_0)H_m(\sqrt{2}y/w_0)\exp\left(-\frac{x^2y^2}{w_0^2}\right) ]归一化常数可以先不写因为最后显示强度图时都要做归一化处理。数值计算的代码如下n 2; m 1; % HG模式阶数 sq2 sqrt(2); wx sq2 * X / w0; wy sq2 * Y / w0; Hx hermite_poly(n, wx); Hy hermite_poly(m, wy); G exp(-(X.^2 Y.^2) / w0^2); E_HG Hx .* Hy .* G; I_HG abs(E_HG).^2; I_HG I_HG / max(I_HG(:)); % 归一化这里有几个关键点值得展开。第一厄米多项式自变量是(\sqrt{2}x/w(z))不是单纯x。这个缩放因子来源于高斯函数指数项的自变量匹配。如果漏掉(\sqrt{2})得到的“HG光束”只是高斯包络和普通厄米多项式的乘积并不是真正的傍轴本征模图案虽然长得相似但节点位置和亮斑间隔都是错的。我在给组里新人审代码时至少见过三次这个错误每次都是光斑看起来“有点怪但说不上哪里怪”。第二强度图是(|E|^2)不是E本身。复数振幅里包含了相位信息直接画abs(E)得到的是振幅分布不是光强。虽然两者的亮暗结构类似但物理含义不同写注释时要说清楚。第三归一化用max(I_HG(:))而不是max(I_HG)。前者找出整个矩阵的最大值后者只会找出每一列的最大值返回一个行向量。如果你后续要做多幅图对比每幅都用各自全局最大值归一化才能保证亮度比例一致。显示强度图时用imagesc配上colormap。我个人比较喜欢用parula这是Matlab默认的配色色觉障碍者也能分辨比老式jet更科学。坐标轴单位习惯转成毫米显示图像标题里注明(n,m)阶数。figure; imagesc(x*1e3, y*1e3, I_HG); axis image; colormap(parula); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(HG(%d,%d) intensity, n, m));axis image保证横纵坐标刻度一致不会因为窗口拉伸让圆形光斑变成椭圆。这个细节看着小但在论文插图里非常影响观感。2.5 LG光束复振幅计算与相位图LG光束在z0处的复振幅计算同样直接。唯一要注意的是奇点处理。公式中有一个因子((\sqrt{2}r/w_0)^{|l|})当(r0)且(l\neq 0)时数值为0所以中心强度为0。但如果网格的中心正好没有采样点比如N是偶数中心位置落在四个像素的交界处数值上中心光强不一定严格为0图案上会出现一个极小值而不是完美的黑洞。解决办法是让N取奇数或者在生成网格后用ifftshift调整坐标原点。我一般在模拟中用N511或513确保原点落在网格点上。LG光束的代码p 0; l 1; % LG模式径向阶数p拓扑荷l alpha abs(l); radial_arg 2 * R.^2 / w0^2; Lp laguerre_poly(p, alpha, radial_arg); E_LG (sqrt(2)*R/w0).^alpha .* Lp .* exp(-R.^2/w0^2) .* exp(-1i*l*Phi); I_LG abs(E_LG).^2; I_LG I_LG / max(I_LG(:));相位图比较复杂。LG光束的相位分布是螺旋状的直接使用angle函数会得到从-π到π的跳变在跳变处会出现一条明显的红蓝分界线这就是相位的2π卷绕。严格来说这不是错误螺旋相位本来就有一个2π相位不连续面实验里干涉测量时也会看到这条亮线。但如果你想让相位图看起来更连续、更像一个完整的螺旋面就需要对相位做解卷绕处理。Matlab里unwrap函数只能沿一维数组解卷绕不能直接处理二维相位图。如果你有Image Processing Toolbox可以用wrapToPi和相位解缠算法如果没有工具箱可以自定义一个简单的二维解缠绕实现或者干脆接受这种锯齿状显示。我自己的做法是显示原始angle相位图但在图注里注明“相位卷绕为[-π,π]区间”这样是规范和诚实的方式。很多期刊的论文插图也这么干。另外要记得画一下LG光束的相位图因为(l±1)的相位图长得很像一片风扇叶旋转方向相反。这是判断拓扑荷符号的最直观方法。模拟时改变l的正负相位图的旋转方向会反过来强度图则几乎不变。这个特性在光学实验中经常用来标定涡旋光的拓扑荷符号。3. 多模式模拟与结果解读3.1 不同阶数HG光束的对比写一个循环就能一口气生成一组HG(0,0)到HG(3,2)的强度图。这里我用subplot把多幅图放到同一个figure里figure(Position, [100 100 1200 800]); idx 1; for n 0 : 3 for m 0 : 2 wx sq2 * X / w0; wy sq2 * Y / w0; Hx hermite_poly(n, wx); Hy hermite_poly(m, wy); Enm Hx .* Hy .* G; Inm abs(Enm).^2; Inm Inm / max(Inm(:)); subplot(4, 3, idx); imagesc(x*1e3, y*1e3, Inm); axis image off; colormap(parula); title(sprintf(HG(%d,%d), n, m)); idx idx 1; end end运行后你会注意到几个规律总阶数(Nnm)相同的HG模式比如HG(2,0)和HG(1,1)它们光斑的整体空间范围差不多因为Gouy相位和光束扩散速度都只依赖(N1)x方向和y方向的亮斑数量独立HG(2,0)在x方向有三个峰在y方向只有一个峰而HG(1,1)在两个方向上各有两个峰。还有一个容易被忽略的性质HG模式在传播过程中会保持自己的形状不变只是整体尺寸缩放和附加一个Gouy相位。这就是为什么它们被称为谐振腔的“本征模”。如果模拟中看到形状随z变化通常是公式里的归一化缩放和相位项没写对。3.2 不同拓扑荷LG光束的结构差异同样方法可以生成p0时l从-3到3的LG模式对比图。你会发现|LG|的强度图案完全相同——l1和l-1的甜甜圈一样大一样亮唯一的区别是相位沿相反方向旋转。这表明LG光束的强度分布只依赖|l|而拓扑荷符号的信息完全编码在相位里。实际操作中检查模拟程序是否正确的一个常用技巧是验证这一点修改l的正负号强度图必须保持不变相位图旋转方向必须反转。如果强度图变了说明代码里用了l而不是|l|或者是相位项exp(-ilφ)里遗漏了负号。p的影响更为微妙。p0时只有一个径向亮环中心暗斑被一个圆环包围。p1时径向方向多了一个零点光强图会呈现同心双环结构外环比内环稍暗。p越大径向零点越多光斑越“碎”。径向零点的位置和Laguerre多项式的根位置一一对应这可以用来验证递推函数是否正确——用roots函数求多项式的根再对比光强暗环的径向位置。3.3 从强度图反推光束结构的方法模拟做多了之后数值图案和物理图像会建立很强的对应关系。看一张HG强度图你能立刻说出n和m分别是几数x方向的暗线数量就是n数y方向的暗线数量就是m。看一张LG强度图数环的数量可以推断p数中心暗斑的大小可以大致估算|l|。这种“直觉”对实验数据分析很有帮助。我经常在模拟同一阶数的HG然后做傅里叶变换得到远场光斑。HG是直角坐标本征解其傅里叶变换仍然是HG只是阶数不变宽度反比变化。LG有类似性质近场涡旋光经过透镜聚焦到焦平面仍然保持涡旋结构拓扑荷不变。在Matlab里这个验证很简单对复振幅E做二维快速傅里叶变换fft2加上fftshift再取强度即可。如果你发现远场光斑中心出现了不该有的亮点通常说明近场涡旋相位没有正确初始化——数值模拟中常见的“伪基础模式”泄漏问题。4. 常见问题与实战排查记录4.1 光斑出现锯齿边界和混叠伪影网格采样密度不够时高阶HG模式的快速振荡项不能被有效采样光斑边缘会出现锯齿形伪影。解决方法是增加N或者先减小计算区域L把有限采样资源集中在光斑核心区域。我一般是先算一阶导数的最大空间频率也就是(H_n)的零点密度最高区域对应的局部条纹间距让采样间隔小于这个间距的三分之一。混叠伪影的另一个来源是区域截断。Hermite多项式在远离中心的区域增长很快但高斯包络会把它压下来。如果计算区域不够大截断边界处振幅不为零fft2之后在远场产生矩形衍射环。解决办法是加一个超高斯窗函数把计算区域边缘的场平滑压到零。window exp(-((X/w0).^8 (Y/w0).^8)); % 超高斯窗口 E_HG E_HG .* window;超高斯阶数8是我经验里比较合适的值它让边缘过渡足够平滑又不会明显改变光斑中心结构。4.2 LG中心强度不为零很多第一次模拟LG光束的人会发现理论上应该中心为0的甜甜圈图案中心却有一个极小的亮点或者暗点不明显。这几乎都是网格奇偶性导致的采样问题。如果N为偶数坐标网格里没有实际点落在(0,0)位置涡旋奇点的强度在四个邻近像素上的平均值不为零导致中心出现假亮度。解决办法很简单把N设为奇数保证x0和y0落在网格点上。同时配合fftshift和ifftshift保持坐标对齐。如果你使用了imagesc显示还需要确认像素边缘的对齐方式有时候微调一下显示范围就能让中心暗点变得干净。4.3 光斑中心相位跳变被误判为错误我收到过好几次这样的提问“为什么我的LG相位图有一条白色的线是不是程序错了”实际上那条线正是螺旋相位的自然结果。由于相位取值范围限制在(-π, π]从π到-π的突变出现在方位角φπ处这就是所谓的割线。模拟LG光束时这条割线必须存在且方向任意如果为了“美观”强行把这条线去掉反而破坏了涡旋相位的拓扑性质。判断相位模拟是否正确的标准不是“没有跳变”而是跳变线的位置和形态是否符合预期——它应该是一条从中心延伸出去的径向线而不是随机的碎片条纹。4.4 高阶LG模式递推发散广义拉盖尔多项式的递推在阶数很高p30且自变量很大时可能出现数值不稳定现象表现为光斑外围出现发散性亮点。这本质上是递推过程中的误差累积。解决办法有两个一是使用更高精度的韦德曼函数可以用符号计算辅助提取高精度初始值二是在物理参数上做限制p超过20的模式在真实激光器中本来也不常见。如果你的研究确实需要高阶模式建议改用正交多项式的高斯求积节点来构造或者利用Matlab的Hankel变换在频域计算。4.5 相位图显示时的颜色映射选择相位是周期量用线性颜色映射显示的时候-π和π附近虽然数值不同但颜色差异很大会让人误以为存在巨大的相位梯度。更合理的显示方式是用圆形颜色映射例如hsv或者自定义的cyclic colormap让首尾颜色自然衔接。Matlab里可以用imagesc(angle(E))配合colormap(hsv)并caxis([-pi pi])体现整个相位周期。hsv在衔接处比parula自然得多特别适合显示涡旋相位。Table: 常见问题速查现象可能原因解决方法光斑边缘锯齿采样率不足增加N或减小L中心出现亮点N为偶数奇点落在网格间隙改用奇数N强度图案矩形方框计算区域截断扩大L或加窗函数LG相位图有割线正常涡旋相位特征无需修改画图时说明强度不对称Hermite递推写错检查.*和递推初值远场中心异常亮近场相位中心不为零检查螺旋相位初始化5. 扩展应用模拟之外的实操心得5.1 从模拟到实验的衔接经验模拟和实验之间总是存在一道鸿沟但Matlab模拟好好利用能显著缩短实验调试时间。我举两个例子。第一个是空间光调制器加载相位图的生成。实验里制作涡旋光通常是在SLM上加载螺旋相位灰度图而这个灰度图可以直接用Matlab生成phi mod(l * Phi, 2*pi); gray_level round(phi / (2*pi) * 255); imwrite(uint8(gray_level), spiral_phase.png);这样生成的灰度图直接上传到SLM驱动软件就能用。l1时相位从0过渡到2πl2时整个圆周内相位走两圈显示为两个“扇叶”。我在实验前一定会先模拟这个相位图的样子确认灰度渐变连续无跳变避免上机才发现图案错误。第二个是透镜聚焦后中心暗斑尺寸的估算。我们知道涡旋光经过透镜后焦平面的甜甜圈半径可以通过数值模拟精确计算。结合透镜焦距f和入射光束尺寸直接用衍射积分公式算一次得到的暗斑直径可以用来选择光电探测器或CCD的拍摄参数避免在实验台上反复试探浪费时间。5.2 从HG到LG的转换模拟HG和LG之间的模式转化在理论上有一个幺正变换矩阵。简单说HG(1,0)和HG(0,1)两个模式的等权重组合可以合成为LG(0,1)模式。这是因为直角坐标的一阶偶模和奇模通过(\frac{1}{\sqrt{2}}(HG_{10}±iHG_{01}))可以得到拓扑荷为±1的LG模。这个操作在实验中可以用柱透镜模式转换器实现在模拟中则只需要复振幅的线性叠加E_HG10 hermite_poly(1, sq2*X/w0) .* hermite_poly(0, sq2*Y/w0) .* G; E_HG01 hermite_poly(0, sq2*X/w0) .* hermite_poly(1, sq2*Y/w0) .* G; E_LG_p0_l1 (E_HG10 1i * E_HG01) / sqrt(2);运行这段代码后你会发现合成的复振幅强度图呈现完美的甜甜圈形状而且中心严格为零。这一步做通过之后对模式变换的理解就不仅仅是公式层面的记忆而是有了数值上的直观验证。很多做光通信模式复用的同行就是靠这种模拟来设计模式复用器的相位分布。5.3 模拟结果的工程复用我习惯把所有模拟函数集中到一个脚本文件里参数放在文件头部用注释标明日期和修改内容。这个习惯帮了大忙因为几个月后回来看代码往往记不清某个参数当时为什么取这个值。比如我在注释里会写“2024-03-12为验证STED显微方案将N从512改为1024原来在l2时中心暗斑轮廓不够光滑”。这些记录对论文复现和组内代码交接都非常有用。对代码做版本管理也很必要。哪怕只是自己一个人用我也推荐用Matlab的Compare功能或者干脆用Git管理.m文件。尤其在修改递推公式或参数单位换算时如果没有版本控制一个隐藏的错误可能在几周后才会被察觉而那时你已经不知道该回退到哪一版了。5.4 坐标网格的隐藏坑单位换算单位换算是这类模拟里最容易出问题但最不容易被注意的环节。我见过有人把波长写在纳米束腰写在毫米位置坐标写在微米结果数值差了好几个数量级光斑画出来完全是一团糊。我的建议是全程统一使用国际单位制波长用米长度用米角度用弧度只在最后显示时转换成毫米或微米。这样公式里的(k2\pi/\lambda)、瑞利距离(z_R\pi w_0^2/\lambda)等物理量才不会出现比例因子错误。如果要做不同波长的对比模拟最好把lambda放在一个独立的变量里不要直接写进公式。这样替换参数时只需修改一行无需逐段排查。这也是代码可维护性的基本功。5.5 最后的调试心得我在调试光学模拟程序时有一个固定流程先算00阶HG确认它和普通高斯光斑完全一致再算高阶模式验证对称性最后算LG验证中心奇点和螺旋相位。这个“从低到高、从简单到复杂”的递进策略帮我快速定位错误是出在公式还是出在代码。如果你一上来就调试HG(5,7)出了问题根本不知道从哪查起。还有一个小技巧是打印模式的能量守恒校验。理想情况下对一个模式做fft2再反变换fft2回来前后的总能量应该完全一致。如果能量有损耗说明计算区域内光场泄漏或边界条件设置有问题。这个校验在编写更复杂的传播算法时特别有用强烈建议在代码里加一个断言检查。模拟HG和LG光束这件事本质上是一个把抽象的数学解变成可视化图像的“翻译”过程。一旦你理解了公式的每一项在物理上对应什么代码就只是实现细节了。我在实际项目中最大的体会是数值模拟最大的价值不在于“画得漂亮”而在于快速建立物理直觉让你在实验之前就已经“看见”了答案。希望这套思路和代码能帮你少走一些弯路。 Matlab模拟HG光束与LG涡旋光：从理论公式到工程应用
从研究生一年级第一次在实验室里看到空间光调制器把一束普通高斯光“拧”成一个甜甜圈形状开始我就对Hermite-GaussianHG和Laguerre-GaussianLG这两种光束模式产生了兴趣。后来做光束整形和光镊实验几乎所有需要模式结构的地方都离不开这两类解。说实话纯解析推导容易把人劝退但Matlab模拟可以很直观地把公式变成眼前跳动的光斑图案。这篇文章就把我自己常用的模拟思路、代码细节和踩坑记录整理出来写给正在做光学模拟、激光模式分析或相关课程设计的同学参考。先说清楚这是个什么东西Hermite-Gaussian光束是傍轴波动方程在直角坐标系下的本征解光斑呈矩形网格状排列Laguerre-Gaussian光束是同一方程在柱坐标系下的本征解带有螺旋相位、中心强度为零就是我们常说的涡旋光。两类模式在激光谐振腔设计、光镊、超分辨显微、光通信模式复用等领域都是核心角色。Matlab模拟的核心价值在于不需要搭真实光路就能研究参数变化对光束行为的影响也能直接验证实验方案里某个相位板或者空间光调制器加载的图案到底该长什么样。尤其适合光学工程、物理专业的学生和刚接触光场调控的工程师上手。下面我从理论公式讲起把代码实现、结果分析和问题排查都过一遍。整个思路是一边解释“为什么这么写”一边给出可以直接抄作业的Matlab代码。1. 理论公式与参数设计思路1.1 HG光束的数学表达与物理直觉Hermite-Gaussian光束在z平面上的复振幅可以写成[ E_{nm}(x,y,z)A\frac{w_0}{w(z)}H_n\left(\frac{\sqrt{2}x}{w(z)}\right)H_m\left(\frac{\sqrt{2}y}{w(z)}\right)\exp\left(-\frac{x^2y^2}{w^2(z)}\right)\exp\left(ik\frac{x^2y^2}{2R(z)}\right)\exp\left(-i(nm1)\psi(z)\right) ]初学者看到这个公式往往直接头皮发麻我拆开给你看。(H_n)和(H_m)是n阶和m阶厄米多项式决定了光场在x方向和y方向上的振荡结构。n0时没有横向振荡n1时有一个过零节点n2时有两个节点以此类推。所以HG光束的强度图案就像一张网格纸暗线对应于厄米多项式的零点。(w(z))是光束在传播距离z处的束半径(w_0)是束腰处的半径。它描述了光束横向尺寸的演化远离束腰时光斑会逐渐发散。(R(z))是波前曲率半径表示等相位面的弯曲程度。束腰处波前是平面R为无穷大远离束腰波前逐渐变成球面。(\psi(z)\arctan(z/z_R))是Gouy相位其中(z_R\pi w_0^2/\lambda)是瑞利距离。这个相位累加项在单纯观察横向光强分布时可以忽略但你如果要模拟干涉或者计算谐振腔模式就必须保留。直观理解HG光束可以这样想拿一块矩形的弹性膜不同振动模式就对应不同的HG阶数。00阶就是普通高斯光斑一个圆点10阶是左右两个瓣中间一条暗线01阶是上下两个瓣11阶是四瓣。阶数越高光斑里的“格子”越多能量分布的范围也越大。模拟这类光束时最关键的参数有三个波长(\lambda)、束腰半径(w_0)、模式阶数((n,m))。另外模拟的横向计算区域大小和采样点数直接决定了图案分辨率和混叠程度这个后面展开说。1.2 LG光束的数学表达与相位奇点Laguerre-Gaussian光束在柱坐标下表示为[ E_{pl}(r,\phi,z)A\frac{w_0}{w(z)}\left(\frac{\sqrt{2}r}{w(z)}\right)^{|l|}L_p^{|l|}\left(\frac{2r^2}{w^2(z)}\right)\exp\left(-\frac{r^2}{w^2(z)}\right)\exp\left(-il\phi\right)\exp\left(ik\frac{r^2}{2R(z)}\right)\exp\left(-i(2p|l|1)\psi(z)\right) ]这里出现了两个新的角色。(L_p^{|l|})是广义拉盖尔多项式p是径向指数决定了光强在径向上有几个亮环。p0时径向没有零点和暗环p1时会出现一个径向暗环光斑呈现“单环套双瓣”的结构。(l)是拓扑荷数也叫涡旋阶数它对应的相位因子(\exp(-il\phi))是LG光束的灵魂。旋转方位角(\phi)一圈相位改变(2\pi l)等相位面是一个螺旋面。中心处因为相位不确定振幅必须为零所以强度图案中心是暗的形成“甜甜圈”。如果你用过空间光调制器产生涡旋光你会发现在中心总有一个黑点那就是相位奇点。拓扑荷(l)的绝对值越大暗斑半径越大甜甜圈越“胖”。p的数值则影响径向次级结构比如p1时中心暗斑外围会再出现一个暗环。HG和LG之间不是孤立的它们都是同一个光学系统在不同坐标系下的完备基。用适当系数将一个方向的HG模式叠加可以合成LG模式。这个关系在复杂光场调控中很有用比如你想生成任意空间结构光束本质上就是在一个基底下做模式展开。Matlab模拟时理解这一点能帮你检查程序结果是否合理——比如00阶HG和00阶LG完全一样都是普通高斯光。1.3 为什么模拟而不是直接实验有人会问公式都有了直接拿激光器和空间光调制器做实验不就行了模拟的意义在哪里一个很现实的原因是成本。搭建一套完整的实验光路需要激光器、透镜组、空间光调制器、CCD相机设备动辄几十万而且每一次调节都要花大量时间。相比之下Matlab模拟只需要一台普通电脑5分钟就能把任意阶数的HG和LG光斑画出来。在项目预研阶段模拟可以快速验证想法这个光束经过透镜聚焦后光斑长什么样如果我在光路上加一个与LG相位共轭的波片能不能把涡旋光还原成高斯光这些问题在代码里修改几个参数就能得到初步答案再决定值不值得上实验。另一个原因是精度可控。实验里总存在像差、噪声、对准误差很难单独观察某一个参数变化的影响。模拟则可以精确控制(w_0)、(\lambda)、模式阶数单独剥离某个变量这是理论推导和实验之间很好的桥梁。我自己习惯把所有模拟结果保存成标准格式的图片和数据方便在写论文或者做汇报时直接引用。2. Matlab数值实现方案2.1 模拟参数的选取依据写Matlab代码之前先把模拟区域和采样参数确定下来。这里的参数选择直接决定结果是否可信我常用的配置如下波长(\lambda)取632.8 nm也就是常用氦氖激光器的波长。这个值在模拟中只影响瑞利距离和传播相位不改变横向光斑形状但后续计算衍射要用。束腰(w_0)取1 mm。这个尺寸和实验室常见的激光器光斑大小接近方便直观对照。横向计算区域L取20 mm × 20 mm边长是束腰的20倍。这个范围足够容纳到高阶模式的整个光斑也不会出现能量“溢出”边界的问题。采样点数N取512。这是一个比较好的平衡点既能看清精细结构计算速度也快。选择合适的计算区域很关键。如果L相对于束腰太小高阶模式的外围旁瓣会被截断强度图案出现明显的矩形边界伪影。如果L太大而N不够中心区域的采样稀疏光斑细节就没了。经验法则是L至少是最高阶模式整体宽度外包络的2到3倍采样间隔要小于最小结构的四分之一周期。坐标网格用meshgrid生成。注意这里我用的是物理坐标单位是米后续显示图像时再转换成毫米方便阅读。lambda 632.8e-9; % 波长单位m w0 1e-3; % 束腰半径单位m z 0; % 观察位置腰斑处 L 20e-3; % 计算区域边长单位m N 512; % 每个方向采样点数 x linspace(-L/2, L/2, N); y x; [X, Y] meshgrid(x, y); [R, Phi] cart2pol(X, Y); % 柱坐标下的径向距离和方位角这里使用cart2pol把直角坐标网格转换成柱坐标网格后面计算LG光束会大量用到径向距离(R)和方位角(\Phi)。很多人会直接写R sqrt(X.^2Y.^2)和Phi atan2(Y,X)效果一样但cart2pol可读性更好一点我习惯用它。束腰处z0时(w(z)w_0)(R(z)\to\infty)Gouy相位(\psi(z)0)公式会大幅简化。实际模拟中我最常在z0处画横向光斑因为这是最典型的观察面也最容易和实验中的CCD采集图像对照。如果要模拟传播过程再引入z循环即可。2.2 厄米多项式的高效递推实现Matlab的符号工具箱里确实有hermiteH函数可以求厄米多项式比如hermiteH(3, x)。但符号计算在512×512的网格上会非常慢而且还要求安装Symbolic Math Toolbox。我自己的项目里从来不用它做二维扫描而是用递推关系实现数值版本。厄米多项式满足递推公式[ H_0(x)1,\quad H_1(x)2x,\quad H_{n1}(x)2xH_n(x)-2nH_{n-1}(x) ]这个递推只需要几次数组乘法和加法比符号计算快几个数量级。我封装成一个函数function H hermite_poly(n, x) % Hermite多项式数值递推 % n: 多项式阶数x: 自变量矩阵 if n 0 H ones(size(x)); elseif n 1 H 2 * x; else H0 ones(size(x)); H1 2 * x; for k 1 : n-1 H2 2 * x .* H1 - 2 * k .* H0; H0 H1; H1 H2; end H H1; end end注意这段代码里2kH0中的k是循环变量代表当前递推步数。厄米多项式递推公式里的系数是2nn就是正在构造的多项式序号所以循环变量k从1跑到n-1时正好对应生成H2之前的那一步。实际用的时候这个函数返回的是和输入x同尺寸的矩阵因为x是meshgrid产生的二维数组ones(size(x))保证每一步都是矩阵运算。这里有一个容易犯的错误递推的初始值H0和H1必须是和x同样大小的矩阵不能用标量1和2*x代替否则矩阵维度对不上。有些初学同学在这里栽跟头报错说“矩阵维度必须一致”其实就是因为初始值没写成矩阵形式。2.3 广义拉盖尔多项式的递推实现广义拉盖尔多项式(L_p^\alpha(x))也有现成的递推公式这里(\alpha|l|)[ L_0^\alpha(x)1 ] [ L_1^\alpha(x)1\alpha-x ] [ L_{k1}^\alpha(x)\frac{(2k1\alpha-x)L_k^\alpha(x)-(k\alpha)L_{k-1}^\alpha(x)}{k1} ]实现的时候需要注意一个细节公式里的x是归一化变量在LG公式里取(x2r^2/w^2(z))。调用递推函数时传入的应当是整个归一化变量矩阵而不是径向距离本身。很多教程直接传r结果画出来的光斑径向结构完全不对就是这个原因。function L laguerre_poly(p, alpha, x) % 广义拉盖尔多项式递推 % p: 多项式阶数, alpha: 参数, x: 自变量矩阵 if p 0 L ones(size(x)); elseif p 1 L 1 alpha - x; else L0 ones(size(x)); L1 1 alpha - x; for k 1 : p-1 L2 ((2*k1alpha-x) .* L1 - (kalpha) .* L0) / (k1); L0 L1; L1 L2; end L L1; end end这段代码里除以(k1)是唯一一个除法操作用的右除/对矩阵来说相当于每个元素都除以同一个标量结果还是同尺寸矩阵。递推过程涉及x与L1的元素级相乘所以不要用而用.这是Matlab新手最常见的坑。如果你把.写成了Matlab会尝试做矩阵乘法矩阵尺寸不对就报错或者更奇怪的是在方阵情况下能“乘”出来但结果完全错误。我见过有人在代码里用*把整个结果算得面目全非还不自知的情况所以这里特意强调。2.4 HG光束复振幅计算与可视化有了递推函数HG光束的复振幅就是一行数学式的直接翻译。在束腰位置z0处(w(z)w0)(R(z))不考虑Gouy相位为0公式简化为[ E_{nm}(x,y)H_n(\sqrt{2}x/w_0)H_m(\sqrt{2}y/w_0)\exp\left(-\frac{x^2y^2}{w_0^2}\right) ]归一化常数可以先不写因为最后显示强度图时都要做归一化处理。数值计算的代码如下n 2; m 1; % HG模式阶数 sq2 sqrt(2); wx sq2 * X / w0; wy sq2 * Y / w0; Hx hermite_poly(n, wx); Hy hermite_poly(m, wy); G exp(-(X.^2 Y.^2) / w0^2); E_HG Hx .* Hy .* G; I_HG abs(E_HG).^2; I_HG I_HG / max(I_HG(:)); % 归一化这里有几个关键点值得展开。第一厄米多项式自变量是(\sqrt{2}x/w(z))不是单纯x。这个缩放因子来源于高斯函数指数项的自变量匹配。如果漏掉(\sqrt{2})得到的“HG光束”只是高斯包络和普通厄米多项式的乘积并不是真正的傍轴本征模图案虽然长得相似但节点位置和亮斑间隔都是错的。我在给组里新人审代码时至少见过三次这个错误每次都是光斑看起来“有点怪但说不上哪里怪”。第二强度图是(|E|^2)不是E本身。复数振幅里包含了相位信息直接画abs(E)得到的是振幅分布不是光强。虽然两者的亮暗结构类似但物理含义不同写注释时要说清楚。第三归一化用max(I_HG(:))而不是max(I_HG)。前者找出整个矩阵的最大值后者只会找出每一列的最大值返回一个行向量。如果你后续要做多幅图对比每幅都用各自全局最大值归一化才能保证亮度比例一致。显示强度图时用imagesc配上colormap。我个人比较喜欢用parula这是Matlab默认的配色色觉障碍者也能分辨比老式jet更科学。坐标轴单位习惯转成毫米显示图像标题里注明(n,m)阶数。figure; imagesc(x*1e3, y*1e3, I_HG); axis image; colormap(parula); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(HG(%d,%d) intensity, n, m));axis image保证横纵坐标刻度一致不会因为窗口拉伸让圆形光斑变成椭圆。这个细节看着小但在论文插图里非常影响观感。2.5 LG光束复振幅计算与相位图LG光束在z0处的复振幅计算同样直接。唯一要注意的是奇点处理。公式中有一个因子((\sqrt{2}r/w_0)^{|l|})当(r0)且(l\neq 0)时数值为0所以中心强度为0。但如果网格的中心正好没有采样点比如N是偶数中心位置落在四个像素的交界处数值上中心光强不一定严格为0图案上会出现一个极小值而不是完美的黑洞。解决办法是让N取奇数或者在生成网格后用ifftshift调整坐标原点。我一般在模拟中用N511或513确保原点落在网格点上。LG光束的代码p 0; l 1; % LG模式径向阶数p拓扑荷l alpha abs(l); radial_arg 2 * R.^2 / w0^2; Lp laguerre_poly(p, alpha, radial_arg); E_LG (sqrt(2)*R/w0).^alpha .* Lp .* exp(-R.^2/w0^2) .* exp(-1i*l*Phi); I_LG abs(E_LG).^2; I_LG I_LG / max(I_LG(:));相位图比较复杂。LG光束的相位分布是螺旋状的直接使用angle函数会得到从-π到π的跳变在跳变处会出现一条明显的红蓝分界线这就是相位的2π卷绕。严格来说这不是错误螺旋相位本来就有一个2π相位不连续面实验里干涉测量时也会看到这条亮线。但如果你想让相位图看起来更连续、更像一个完整的螺旋面就需要对相位做解卷绕处理。Matlab里unwrap函数只能沿一维数组解卷绕不能直接处理二维相位图。如果你有Image Processing Toolbox可以用wrapToPi和相位解缠算法如果没有工具箱可以自定义一个简单的二维解缠绕实现或者干脆接受这种锯齿状显示。我自己的做法是显示原始angle相位图但在图注里注明“相位卷绕为[-π,π]区间”这样是规范和诚实的方式。很多期刊的论文插图也这么干。另外要记得画一下LG光束的相位图因为(l±1)的相位图长得很像一片风扇叶旋转方向相反。这是判断拓扑荷符号的最直观方法。模拟时改变l的正负相位图的旋转方向会反过来强度图则几乎不变。这个特性在光学实验中经常用来标定涡旋光的拓扑荷符号。3. 多模式模拟与结果解读3.1 不同阶数HG光束的对比写一个循环就能一口气生成一组HG(0,0)到HG(3,2)的强度图。这里我用subplot把多幅图放到同一个figure里figure(Position, [100 100 1200 800]); idx 1; for n 0 : 3 for m 0 : 2 wx sq2 * X / w0; wy sq2 * Y / w0; Hx hermite_poly(n, wx); Hy hermite_poly(m, wy); Enm Hx .* Hy .* G; Inm abs(Enm).^2; Inm Inm / max(Inm(:)); subplot(4, 3, idx); imagesc(x*1e3, y*1e3, Inm); axis image off; colormap(parula); title(sprintf(HG(%d,%d), n, m)); idx idx 1; end end运行后你会注意到几个规律总阶数(Nnm)相同的HG模式比如HG(2,0)和HG(1,1)它们光斑的整体空间范围差不多因为Gouy相位和光束扩散速度都只依赖(N1)x方向和y方向的亮斑数量独立HG(2,0)在x方向有三个峰在y方向只有一个峰而HG(1,1)在两个方向上各有两个峰。还有一个容易被忽略的性质HG模式在传播过程中会保持自己的形状不变只是整体尺寸缩放和附加一个Gouy相位。这就是为什么它们被称为谐振腔的“本征模”。如果模拟中看到形状随z变化通常是公式里的归一化缩放和相位项没写对。3.2 不同拓扑荷LG光束的结构差异同样方法可以生成p0时l从-3到3的LG模式对比图。你会发现|LG|的强度图案完全相同——l1和l-1的甜甜圈一样大一样亮唯一的区别是相位沿相反方向旋转。这表明LG光束的强度分布只依赖|l|而拓扑荷符号的信息完全编码在相位里。实际操作中检查模拟程序是否正确的一个常用技巧是验证这一点修改l的正负号强度图必须保持不变相位图旋转方向必须反转。如果强度图变了说明代码里用了l而不是|l|或者是相位项exp(-ilφ)里遗漏了负号。p的影响更为微妙。p0时只有一个径向亮环中心暗斑被一个圆环包围。p1时径向方向多了一个零点光强图会呈现同心双环结构外环比内环稍暗。p越大径向零点越多光斑越“碎”。径向零点的位置和Laguerre多项式的根位置一一对应这可以用来验证递推函数是否正确——用roots函数求多项式的根再对比光强暗环的径向位置。3.3 从强度图反推光束结构的方法模拟做多了之后数值图案和物理图像会建立很强的对应关系。看一张HG强度图你能立刻说出n和m分别是几数x方向的暗线数量就是n数y方向的暗线数量就是m。看一张LG强度图数环的数量可以推断p数中心暗斑的大小可以大致估算|l|。这种“直觉”对实验数据分析很有帮助。我经常在模拟同一阶数的HG然后做傅里叶变换得到远场光斑。HG是直角坐标本征解其傅里叶变换仍然是HG只是阶数不变宽度反比变化。LG有类似性质近场涡旋光经过透镜聚焦到焦平面仍然保持涡旋结构拓扑荷不变。在Matlab里这个验证很简单对复振幅E做二维快速傅里叶变换fft2加上fftshift再取强度即可。如果你发现远场光斑中心出现了不该有的亮点通常说明近场涡旋相位没有正确初始化——数值模拟中常见的“伪基础模式”泄漏问题。4. 常见问题与实战排查记录4.1 光斑出现锯齿边界和混叠伪影网格采样密度不够时高阶HG模式的快速振荡项不能被有效采样光斑边缘会出现锯齿形伪影。解决方法是增加N或者先减小计算区域L把有限采样资源集中在光斑核心区域。我一般是先算一阶导数的最大空间频率也就是(H_n)的零点密度最高区域对应的局部条纹间距让采样间隔小于这个间距的三分之一。混叠伪影的另一个来源是区域截断。Hermite多项式在远离中心的区域增长很快但高斯包络会把它压下来。如果计算区域不够大截断边界处振幅不为零fft2之后在远场产生矩形衍射环。解决办法是加一个超高斯窗函数把计算区域边缘的场平滑压到零。window exp(-((X/w0).^8 (Y/w0).^8)); % 超高斯窗口 E_HG E_HG .* window;超高斯阶数8是我经验里比较合适的值它让边缘过渡足够平滑又不会明显改变光斑中心结构。4.2 LG中心强度不为零很多第一次模拟LG光束的人会发现理论上应该中心为0的甜甜圈图案中心却有一个极小的亮点或者暗点不明显。这几乎都是网格奇偶性导致的采样问题。如果N为偶数坐标网格里没有实际点落在(0,0)位置涡旋奇点的强度在四个邻近像素上的平均值不为零导致中心出现假亮度。解决办法很简单把N设为奇数保证x0和y0落在网格点上。同时配合fftshift和ifftshift保持坐标对齐。如果你使用了imagesc显示还需要确认像素边缘的对齐方式有时候微调一下显示范围就能让中心暗点变得干净。4.3 光斑中心相位跳变被误判为错误我收到过好几次这样的提问“为什么我的LG相位图有一条白色的线是不是程序错了”实际上那条线正是螺旋相位的自然结果。由于相位取值范围限制在(-π, π]从π到-π的突变出现在方位角φπ处这就是所谓的割线。模拟LG光束时这条割线必须存在且方向任意如果为了“美观”强行把这条线去掉反而破坏了涡旋相位的拓扑性质。判断相位模拟是否正确的标准不是“没有跳变”而是跳变线的位置和形态是否符合预期——它应该是一条从中心延伸出去的径向线而不是随机的碎片条纹。4.4 高阶LG模式递推发散广义拉盖尔多项式的递推在阶数很高p30且自变量很大时可能出现数值不稳定现象表现为光斑外围出现发散性亮点。这本质上是递推过程中的误差累积。解决办法有两个一是使用更高精度的韦德曼函数可以用符号计算辅助提取高精度初始值二是在物理参数上做限制p超过20的模式在真实激光器中本来也不常见。如果你的研究确实需要高阶模式建议改用正交多项式的高斯求积节点来构造或者利用Matlab的Hankel变换在频域计算。4.5 相位图显示时的颜色映射选择相位是周期量用线性颜色映射显示的时候-π和π附近虽然数值不同但颜色差异很大会让人误以为存在巨大的相位梯度。更合理的显示方式是用圆形颜色映射例如hsv或者自定义的cyclic colormap让首尾颜色自然衔接。Matlab里可以用imagesc(angle(E))配合colormap(hsv)并caxis([-pi pi])体现整个相位周期。hsv在衔接处比parula自然得多特别适合显示涡旋相位。Table: 常见问题速查现象可能原因解决方法光斑边缘锯齿采样率不足增加N或减小L中心出现亮点N为偶数奇点落在网格间隙改用奇数N强度图案矩形方框计算区域截断扩大L或加窗函数LG相位图有割线正常涡旋相位特征无需修改画图时说明强度不对称Hermite递推写错检查.*和递推初值远场中心异常亮近场相位中心不为零检查螺旋相位初始化5. 扩展应用模拟之外的实操心得5.1 从模拟到实验的衔接经验模拟和实验之间总是存在一道鸿沟但Matlab模拟好好利用能显著缩短实验调试时间。我举两个例子。第一个是空间光调制器加载相位图的生成。实验里制作涡旋光通常是在SLM上加载螺旋相位灰度图而这个灰度图可以直接用Matlab生成phi mod(l * Phi, 2*pi); gray_level round(phi / (2*pi) * 255); imwrite(uint8(gray_level), spiral_phase.png);这样生成的灰度图直接上传到SLM驱动软件就能用。l1时相位从0过渡到2πl2时整个圆周内相位走两圈显示为两个“扇叶”。我在实验前一定会先模拟这个相位图的样子确认灰度渐变连续无跳变避免上机才发现图案错误。第二个是透镜聚焦后中心暗斑尺寸的估算。我们知道涡旋光经过透镜后焦平面的甜甜圈半径可以通过数值模拟精确计算。结合透镜焦距f和入射光束尺寸直接用衍射积分公式算一次得到的暗斑直径可以用来选择光电探测器或CCD的拍摄参数避免在实验台上反复试探浪费时间。5.2 从HG到LG的转换模拟HG和LG之间的模式转化在理论上有一个幺正变换矩阵。简单说HG(1,0)和HG(0,1)两个模式的等权重组合可以合成为LG(0,1)模式。这是因为直角坐标的一阶偶模和奇模通过(\frac{1}{\sqrt{2}}(HG_{10}±iHG_{01}))可以得到拓扑荷为±1的LG模。这个操作在实验中可以用柱透镜模式转换器实现在模拟中则只需要复振幅的线性叠加E_HG10 hermite_poly(1, sq2*X/w0) .* hermite_poly(0, sq2*Y/w0) .* G; E_HG01 hermite_poly(0, sq2*X/w0) .* hermite_poly(1, sq2*Y/w0) .* G; E_LG_p0_l1 (E_HG10 1i * E_HG01) / sqrt(2);运行这段代码后你会发现合成的复振幅强度图呈现完美的甜甜圈形状而且中心严格为零。这一步做通过之后对模式变换的理解就不仅仅是公式层面的记忆而是有了数值上的直观验证。很多做光通信模式复用的同行就是靠这种模拟来设计模式复用器的相位分布。5.3 模拟结果的工程复用我习惯把所有模拟函数集中到一个脚本文件里参数放在文件头部用注释标明日期和修改内容。这个习惯帮了大忙因为几个月后回来看代码往往记不清某个参数当时为什么取这个值。比如我在注释里会写“2024-03-12为验证STED显微方案将N从512改为1024原来在l2时中心暗斑轮廓不够光滑”。这些记录对论文复现和组内代码交接都非常有用。对代码做版本管理也很必要。哪怕只是自己一个人用我也推荐用Matlab的Compare功能或者干脆用Git管理.m文件。尤其在修改递推公式或参数单位换算时如果没有版本控制一个隐藏的错误可能在几周后才会被察觉而那时你已经不知道该回退到哪一版了。5.4 坐标网格的隐藏坑单位换算单位换算是这类模拟里最容易出问题但最不容易被注意的环节。我见过有人把波长写在纳米束腰写在毫米位置坐标写在微米结果数值差了好几个数量级光斑画出来完全是一团糊。我的建议是全程统一使用国际单位制波长用米长度用米角度用弧度只在最后显示时转换成毫米或微米。这样公式里的(k2\pi/\lambda)、瑞利距离(z_R\pi w_0^2/\lambda)等物理量才不会出现比例因子错误。如果要做不同波长的对比模拟最好把lambda放在一个独立的变量里不要直接写进公式。这样替换参数时只需修改一行无需逐段排查。这也是代码可维护性的基本功。5.5 最后的调试心得我在调试光学模拟程序时有一个固定流程先算00阶HG确认它和普通高斯光斑完全一致再算高阶模式验证对称性最后算LG验证中心奇点和螺旋相位。这个“从低到高、从简单到复杂”的递进策略帮我快速定位错误是出在公式还是出在代码。如果你一上来就调试HG(5,7)出了问题根本不知道从哪查起。还有一个小技巧是打印模式的能量守恒校验。理想情况下对一个模式做fft2再反变换fft2回来前后的总能量应该完全一致。如果能量有损耗说明计算区域内光场泄漏或边界条件设置有问题。这个校验在编写更复杂的传播算法时特别有用强烈建议在代码里加一个断言检查。模拟HG和LG光束这件事本质上是一个把抽象的数学解变成可视化图像的“翻译”过程。一旦你理解了公式的每一项在物理上对应什么代码就只是实现细节了。我在实际项目中最大的体会是数值模拟最大的价值不在于“画得漂亮”而在于快速建立物理直觉让你在实验之前就已经“看见”了答案。希望这套思路和代码能帮你少走一些弯路。