简介一个基于MATLAB模拟涡旋光束的源码包面向光学、物理相关专业的研究者与高年级学生。涡旋光束凭借轨道角动量特性在量子通信、光镊与粒子操控中有广泛应用该程序通过构造螺旋相位函数并计算复振幅可直观呈现不同阶数涡旋光的光强分布与相位结构适合作为理解OAM物理本质的入门工具。压缩包体积仅30KB共含3个文件1个m主程序负责生成螺旋相位、执行二维傅里叶变换并绘制图像2个jpg为典型参数下的模拟结果图可帮助校对运行效果。目前已有894人学习下载配套代码逻辑清晰涵盖相位构造、振幅计算、fft2频域处理和图像显示等关键环节读者可自行修改涡旋阶数、波长等参数快速观察光场变化为进一步开展量子通信或微粒操控相关科研打下基础。1. 涡旋光束在 MATLAB 里怎么入手先想清楚拓扑荷和相位屏涡旋光束的光场携带一个螺旋相位因子 exp(ilθ)l 叫拓扑荷表示绕光轴一圈相位改变 2πl。这个看起来只是改了一下波前的操作实际把轨道角动量 OAM 也带上了OAM 复用通信、光镊和超分辨显微都在这上面做文章。在 MATLAB 里做涡旋光束仿真目标不是把一段现成代码跑出几张彩虹图而是能控制相位屏、拓扑荷、传播距离和采样条件让每一张图的特征都对得上物理图像。如果你手上拿着从某个 rar 包里解出来的 vortex beam 程序先别急着整套改参数把下面这套建模、实现和验证路径过一遍大部分包里的坑都能在这些步骤里定位到。2. 涡旋光束复场建模从螺旋相位因子到 LG 模的 MATLAB 实现2.1 为什么是 exp(1iltheta)纯相位模型与 LG 模的取舍涡旋光束 MATLAB 仿真的最小描述是复振幅 E(r,θ) A(r) exp(i l θ)。A(r) 是振幅包络纯相位涡旋把包络设成常数或高斯包络适合做空间光调制器 SLM 的相位掩膜仿真LG拉盖尔-高斯模则把 A(r) 写成带径向节数的解析形式更适合传播演化和 OAM 谱分析。LG_p^l 在 p0 时可简化为 A(r) C (√2 r/w0)^|l| exp(-r²/w0²)。注意根号 2 这个系数很多网上程序会漏掉虽然看不出来但它直接影响暗核尺寸和干涉条纹的叉形位置。拿到 vortex beam.rar 后第一步不是运行而是先定位它用的是纯相位模型还是 LG 模。两者的区别见下表。模型强度分布适用场景典型写法纯相位涡旋均匀或高斯包络中心不一定是零SLM 掩膜、干涉检测exp(1i*l*Theta)LG 模p0中心暗环半径与 w0、l相关判断方法很简单把强度图截一条过中心的横线看最小值位置。纯相位模型中心可能不明显变暗而 LG 模一定有一条接近零的暗带。2.2 第一个能跑起来的涡旋光束 MATLAB 程序LG 模生成与相位图下面这个脚本直接生成拓扑荷 l2 的 LG 模并同时输出强度图和相位图。lambda 632.8e-9; % 波长 632.8nm常用 He-Ne 激光 w0 1e-3; % 束腰半径 1mm l 2; % 拓扑荷整数可正可负 N 512; % 网格像素数建议偶数 L 5e-3; % 计算区域边长 5mm x linspace(-L/2, L/2, N); [X, Y] meshgrid(x, x); [Theta, R] cart2pol(X, Y); % 注意返回顺序第一个是角度 rho R / w0; LG_amp (sqrt(2)*rho).^abs(l) .* exp(-rho.^2); LG LG_amp .* exp(1i * l * Theta); figure subplot(1, 2, 1); imagesc(x*1e3, x*1e3, abs(LG).^2); axis image; colormap hot; title(强度图); subplot(1, 2, 2); imagesc(x*1e3, x*1e3, angle(LG)); axis image; caxis([-pi pi]); colormap hsv; title(相位图);cart2pol返回结果是[Theta, R]而不是[R, Theta]这个顺序写反会出现奇怪的镜像相位。angle显示的相位被包裹到 [-π, π]所以拓扑荷越大相位图上的黑白跳变线越多一条表示绕一圈相位回绕一次。强度图虽然视觉上和实验接近但不能作为拓扑荷判据必须看相位图或干涉图。2.3 像素与奈奎斯特条件为什么 l 一大相位图全糊掉涡旋相位在半径 r 处的横向相位梯度约为 l/r越靠近中心梯度越大。要分辨这个梯度每个像素上的相位变化不能超过 π否则相位包裹线会断裂强度环也会出现锯齿。实际工程判断可以用一个简单经验网格间距 h L / N要求 w0 / h |l|。比如 N512、L5mm、w01mm 时 h≈9.8μmw0/h≈102那么 |l| 取到 4~5 都安全。想快速确认当前参数下能支持多大拓扑荷可以运行下面两行h L / N; fprintf(w0/h %.1f建议 |l| %d\n, w0/h, floor(w0/h));这个经验式没有严格推导应付课程设计和科研预研已经足够。当 l 超过该范围时常见做法不是盲目加大 N而是先缩小 L把计算区域收窄到 5w0 以内再升 N 到 1024。2.4 常见第一个误用把强度图当成涡旋判据很多人第一次运行会拿imagesc(abs(LG))看如果看到中心不黑就认为代码错了。实际上abs取模得到的是振幅要平方才是光强。在纯相位模型中没有 (√2 r/w0)^|l| 这个因子中心区域仍然有非零强度看起来像亮斑这是模型假设而不是 bug。判断有没有产生涡旋核心证据是相位图中心存在 2πl 的相位缠绕和位错。最直接的做法是输出angle(LG)并绕中心走一圈观察相位是否递增或递减 2πl。3. 涡旋光束程序主干相位屏、高斯载波与叉形干涉3.1 相位屏的作用把拓扑荷变成可加工的数字掩膜SLM 或相位掩膜的灰度对应相位延迟量在 MATLAB 里就是一个二维数组值域在 [-π, π] 或 [0, 2π)。相位屏的质量直接决定后面干涉和传播的可靠性所以值得单独拿出来讲。生成相位屏最直接的方式是wrapToPi(l * Theta)但wrapToPi属于 Mapping Toolbox不是所有廉价的 MATLAB 版本都有。更通用的写法是mod(l * Theta pi, 2*pi) - pi效果完全一致。l1 时相位图有一条割线l2 有两条这是相位包裹的自然结果不是程序错误。3.2 用相位屏混合高斯载波生成可观测的叉形干涉图干涉图是验证涡旋光的标准手段也是 vortex beam.rar 里最常见的输出。下面的代码生成一个 l2 的叉形干涉图。h L / N; w0 2e-3; % 载波束腰比相位屏的束腰稍宽 E_gauss exp(-R.^2 / w0^2); % 高斯载波 phase_mask mod(l * Theta pi, 2*pi) - pi; E_vortex E_gauss .* exp(1i * phase_mask); % 倾斜参考光形成横向条纹 tilt_freq 1 / (16 * h); % 参考光空间频率周期16像素 E_ref exp(1i * 2 * pi * tilt_freq * X); I_dual abs(E_vortex E_ref).^2; figure imagesc(x*1e3, x*1e3, I_dual); axis image; colormap gray; title(叉形干涉图);exp(1i*phase_mask)与直接写exp(1i*l*Theta)等价但走相位屏这条路径后后续可以替换成实测畸变波前或随机湍流相位更贴近真实实验。tilt_freq用像素周期 16 表示条纹密度适中叉形分叉清晰可数。分叉方向能判断拓扑荷符号l0 且参考光向右倾斜时叉口朝上l0 时叉口朝下。3.3 涡旋光束程序的四个关键参数速查表下面的参数表覆盖了绝大多数涡旋光束 MATLAB 程序需要反复调整的变量。参数符号典型值主要影响网格点数N512 或 1024相位分辨率和频域格子密度计算域边长L5w0 ~ 10w0边缘能量泄漏和频域采样间隔束腰半径w01e-3 ~ 2e-3暗核尺寸和远场环半径拓扑荷l-3 ~ 3叉形数量、OAM 谱峰值位置一个常用策略是先把 L 固定到 5w0调好 N 和 w0最后再把传播距离放进来。不要同时改四个参数否则一旦出现异常你很难判断是哪个变量引起的。传播距离看瑞利距离 zR π w0² / λw01mm、λ632.8nm 时 zR 约 4.96m做桌面实验距离时一般都远小于 zR强度环半径几乎不随距离变化。3.4 程序包里常见的中心偏移问题很多网上版本用linspace(0, L, N)建坐标系涡旋奇点落到网格角点而不是中心像素结果强度环看起来像缺了一角。用x linspace(-L/2, L/2, N)、N 取偶数时0 不会落在任何像素上奇点被四个像素平分表现是暗核中心出现一个微小亮十字。我一般会单独处理如果以像素为准就采用(-N/2:N/2-1) * h如果以物理坐标为准则用(-N/20.5:N/2-0.5) * h做半个像素偏移。关键是相位屏、高斯载波和后面的传播函数必须统一使用同一套坐标定义混用坐标系是很多 rar 包里传播结果不对称的真正原因。验证奇点位置的办法打印[~, idx] min(abs(x));得到中心列索引再观察相位图中心四个像素是否有逆时针或顺时针的螺旋次序。如果强度环沿 x 或 y 方向不对称优先检查坐标定义。4. 传播仿真与拓扑荷检测角谱法、叉数判据与 OAM 谱分解4.1 角谱法传播涡旋光束程序里最核心的频域循环涡旋光的横向分辨和 OAM 谱不仅取决于初始面分辨率还取决于传播算子是否写对。角谱法是做涡旋光传播最常见的方案公式是一个傅里叶变换乘频域传播因子再做逆变换实现门槛低物理意义直观。zD 0.3; % 传播距离 0.3m k 2 * pi / lambda; fx (-N/2:N/2-1) * (1/L); % 频域坐标 [FX, FY] meshgrid(fx, fx); H_AS zeros(N); % 预分配 idx_valid (FX.^2 FY.^2) (1/lambda)^2; % 排除倏逝波 H_AS(idx_valid) exp(1i * k * zD * sqrt(1 - lambda^2 * (FX(idx_valid).^2 FY(idx_valid).^2))); Ef fft2(fftshift(E_vortex)); % 空间域原点移到数组中心再做变换 E_z ifft2(Ef .* ifftshift(H_AS));频域格点间隔是 1/L原点在数组中心所以坐标用(-N/2:N/2-1)这种对称形式。平方根里小于零的频域分量对应倏逝波必须清零否则逆变换会把指数增长带回来。fftshift和ifftshift的顺序不能反写反的表现是传播后的场整体旋转 180° 或者出现中心镜像。传播距离小于 0.1 zR 时强度环半径几乎不变这是检查代码是否正确的一道天然测试。如果要在循环里做多步传播把上面三步封装成propagate(E, zD, lambda, L)这样的函数每次只传复振幅和距离避免重复分配频域变量。4.2 传播后看什么暗环半径变化与模式保持传播后最先检查的是中心暗核是否仍在。涡旋的相位奇点会在传播过程中保持但数值色散会让暗核变形。把abs(E_z).^2画出来过中心取横截面看最小值是否低于周围峰值的 10%。如果中心已经是局部最大值说明采样不足或频域掩膜被截断太狠。另一个检查项是暗环半径。近场范围内暗环半径约等于 w0 sqrt(|l|/2)远场时随传播距离近似线性扩散。如果看到暗环收缩或者出现多个同心环大概率是初始场里混入了高阶径向模式或者频域网格间隔设置过大导致高频成分被截断。4.3 用 OAM 谱分解把涡旋光“称重”干涉图能看符号OAM 谱才能看纯度。做法是把接收到的复振幅在一个环形掩膜上投影到 exp(-imθ)得到各个拓扑荷分量的能量占比。这段代码是涡旋光程序包里最重要的验收工具。m_vals -5:5; spec zeros(size(m_vals)); annulus (R 0.5*w0) (R 1.5*w0); % 环形采样带 for i 1:numel(m_vals) m m_vals(i); proj sum(E_z(annulus) .* exp(-1i * m * Theta(annulus))); spec(i) abs(proj)^2; end bar(m_vals, spec / max(spec)); xlabel(拓扑荷 m); ylabel(归一化投影能量);投影必须限定在环带上不能全平面做。全平面投影会因为高斯振幅分布的径向形状把一部份能量泄漏到邻近拓扑荷谱变得又宽又平。环带内外半径取 0.5w0 到 1.5w0 是最常见的起始值传播后的光束展宽了可以把上限放宽到 2.5w0。如果谱里出现明显的主峰之外的杂散峰先检查初始相位屏的割线是否落在某个 m 分量上。4.4 一个常被忽略的符号问题正负拓扑荷与镜像MATLAB 的cart2pol在默认坐标下角度 θ 从 x 轴正向开始逆时针增长。把 l 改成 -l强度图不变相位图和干涉图整体翻转。如果 rar 包里的干涉图叉形方向和你自己实验的语义反了常见原因是有人把atan2(Y, X)的坐标写成了atan2(X, Y)等价于做了一次镜像。处理这个问题的建议是始终在代码里显式写[Theta, R] cart2pol(X, Y);然后找一个已知符号的 l 生成干涉图和实验图对比叉口方向。方向对上号以后再开始改其他参数。5. 拿到涡旋光束程序包后先改的三个参数与控制台自检5.1 先改这三个参数不要动模型解压 vortex beam.rar 后先不碰物理结构只调三个量就能解决大部分异常。参数改法对应现象N512 → 1024相位图糊成一团干涉分叉不清晰L5e-3 → 8e-3远场强度环被计算边界吞掉w01e-3 → 2e-3暗核像素太少看不到完整光环这三个参数相互耦合建议一次只改一个每次改动后重新跑相位图和干涉图观察变化趋势是否符合预期。像素糊掉先升 N环被切断先放大 L暗核看不清再加大 w0。顺序反了你会同时遇到两个方向的异常很难定位。5.2 控制台三行快速验收代码下面三行能在不改动其他代码的情况下确认程序是否真的在算涡旋光。assert(~isreal(E_vortex), 复振幅被拍成实数矩阵了); phase_ramp diff(unwrap(angle(E_vortex(N/21, :)))); fprintf(横向相位变化范围 %.3f rad\n, max(abs(phase_ramp))); [~, pk] max(spec); fprintf(OAM 谱主峰 m %d\n, m_vals(pk));第一条检查在最常见的问题解压后用load读取数据时变量已经被保存成强度图而不是复振幅。第二条看横向相位是不是近似线性如果是一条乱跳的线说明坐标定义错了。第三条确认OAM 谱峰值落位正确这是整个涡旋光束程序能继续往下用的底线。5.3 固定一组基线参数把输出保存成 PNG项目答辩或者传包给别人之前我一般会把 l2、zD0.3m、N512 这组参数固定住输出强度图、相位图、叉形干涉图和 OAM 谱四张图。之后换机器、换 MATLAB 版本只要重跑一遍脚本、对比这四张图的特征就能快速判断计算结果是否退化。对采样和坐标的所有改动最终都会收敛到这份基线图的形态上。本文还有配套的精品资源点击获取
