简介拉盖尔-高斯光束的MATLAB仿真资源围绕p01模态展开面向激光物理、光学工程及量子光学领域的学生与科研人员用于解决高阶模式光束难以直观构建与分析的问题。压缩包共2个文件一个MATLAB脚本负责生成拉盖尔-高斯光束的复数光场并完成强度计算另一张图像则展示对应强度分布清晰呈现径向亮环与中心暗斑甜甜圈形状资源整体仅66KB轻量便捷。已有2190人学习下载学习热度可观借助MATLAB的矩阵运算与可视化能力读者可自由调整模式参数观察分布演化从而深入理解轨道角动量、径向阶数与相位奇点的物理意义为光学仿真入门提供直观范例。无论是课程演示、毕业设计还是课题预研都能快速获得可复现的仿真结果为光通信、量子光学或粒子操控等研究提供扎实的仿真基础。1. LG01 模式在 MATLAB 里到底算什么——先解释一下你下载到的东西从 LG01.rar 里解压出来通常只有两个文件LG01.m 和 01.png。很多人在 MATLAB 里跑完脚本看到一张中心带黑孔的“甜甜圈”强度图就收工了但真正值得搞明白的是这个黑孔不是数值误差而是光束携带轨道角动量的直接证据。按常见的 LG_{p,l} 命名习惯LG01 的意思是径向指数 p0、螺旋指数 l1也就是单位拓扑荷涡旋光。相比普通高斯光束它多了螺旋相位项 exp(-iφ)这直接决定中心场的相位不确定于是强度必须为零。这篇文章会把 LG01 从公式推导到 MATLAB 实现逐层拆开顺带讲清楚版本兼容、网格采样和验证方法既能当 matlab 可视化大学物理的参考资料也能直接当 matlab 图像处理大作业的仿真模板。2. LG 光束的数学结构拉盖尔多项式与轨道角动量如何耦合LG01 不是人为拼凑出来的波函数而是波动方程在圆柱坐标下的本征解。理解这一点你在 MATLAB 里写的每一行代码都会对应明确的物理量反之只抄公式不关心归一化画出来的图永远和 01.png 有偏差。2.1 从高斯光束到 LG 模式p 和 l 是两个独立自由度高斯光束在直角坐标里可以写成厄米多项式与高斯函数乘积的叠加这类解叫厄米-高斯模。在圆对称的光纤、谐振腔或柱对称光学系统中更自然的基函数是拉盖尔-高斯模。它把横向分布拆成径向部分和角向部分径向部分由拉盖尔多项式决定角向部分是一个纯相位因子 exp(-ilφ)。p 代表径向节线数也叫径向模式阶数l 代表角向绕轴一整圈积累的相位是 2πl 的整数倍。这个 l 在量子光学里就是每个光子携带的轨道角动量 lħ在经典光学里对应螺旋波前。对 LG01 来说p0 表示径向方向没有节线沿半径从中心向外先增后减只有一个环。l1 表示相位沿方位角旋转一圈变化 2π也就是拓扑荷为 1。这种模态和 p1、l0 的 LG10 看起来都是空心环但成因完全不同LG01 的空心来自相位奇点LG10 的空心来自拉盖尔多项式在 r0 处的因子过零。判断一个空心束到底是哪一类只需要观察它的相位图——有涡旋的一定存在绕轴 2π 的相位跳变。2.2 LG01 在 z0 处的复振幅公式取束腰位置 z0忽略曲率项和 Gouy 相位LG 模式的一般形式可以写成$$U_{p,l}(r,\phi) C_{p,l} \cdot \left(\frac{\sqrt{2}r}{w_0}\right)^{|l|} \cdot L_p^{|l|}\left(\frac{2r^2}{w_0^2}\right) \cdot \exp\left(-\frac{r^2}{w_0^2}\right) \cdot e^{-il\phi}$$其中 L_p^{|l|} 是广义拉盖尔多项式w0 是束腰半径。对 p0、l1拉盖尔多项式 L_0^1(x)1公式直接简化成$$U_{01}(r,\phi) \frac{2}{\sqrt{\pi}, w_0^2}, r, e^{-r^2/w_0^2}, e^{-i\phi}$$注意 r 前面只有一次方这是因为 |l|1。这个 r 因子是涡旋光的标志它保证 r0 处振幅严格为 0不会因为坐标网格对称性而出现残留亮点。强度 |U01|² 的最大值位置满足 d|U|²/dr0解出来是 r w0/√2约为束腰的 0.707 倍。这个值非常有用后面在 MATLAB 里验证 01.png 的环形半径时可以直接对照。实际写代码时我建议先用上面这个 p0 的解析形式做版本等图对了再改成通用形式。因为 p0 时没有拉盖尔多项式的数值抖动如果环位置、中心暗斑、环宽都正确说明你的网格和坐标写法没问题这时再把一般公式加进去出问题就知道是多项式参数的问题。2.3 相位奇点为什么中心越“空”越说明螺旋是正确的涡旋光束中心强度为零常被初学者误以为是强度被屏蔽了实际原因是相位在中心处没有定义。想象沿 r0 绕一圈相位从 0 连续变化到 2π任何一点都要给一个确定值但对同一个物理点来说相位必须是单值的。唯一能自洽的结果是轴上点强度为 0场在这里“消失”问题自然回避。这就是相位奇点也叫涡旋核。在 MATLAB 里这个奇点会在 angle() 函数画相位图时表现为 2π 到 -2π 的跳变线从中心向外延伸一条线。这不是 bug是拓扑荷为 1 的直接证明。如果这条跳变线出现在光束中心以外的其他位置说明生成光束的螺旋相位中心没有对准网格原点即坐标 Phi atan2(Y, X) 的奇点和波前涡旋中心错位仿真结果会整体偏移。2.4 LG_{p,l} 命名约定的坑不同教材对模式名的写法不完全一致。有的写成 LG_{l,p}有的写成 LG_{p,l}有的直接用 LGpl 两个数字连写比如 LG01、LG10。你手里的资源如果标的是 LG01在 MATLAB 代码里几乎可以确定是 p0、l1。但如果在论文里看到 LG_{2,1}一定要先看这个作者前面有没有定义过下标顺序我见过不止一次在模式复用仿真里把 p 和 l 写反导致整个正交性分析全错。3. MATLAB 实现 LG01网格、复数场与可视化这一章给出可以直接跑的 MATLAB 代码。核心思路是用 meshgrid 建立直角坐标网格再用 cart2pol 转成极坐标生成复数光场后离散采样画图。我在写这段代码时特意不依赖任何光学工具箱只用基础函数保证任何版本的 MATLAB 都能运行。3.1 建立径向和角向网格clear; clc; lambda 632.8e-9; % 氦氖激光波长 632.8 nm w0 0.5e-3; % 束腰半径 0.5 mm Lx 4e-3; % 仿真窗口边长 4 mm N 1024; % 网格数建议 2 的幂次 x linspace(-Lx/2, Lx/2, N); [X, Y] meshgrid(x, x); [Phi, R] cart2pol(X, Y);这段代码的要点是窗口边长要取束腰的 6 到 10 倍太小会把环形截断太大则中心区域采样不足。R 的范围从 0 到约 2.8 mmLG01 的环峰值在 0.35 mm 左右折算到网格上大约有 90 个像素跨越环半径足够画出平滑的甜甜圈。N1024 时 X、Y、R、Phi 各占约 8 MB 内存四个矩阵加起来 32 MB普通电脑无压力。如果还想加快速度N512 时环半径也有约 45 个像素视觉上差别不大。3.2 LG01 直接生成代码不依赖任何工具箱E01 2 / (sqrt(pi) * w0^2) ... .* R .* exp(-R.^2 / w0^2) ... .* exp(-1i * Phi); I01 abs(E01).^2; I01 I01 / max(I01(:)); % 归一化方便显示 figure(Name,LG01 intensity); imagesc(x*1e3, x*1e3, I01); axis image; colormap gray; xlabel(x / mm); ylabel(y / mm); title(LG01 intensity (p0, l1));注意最后 exp(-1i*Phi) 前面的负号。角向相位因子在不同文献里有 e^{-ilφ} 和 e^{ilφ} 两种约定MATLAB 里用 exp(-1i * ell * Phi) 对应最常见的光学教材。如果改用正号产生的涡旋旋转方向会反过来强度图看不出区别但相位图会呈现反向螺旋。对单纯的仿真显示正负号无伤大雅如果后面要和实验全息图叠加干涉必须提前统一符号约定否则做出来的全息图会生成拓扑荷相反的涡旋。这里的 R 矩阵已经包含了 r^1 因子所以 R0 处 R.*exp(...) 严格为 0不会产生 NaN。直接对强度图看环峰位置在 r w0/sqrt(2) ≈ 0.354 mm。可以用下面的代码沿径向截线验证mid ceil(N/2); profile I01(mid, mid:end); r_axis (0:N-mid) * (Lx/N) * 1e3; % 从中心到边缘的距离单位 mm [Peak, idx] max(profile); fprintf(峰位半径: %.4f mm\n, r_axis(idx));如果输出的峰位半径接近 0.35 mm说明网格和公式都正确。这一步是整个仿真里最重要的自检也是我判断一份 LG01 资源代码是否靠谱的第一步。3.3 推广到任意 p 和 l 的通用版本LG01 只是特例实际做模式复用或轨道角动量通信仿真时通常需要生成任意 LG_{p,l}。这时要用广义拉盖尔多项式。MATLAB 从 R2017a 开始提供 laguerreL 函数属于符号数学工具箱。如果你只有基础版 MATLAB直接调用会报错 Undefined function laguerreL。function Epl lg_mode(p, ell, R, Phi, w0) u sqrt(2) * R / w0; Lp laguerreL(p, abs(ell), u.^2); Epl sqrt(2 * factorial(p) / (pi * factorial(p abs(ell)))) / w0 ... .* u.^abs(ell) .* Lp .* exp(-u.^2 / 2) .* exp(-1i * ell * Phi); endlaguerreL 接受矩阵输入时会逐个元素计算但对 1024×1024 的矩阵符号计算速度很慢而且结果类型是 sym直接绘图的 imagesc 不一定能处理。建议在调用后加一行 double()Epl double(Epl);。没有工具箱的话可以先用第 3.2 节的 p0 简化公式p0 时 L_0^abs(ell)1并不需要多项式真正需要拉盖尔多项式的一般是高阶 p比如 p1 时 L_1^m(x)1m-x也可以手写。在同一个窗口里叠加比较 LG00、LG10、LG01 三个模式的截面是理解这三个模式差别的快捷方式。LG00 中心亮、单峰LG10 中心暗但相位没有螺旋LG01 中心暗且有螺旋相位。三张图放在一起旋涡光的“拓扑荷”概念就直观了。3.4 强度图、相位图和三维曲面图的绘制细节强度图可以直接 imshow但相位图必须用 angle 函数而不是 abs且不能直接取实部或虚部画否则看不出 2π 跳变规律。figure(Name,LG01 phase); imagesc(angle(E01)); axis image; colormap hsv; colorbar; title(phase of LG01); figure(Name,LG01 3D); surf(x*1e3, x*1e3, I01, EdgeColor, none); view(45, 60); xlabel(x/mm); ylabel(y/mm); zlabel(intensity);相位图应该看到以中心为原点的螺旋条纹且沿绕中心的一周颜色从红色到蓝色再到红色渐变一圈刚好是 2π。用 hsv colormap 的原因正是它的色环首尾相接对应角度周期性。如果你用 jet 或 parula看到的跳变线会不明显。三维曲面图中如果直接用 surf 画原始强度会在中心生成一个很高的尖刺形状因为surf默认以矩阵行列作为 x/y 坐标物理坐标不同轴会拉伸变形。上面代码里显式传入 x 和 y 物理坐标view(45,60) 给了一个合适的倾斜视角。4. 运行调试与参数调整从 R2023b 到新版 MATLAB 的兼容细节下载的 LG01.rar 里的脚本在不同版本 MATLAB 下表现会有差异。下面把常见问题、参数表、边界情况列清楚按顺序排查基本能解决。4.1 解压后最常见的四类报错第一类是 laguerreL 未定义这属于缺少符号数学工具箱。解决方法已经在上一章给出p0 直接展开p1 手写一次多项式。第二类是内存不足。N 被改成 2048 以上时X、Y、R、Phi 和 E 五个 double 矩阵就要 160 MB如果电脑内存不大MATLAB 会提示 Out of memory。LG01 仿真 512 到 1024 网格完全够用没必要追求高分辨率。第三类是图像显示全黑或全白。这通常是没对强度做归一化。LG01 中心强度为 0环峰值很小如果直接 imshow(I01) 会自动把最小值映射成黑色最大值映射成白色但乘了 1e-4 之类的系数动态范围错位。写 I01 I01 / max(I01(:)) 即可。第四类是相位图出现随机雪花点。这往往出现在 R 很小的区域R 虽不为 0 但接近网格分辨率angle(E) 的值受浮点误差主导。LG01 的强度在高斯因子衰减下本来就趋近于 0相位无意义。掩盖方式是只显示强度大于阈值 0.01 的位置即对相位图乘以一个掩膜mask I01 0.01; phase_plot angle(E01) .* mask; phase_plot(phase_plot 0) NaN; imagesc(phase_plot);把强度为零区域设成 NaNcolormap 会把 NaN 渲成空白视觉上比一大片雪花点干净。4.2 网格采样与窗口宽度怎么搭配合适网格参数是这只仿真最容易被忽略的部分。窗口太长LG01 环形只占一个小点窗口太短高斯尾部被硬截断傅里叶传播后会出现方形衍射纹。我常用的组合是窗口边长取 8w0N 取 1024。这样空间采样步长 Δx 8w0/1024 ≈ 0.0078w0最大可记录空间频率对应离散网格的高频截止在传播仿真里能覆盖到足够的角谱范围。如果只画强度图窗口 6w0、N 512 也够如果要插入拉盖尔多项式的高阶项比如 p5建议 N1024避免环与环之间的条纹因采样不足而混叠。4.3 参数调整参考表参数推荐范围对 LG01 的影响参考依据w00.2–2 mm决定环半径 rw0/√2 和高斯衰减速度可从 01.png 中测量环径反推lambda400–1550 nm只影响传播相位不影响 z0 强度632.8 nm 最常见N512–2048过低导致环不圆过高内存翻倍1024 是均衡点Lx6w0–10w0截断误差与视觉占比的平衡8w0 最稳归一化方式峰值归一化或能量归一化对比相对强度时用峰值积分能量时用全套因子与 01.png 比较时用峰值检测环半径是判断代码与 01.png 是否匹配的关键步骤。写一个自动估半径的函数先找强度最大值所在像素位置然后计算该点到中心距离再乘以 Δx。如果 p 和 l 为 0、1实测半径应在 w0/√2 的 5% 误差范围内。超过 10% 说明网格中心偏了或者用了 meshgrid 而不是 ndgrid 导致 x/y 方向尺度不一致。4.4 用 01.png 做结果对照时的注意事项01.png 是一张静态栅格图压缩成 png 时有损但强度分布基本保留。和它对比时最稳妥的做法不是像素级比大小而是比两个特征中心是否出现零值单像素暗区以及沿径向的强度曲线是否在约 0.707w0 处取最大值。因为不同脚本的归一化方式不同峰值强度被缩放到 0~255 的范围直接对比像素值没有意义。如果发现你的仿真图和 01.png 的暗环粗细不一致优先确认 w0 取了多少。常见问题是 w01 但没写单位导致和取 0.5e-3 的结果差 2000 倍。这里建议直接在文件头部注释里写清单位避免后续换机器重新运行时出现相同的坑。5. 从 LG01 延伸到轨道角动量通信复用、解调与传播验证LG01 的价值在于它的 l1 可以作为一种正交信道在自由空间光通信领域实现模分复用。这一章看两个更实用的扩展多个 LG 模式叠加和通过角谱法观察传播后拓扑荷是否保持。5.1 叠加多个 LG 模式做 OAM 复用样例E00 lg_mode(0, 0, R, Phi, w0); % 高斯模式 E01 lg_mode(0, 1, R, Phi, w0); % 涡旋模式 E10 lg_mode(1, 0, R, Phi, w0); % 径向高阶模式 E_mix E00 E01 E10; % 三路复用 I_mix abs(E_mix).^2; I_mix I_mix / max(I_mix(:)); imagesc(I_mix); axis image; colormap gray;上面代码用三种模式叠加得到的光场已不再是一个干净的环中心不会完全为 0。这是因为 E00 和 E10 都有非零的轴上分量混合后轴上强度取决于它们之间干涉相位。这个现象在 OAM 复用实验里很常见接收端不是直接从强度图上读出模式而是用匹配滤波做相关解调coeff01 sum(sum(conj(lg_mode(0,1,R,Phi,w0)) .* E_mix)) * (Lx/N)^2;coeff01 是接收场与 LG01 模板的内积。由于 LG 模式在连续域正交离散网格上只要采样够细不同模式之间的 inner product 会接近 0而同一模式的内积接近 1。这个系数就是 OAM 信道的解调输出。你可以把该项作为判别 LG01 存在与否的依据——比看强度图精确得多。5.2 用角谱法传播一段距离验证涡旋核保留LG01 在自由空间传播时涡旋结构不会被抹掉但半径会随距离增大。用角谱法可以观察这一过程z 2e-3; % 传播距离 2 m fx (-N/2:N/2-1) / Lx; [FX, FY] meshgrid(fx, fx); H exp(1i * 2*pi/lambda * z .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); Ef ifft2(fft2(E01) .* H); I_far abs(Ef).^2;这段代码的核心是传递函数 H它对方位角的傅里叶谱乘上不同的相位等效于自由空间衍射。注意 sqrt 内可能出现负值说明该高频分量是倏逝波在远场忽略即可。用上面这组参数时2 m 后 LG01 环半径会从 0.35 mm 涨到约 1 mm 以上但中心仍然是零强度点。用上一章掩膜方法画相位图仍能看到从中心发出的 2π 跳变线证明传播不改变拓扑荷。5.3 用某个核心指标检验仿真与实验一致性光学实验里测量涡旋光最常用的是环形场干涉法也就是让 LG01 和一束平面波倾斜干涉干涉条纹在涡旋中心会出现一个分叉点。在 MATLAB 里加一束倾斜平面波 do exactly thisphase_tilt exp(1i * 2*pi/lambda * 0.1 * X); % 沿 x 轴加 0.1 rad/um 波矢 I_int abs(E01 phase_tilt).^2; I_int I_int / max(I_int(:)); imagesc(I_int); axis image; colormap gray;得到的条纹图中心会有一个明显的叉延叉的位置就是涡旋核心。把这个仿真图和实验干涉照片放一起对比叉的方向和条纹间距就能判断实际生成的涡旋是 l 还是 -l。这样一套方法比单纯看 01.png 的“中间有个洞”要可验证得多也是把这条 LG01 MATLAB 例程从演示代码升格成实验预判工具的关键一步。本文还有配套的精品资源点击获取