简介本资源是一套面向光学工程初学者与科研人员的MATLAB实践工具包聚焦涡旋光束生成原理与超表面相位调控建模解决轨道角动量光场设计、相位分布仿真及纳米光学器件模拟等核心问题。压缩包共10个文件9个.m脚本1个Word文档总大小仅18KB轻量高效.m文件涵盖涡旋相位模板构建含不同拓扑荷l设置、傅里叶域传播模拟、超表面等效相位响应建模及衍射图样可视化等关键环节配套docx文档系统梳理了相位变化物理机制、涡旋光数学表达与MATLAB实现逻辑。已有880人学习下载资源结构清晰、代码模块化强每段脚本均对应明确光学功能如zone plate.m实现菲涅耳波带片调制Untitled系列文件分别演示基础涡旋、叠加态及超表面耦合效应便于读者逐层理解、调试验证并快速迁移至实际超表面设计任务。 从“用MATLAB看涡旋光的相位变化”说起这半年陆续帮几个课题组看过类似的问题发现大家的困惑基本集中在三个地方——涡旋光相位不会写、超表面相位不会叠加、传播之后的相位不会分析。与其一次一次回复不如把整套流程完整拆开讲一遍。这篇主要聊的是怎么用MATLAB构建一束涡旋光怎么在“超表面”这个薄层上施加相位变化以及怎么把相位变化的结果用直观的图呈现出来。内容适合刚接触涡旋光或者超表面仿真的人也适合那些已经在用FDTD/CST做单元结构、想先用MATLAB快速验证系统级相位分布的人。整套代码以相位计算为核心不依赖商用光电仿真软件也不需要额外工具箱纯手写就可以跑起来。1. 先搞清楚涡旋光、相位变化、超表面这三者的关系1.1 涡旋光到底“旋”在哪里涡旋光的核心不是强度分布多好看而是它携带一个螺旋波前。普通的平面波波前是平面传播方向就是波矢方向涡旋光的波前则是一个螺旋面光在传播时相位绕着传播轴持续旋转中心位置因为相位不确定所以强度天然为零这就是常说的“暗核”。在数学上涡旋光的相位可以写成φ l·θ这里的l是拓扑荷θ是方位角。拓扑荷取整数时相位绕中心一圈变化2πl。l为正号说明螺旋方向是逆时针负号则是顺时针。这个螺旋相位会带来轨道角动量OAM这也是涡旋光这几年在光通信、微粒操控、超分辨成像里被反复研究的原因。用MATLAB写这个相位很简单核心就是atan2(Y,X)[X, Y] meshgrid(x, x); theta atan2(Y, X); l 1; phase_vortex l * theta;这里atan2得到的就是从-π到π的方位角乘上拓扑荷之后就得到完整的螺旋相位分布。要注意l2时相位变化是4π肯定会在显示图上看到两圈颜色跳变这不是错误是正常的相位缠绕。1.2 超表面在这里充当什么角色超表面metasurface简单理解就是一层亚波长厚度的“相位画笔”上面排列着很多亚波长的纳米结构单元每个单元都可以独立改变入射光的相位。传统光学用透镜、棱镜的厚度差延迟光程超表面则靠单元的几何尺寸、旋转角度或材料参数在亚波长尺度内改变局部相位。在我们的仿真语境里超表面通常被抽象成一个纯相位屏即在某个z位置上给光场乘上一个复振幅exp(iφ_meta(x,y))。这种“薄屏近似”在绝大多数系统级分析中足够用因为它关注的是相位分布的整体行为而不是单个纳米柱的散射细节。相位变化这个关键词在这里就有两个层面的意思。第一层是超表面自身产生的相位型分布比如螺旋相位、梯度相位、聚焦相位第二层是涡旋光透过超表面之后总相位场发生的变化。两者叠加本质上就是两个相位面的干涉与重构。1.3 为什么这套仿真值得用MATLAB做市面上能做光场传播的软件不少比如VirtualLab、Zemax、Lumerical但它们的学习成本和授权成本都不低。用MATLAB的好处在于相位构建、传播计算、图像分析、参数扫描可以在同一个环境里完成而且每一步都能直观看到中间结果。尤其是当你需要快速验证“超表面排布方式是否能把高斯光转成涡旋光”“不同拓扑荷叠加后相位长什么样”“传播一段距离后相位演化是否稳定”这类问题时MATLAB写个脚本通常几分钟就能出图比先在CAD里建模再丢进FDTD仿真要快几个量级。当然FDTD那种全波仿真是另一回事它是用来做单元级验证的MATLAB和FDTD的正确关系是“先系统级快速验证再单元级精确仿真”。这也就是标题里“phasechange”的准确价值相位变化不是一个抽象概念而是可以用几十行代码定量计算、可视化、分析的具体量。2. 相位构建的数值细节网格、坐标、采样精度2.1 网格划分是相位仿真的地基相位仿真最容易被忽略的就是网格设置。网格的物理尺寸和采样点数直接决定后续所有计算的有效性。推荐这样设置参数lambda 632.8e-9; % He-Ne激光波长 L 4e-3; % 仿真区域边长 4mm N 1024; % 网格点数 dx L / N; % 空间采样间隔 x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x);这里有几个原则值得说明。仿真区域边长L要远大于光束束腰否则光场会被边界截断采样点数N取2的幂1024、2048那种这样不管是用FFT还是角谱传播都方便。采样间隔dx必须满足奈奎斯特采样条件特别是超表面相位梯度很陡的时候如果dx过大高频相位变化会被直接抹平。我见过最典型的问题就是有人把L设得过大比如10mm但N只有256结果一个涡旋相位圆环上就十几个像素画出来的相位图跟马赛克一样中心暗核根本看不出来。经验值是束腰约1mm时仿真区域4mm配1024点是比较稳的组合最小特征尺寸大约4微米足够分辨绝大多数亚波长结构相位梯度的宏观效果。2.2 构建标准的拉盖尔-高斯涡旋光束很多教程会说“涡旋光就是exp(ilθ)”这只是相位因子。光束还得有强度包络否则直接用平面波乘exp(ilθ)也能得到涡旋相位但强度分布太均匀不贴近实际实验。最常用的是拉盖尔-高斯模式LG模式其复振幅化简到二维截面可以写成w0 1e-3; % 束腰半径 r sqrt(X.^2 Y.^2); theta atan2(Y, X); l 2; % 拓扑荷 E_lg (sqrt(2)*r/w0).^abs(l) ... .* exp(-r.^2/w0^2) ... .* exp(1i*l*theta);这个公式里的(sqrt(2)*r/w0).^abs(l)决定了高阶涡旋的强度环会随着拓扑荷增大而向外扩展而exp(-r.^2/w0^2)负责把整体能量约束在束腰范围内。实际用的时候还可以加上曲率因子和Gouy相位但对相位变化分析来说上面这个简化版本已经足够。为什么要强调用LG模式而不是纯exp(ilθ)因为纯相位项在r趋近零的地方不存在强度为零的约束数值上可能因为相位不连续导致暗核不够暗。LG模式在r0处自然为零符合真实涡旋光的物理特性。2.3 超表面相位屏的三种常见写法超表面的相位设计按功能可以分成几类在MATLAB里都需要写成二维相位矩阵。第一种是螺旋相位用于把高斯光转换成涡旋光或改变涡旋光拓扑荷l_meta 3; phase_spiral l_meta * atan2(Y, X);第二种是梯度相位让光束发生偏转对应超表面的“广义斯涅尔定律”常见于分束器设计period 500e-6; % 超表面单元周期对应的宏观相位周期 phase_gradient 2*pi/period * X;第三种是聚焦相位相当于超表面透镜用于汇聚光场f 10e-3; % 焦距 phase_lens -k * (X.^2 Y.^2) / (2*f);实际项目中超表面上的相位分布往往是几种基本相位的复合。比如既要偏转又要涡旋就写成θ_meta l_meta·θ 2π/Λ·x。这种线性叠加在物理上是合理的因为超表面每个单元只能给出一个确定的相位延迟这个相位延迟可以同时满足多个功能需求只要空间采样足够。要注意的是写成相位矩阵之后真正作用在光场上时要用exp(1i * phase)构造复振幅调制项。这一步别写成phase .* E那是错误的相位必须放在指数里。3. 把相位变化动态演出来传播与可视化3.1 相位叠加的核心操作入射涡旋光经过超表面后的出射光场在薄屏近似下就是两串复振幅逐点相乘E_out E_lg .* exp(1i * phase_meta);这行代码是整篇文章里最核心的东西。它的物理含义是超表面上每个位置对入射光的相位做一次局部修改透射光的等相位面因此被重新塑形。如果phase_meta本身是螺旋相位就会改变入射涡旋光的拓扑荷如果phase_meta是梯度相位就会给整个光束叠加一个横向波矢偏折如果phase_meta是透镜相位就会改变波前曲率。相位叠加之后最好在频域里看一眼角谱这样可以快速判断出射光是否混入了不需要的高频分量。尤其是梯度相位叠加涡旋相位时频谱中心会往一侧偏移偏折角度跟梯度周期成反比。3.2 用角谱法传播到任意距离只看零距离处的相位分布是不够的实际系统里探测器或者下一级光学元件都离超表面有一段距离。这里推荐用角谱法Angular Spectrum Method做光场传播因为它物理意义清晰MATLAB实现也方便。角谱法的思路是把空间光场分解成一系列平面波每个平面波传播一段距离后在频域乘上一个传播相位因子再逆变换回空间域。z 100e-3; % 传播距离 10cm k 2*pi / lambda; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); H exp(1i * k * z .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); E_far ifft2(ifftshift(fftshift(fft2(E_out)) .* H));这里面有几个容易出错的地方单独说一下。第一频域坐标的构造必须与空间域匹配。fftshift(fft2(E_out))之后零频在矩阵中心频域坐标也必须从负值到正值排列-N/2 : N/2-1而且频率间隔是1/L其中L是空间域仿真区域的边长。第二sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)里面如果出现负数表示倏逝波可以直接令其为零否则数值上会指数发散。虽然绝大多数传播距离下这个项不会碰到但当你把z设得极小或者网格设得很大时一定要加一个判断mask (1 - (lambda*FX).^2 - (lambda*FY).^2) 0; H zeros(N); H(mask) exp(1i * k * z .* sqrt(1 - (lambda*FX(mask)).^2 - (lambda*FY(mask)).^2));第三角谱法的适用距离很宽近场远场都行但注意它要求传播路径上是均匀介质不能穿过透镜或其他折射率突变界面。要模拟经过透镜后再传播就得拆成三段透镜前传播、透镜相位作用、透镜后传播。3.3 相位图、强度图、截面图组合出效果相位图是判断涡旋光质量最直观的方式。MATLAB里可以用imagesc显示相位figure; imagesc(x*1e3, x*1e3, angle(E_far)); axis xy; axis square; colormap(hsv); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(远场相位分布);这里用hsv色图而不是默认parula因为hsv在±π处是同一颜色首尾衔接特别适合显示周期性相位。如果用parula你会看到一条明显的颜色断层那不是物理问题是色图不匹配造成的视觉假象。强度图用imagesc(abs(E_far).^2)显示为了看清暗核通常会把强度做归一化再取对数或者直接叠加一个对比度增强I abs(E_far).^2; I I / max(I(:)); imagesc(x*1e3, x*1e3, log10(I 1e-4));相位变化不光要看二维分布还可以提取沿某个圆的相位值来验证拓扑荷。比如在以光轴为中心、半径r0处采样一圈相位unwrap之后如果线性变化2πl就说明拓扑荷确实是l。这个验证方法实验和仿真都通用。r0 0.5e-3; % 采样子圆半径 th_sample linspace(0, 2*pi, 360); x_sample r0*cos(th_sample); y_sample r0*sin(th_sample); phase_sample interp2(X, Y, angle(E_far), x_sample, y_sample, spline); phase_unwrap unwrap(phase_sample);如果你看到unwrap之后相位曲线不是单调上升或下降的斜线而是来回起伏说明这一圈相位不纯净可能是叠加了别的像差或者涡旋中心没对准网格原点。4. 实操过程一整套涡旋光超表面相位变化仿真4.1 参数初始化我习惯把所有可能需要调整的量集中放在脚本最前面方便扫描参数。这里给出一份可直接跑通的基础配置。clear; clc; close all; % 光场参数 lambda 632.8e-9; % 波长 (m) w0 1e-3; % 束腰半径 (m) l_inc 1; % 入射涡旋拓扑荷 % 超表面参数 l_meta 2; % 超表面叠加的拓扑荷 period_g 1.2e-3; % 梯度相位周期 (m)可选 f_meta 20e-3; % 超表面透镜焦距 (m)可选 % 仿真参数 L 4e-3; % 仿真区域边长 (m) N 1024; % 网格点数 z_prop 50e-3; % 传播距离 (m) dx L / N; x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); r sqrt(X.^2 Y.^2); theta atan2(Y, X); k 2*pi / lambda;把拓扑荷、焦距、周期都写在顶部的好处是后面想扫一个参数曲线时直接套一层for循环就行不用改核心代码。4.2 构建光场与超表面相位入射光场用LG模式超表面相位做一个“涡旋偏转”的复合设计这在很多光束操纵场景里是实用组合。% 入射涡旋光场 E_inc (sqrt(2)*r/w0).^abs(l_inc) .* exp(-r.^2/w0^2) .* exp(1i * l_inc * theta); % 超表面复合相位 phase_meta l_meta * theta 2*pi/period_g * X; % 出射光场 E_out E_inc .* exp(1i * phase_meta);复合相位的意思是这个超表面同时做两件事把入射的l1涡旋变成l3涡旋并让整个光束朝x方向偏转一个角度。偏转角理论上等于arcsin(λ/period_g)代入上面的参数大约是0.03度偏转量很小但在远场图上可以看出来。如果只关心纯涡旋变化把梯度项去掉就行。透镜聚焦项也一样需要就加不需要就删整个脚本是模块化的。4.3 传播与相位提取传播用角谱法按第三节的代码封装成一个函数会更方便。function E_prop asm_propagate(E_in, dx, lambda, z) N size(E_in, 1); L N * dx; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); k 2*pi / lambda; H exp(1i * k * z .* sqrt(max(0, 1 - (lambda*FX).^2 - (lambda*FY).^2))); E_prop ifft2(ifftshift(fftshift(fft2(E_in)) .* H)); end调用方式很简单E_far asm_propagate(E_out, dx, lambda, z_prop);然后用前面提到的方法分别画入射相位、超表面相位、出射相位、传播后相位四张图排成2×2一眼就能看出相位变化链条。4.4 一组完整结果判读画出图之后建议按下面顺序检查结果。第一看入射光强度图中心应该有暗核暗核尺寸大约在波长量级周围是亮环。第二看出射光的相位分布如果超表面只加螺旋相位那么相位图绕中心一圈应该有l_incl_meta个2π跳变。如果加了梯度项相位图上会叠加一系列平行斜条纹这是偏转相位的表现。第三看传播后强度图涡旋中心依然保持暗核说明相位奇点稳定。如果暗核被填充大概率是角谱法里没做倏逝波滤波或者网格采样不够导致高频信息混叠。第四看提取的圆周相位unwrap之后斜率对应的拓扑荷数应该等于l_incl_meta。这是量化验证比肉眼看图靠谱得多。5. 常见问题与排查技巧实录5.1 相位图看起来一团乱全是细条纹出现这种情况通常是相位值的2π跳变太密导致imagesc的色图断层密集。在显示相位之前可以用wrapToPi对相位做缠绕处理或者直接用angle(exp(1i*phase))把相位压到[-π, π]区间。如果是显示传播后的相位由于光束整体有一个很强的平面波因子相位随空间快速振荡这种情况不适合直接看总相位应该先减去参考平面波相位比如k*z或者k*x*sinθ只看相对相位变化。5.2 中心暗核消失涡旋变成实心亮斑从仿真角度说原因八成是网格中心点没有精确落在坐标原点上。x (-N/2 : N/2-1) * dx在N为偶数时坐标点并不包含零导致涡旋中心落在四个像素的交叉点上相位奇点在离散网格中无法正确表达。解决办法是让N为奇数或者把坐标稍微平移半个像素x (-(N-1)/2 : (N-1)/2) * dx;这样中心像素就是0暗核不会被抹平。另一个可能原因是入射光场里加了过大的束腰w0导致暗核半径小于采样间隔中心看起来就是亮的。可以把w0调小或者把仿真区域L调大让暗核占据5个像素以上。5.3 角谱传播后能量明显减少角谱法本身是能量守恒的在无损耗条件下能量减少通常是因为边界截断。光场传播时会在空间上逐渐扩展如果仿真区域L不够大光场边缘会碰到边界FFT的周期性会把泄漏过来的场叠加到另一边造成数值误差。解决方法是把仿真区域设为光束尺寸的3到5倍或者像实际仿真里常用的做法一样在初始光场上乘一个超高斯窗口来抑制边缘衍射。还有一种情况是传播距离太长远场光斑直径超过L这种时候应该减小z_prop或者改用傅里叶变换法Fraunhofer衍射直接计算远场分布。5.4 拓扑荷提取不准unwrap后斜率不是整数这个问题的根源通常是采样本圆穿过涡旋中心导致相位跳变不连续。采样半径不能选太大也不能选太小太大会进入强度很低的区域相位噪声大太小会切割到暗核区域相位奇点附近相位本身就不连续。经验做法是把采样半径设为亮环的峰值位置大概在LG模式的半径处或者选束腰的0.7到1.0倍。另外interp2默认的线性插值在相位跳变边缘会出问题建议用spline插值并且先对复振幅插值再取相位角而不是直接对相位矩阵插值E_complex_sample interp2(X, Y, E_far, x_sample, y_sample, spline); phase_sample angle(E_complex_sample);这样处理之后unwrap的稳定性会好很多。5.5 复合相位扫描时结果不理想如果你在做参数扫描比如连续改变l_meta要注意绕回问题。相位从π到-π的跳变在数值上表现为一个很大的负向跳变unwrap可以解决但在做乘法叠加时可能会引入意想不到的伪影。建议所有相位操作都用复指数形式进行也就是把相位一直放在exp(1i*φ)里面等到最后需要显示相位时才取angle()。5.6 仿真配置问题速查表现象可能原因排查与解决相位图条纹过密相位周期性太强或显示方式不对使用wrapToPi压缩到[-π,π]或显示相对相位中心暗核消失N为偶数导致中心不在网格点改用奇数N或加半像素平移传播后能量减少边界截断或z过长增大L或改Fraunhofer远场法拓扑荷提取非整数采样半径不合适或插值方法不当在强度峰值半径处采样对复振幅插值再取相位低频背景太强入射高斯分量未完全转化为涡旋分量检查超表面相位幅值是否精确为2π整数倍周期图中有很强网格状条纹FFT采样混叠提高N或增大L检查频域坐标是否与空间域匹配6. 超表面相位设计的几个实用心得6.1 相位量化不一定非要8级有的同学一开始设计超表面单元就会陷入“8台阶量化”的思维定式因为传统的二元光学元件确实需要量化台阶。但这里要分清楚台阶量化是加工工艺限制不是物理必需。在MATLAB初步验证时完完全全可以先用连续相位做仿真确认系统级效果没问题再量化。如果量化后效果变差再考虑用一些优化算法来补偿量化误差。我实际做过一个对比连续相位、4台阶、8台阶、16台阶当单元的相位调控范围小于2π时比如只调控0.6π4台阶和8台阶的结果几乎没有差别。量化的主要损失发生在相位梯度很陡的区域那里相邻像素的相位跳变大量化台阶把本来就稀疏的采样信息又砍了一刀。6.2 孔径效应对相位变化的影响超表面不可能无限大但它有边界。这种有限孔径会给涡旋光引入衍射效应反映在远场就是强度环展宽、暗核边缘弥散。仿真的时候想模拟真实的有限孔径超表面只要在超表面相位上乘一个圆形孔径函数aperture_r 1.5e-3; aperture_mask double(r aperture_r); E_out E_inc .* exp(1i * phase_meta) .* aperture_mask;孔径边缘会导致高频衍射如果传播距离较长你会发现强度分布会在孔径边缘出现菲涅尔环纹这是正常现象。想要压低这个效应可以把硬边改成软边超高斯渐变但物理上实验加工出来的超表面基本都是硬边所以这个衍射效应其实也限制着实际系统的性能。6.3 别忘了偏振自由度涡旋光与超表面的结合很多时候还会涉及偏振。最典型的就是几何相位Pancharatnam-Berry相位超表面它的相位调控依赖单元的旋转角入射圆偏振光经过后转换为反向圆偏振光同时携带2倍旋转角的相位延迟。如果你研究的超表面是基于这类机制那么MATLAB仿真里就不能只用标量相位要考虑琼斯矩阵。做偏振相关仿真时推荐用琼斯矢量形式E_in_xy [E_x; E_y]; % 2xN矩阵两个偏振分量 J_meta [cos(alpha), -sin(alpha); sin(alpha), cos(alpha)]; % 单元旋转计算过程会比标量形式稍微复杂但思路一致每个像素位置用单元旋转角α生成一个局域琼斯矩阵作用在入射偏振态上然后分别传播两个正交分量。这样既能分析强度也能分析偏振转换效率。6.4 从远场反推超表面相位有一个前提很多人想做逆向设计给定想要的远场涡旋光反推超表面的相位分布。这个概念没错但有一个前提——超表面只能提供一个2π范围内的相位且这个相位对波前的作用是局部的。简单说它只能改“相位面”的形状不能凭空改变光场的振幅分布。所以逆推的时候如果目标光场的振幅分布和入射光差异很大单靠纯相位超表面做不好这时候要么加入振幅调控这会牺牲效率要么接受部分能量转换不理想。这个限制在做涡旋光生成时尤其明显高斯光转涡旋光效率天然受制于模式匹配。6.5 MATLAB仿真流程图在脑中的正确顺序把整套仿真浓缩成一句话光源构建、超表面相位叠加、角谱传播、相位/强度分析。这个顺序千万别乱。如果你先传播再叠加超表面相位等于把超表面放在传播路径末端物理场景完全不同。7. 从相位仿真到实际器件的衔接建议仿真不等于实际两者之间还有一大截路要走。我建议做完MATLAB系统级仿真之后再按这个顺序往下推进首先做单元级全波仿真CST或Lumerical FDTD验证你设计的超表面单元确实能实现所需相位的2π覆盖。单元结构的具体尺寸和材料参数在这个阶段确定仿真周期一般为亚波长微米量级。其次把单元仿真得到的相位-尺寸/旋转角关系表导出映射回MATLAB。这个映射很关键因为单元仿真出来的相位响应通常不是线性的直接按理想相位排布单元会导致实际相位误差产生明显的散射噪声。然后在MATLAB里做一个“真实相位”的二次验证用映射后的实际相位矩阵替换理想的phase_meta再次执行传播仿真。这一步能暴露很多单元设计问题比如相位覆盖不足、相邻单元耦合导致相位偏差等。最后才考虑加工和光学测量。到实验阶段你手上已经有一套“理想仿真→真实参数仿真”的对照数据测试结果出来也知道该怎么分析。这个流程很稳能避免那种“软件仿真看着完美加工出来完全不能用”的尴尬情况。我个人在实际操作中的体会是MATLAB最容易被低估的价值不是计算本身而是它让你在不到半小时内把整套物理逻辑过了一遍。相位构建、叠加、传播、分析每一步都在逼着你回答“这一步的物理意义是什么”等这些概念在脑子里彻底清晰了再上手任何专业商用软件都很快。分享一个小技巧收尾在你第一次跑通仿真后把入射拓扑荷改成0、把超表面改成纯螺旋相位对比一下出射光和入射光的相位分布你会在图上直观看到“一个平面波如何被超表面卷成涡旋”的整个过程。这个对比图特别适合用来检查自己的仿真链路是否正确也非常适合写论文时放给审稿人看——一张图解释清楚相位变化的核心机理。本文还有配套的精品资源点击获取 MATLAB涡旋光与超表面相位变化仿真：从相位构建到角谱传播
简介本资源是一套面向光学工程初学者与科研人员的MATLAB实践工具包聚焦涡旋光束生成原理与超表面相位调控建模解决轨道角动量光场设计、相位分布仿真及纳米光学器件模拟等核心问题。压缩包共10个文件9个.m脚本1个Word文档总大小仅18KB轻量高效.m文件涵盖涡旋相位模板构建含不同拓扑荷l设置、傅里叶域传播模拟、超表面等效相位响应建模及衍射图样可视化等关键环节配套docx文档系统梳理了相位变化物理机制、涡旋光数学表达与MATLAB实现逻辑。已有880人学习下载资源结构清晰、代码模块化强每段脚本均对应明确光学功能如zone plate.m实现菲涅耳波带片调制Untitled系列文件分别演示基础涡旋、叠加态及超表面耦合效应便于读者逐层理解、调试验证并快速迁移至实际超表面设计任务。 从“用MATLAB看涡旋光的相位变化”说起这半年陆续帮几个课题组看过类似的问题发现大家的困惑基本集中在三个地方——涡旋光相位不会写、超表面相位不会叠加、传播之后的相位不会分析。与其一次一次回复不如把整套流程完整拆开讲一遍。这篇主要聊的是怎么用MATLAB构建一束涡旋光怎么在“超表面”这个薄层上施加相位变化以及怎么把相位变化的结果用直观的图呈现出来。内容适合刚接触涡旋光或者超表面仿真的人也适合那些已经在用FDTD/CST做单元结构、想先用MATLAB快速验证系统级相位分布的人。整套代码以相位计算为核心不依赖商用光电仿真软件也不需要额外工具箱纯手写就可以跑起来。1. 先搞清楚涡旋光、相位变化、超表面这三者的关系1.1 涡旋光到底“旋”在哪里涡旋光的核心不是强度分布多好看而是它携带一个螺旋波前。普通的平面波波前是平面传播方向就是波矢方向涡旋光的波前则是一个螺旋面光在传播时相位绕着传播轴持续旋转中心位置因为相位不确定所以强度天然为零这就是常说的“暗核”。在数学上涡旋光的相位可以写成φ l·θ这里的l是拓扑荷θ是方位角。拓扑荷取整数时相位绕中心一圈变化2πl。l为正号说明螺旋方向是逆时针负号则是顺时针。这个螺旋相位会带来轨道角动量OAM这也是涡旋光这几年在光通信、微粒操控、超分辨成像里被反复研究的原因。用MATLAB写这个相位很简单核心就是atan2(Y,X)[X, Y] meshgrid(x, x); theta atan2(Y, X); l 1; phase_vortex l * theta;这里atan2得到的就是从-π到π的方位角乘上拓扑荷之后就得到完整的螺旋相位分布。要注意l2时相位变化是4π肯定会在显示图上看到两圈颜色跳变这不是错误是正常的相位缠绕。1.2 超表面在这里充当什么角色超表面metasurface简单理解就是一层亚波长厚度的“相位画笔”上面排列着很多亚波长的纳米结构单元每个单元都可以独立改变入射光的相位。传统光学用透镜、棱镜的厚度差延迟光程超表面则靠单元的几何尺寸、旋转角度或材料参数在亚波长尺度内改变局部相位。在我们的仿真语境里超表面通常被抽象成一个纯相位屏即在某个z位置上给光场乘上一个复振幅exp(iφ_meta(x,y))。这种“薄屏近似”在绝大多数系统级分析中足够用因为它关注的是相位分布的整体行为而不是单个纳米柱的散射细节。相位变化这个关键词在这里就有两个层面的意思。第一层是超表面自身产生的相位型分布比如螺旋相位、梯度相位、聚焦相位第二层是涡旋光透过超表面之后总相位场发生的变化。两者叠加本质上就是两个相位面的干涉与重构。1.3 为什么这套仿真值得用MATLAB做市面上能做光场传播的软件不少比如VirtualLab、Zemax、Lumerical但它们的学习成本和授权成本都不低。用MATLAB的好处在于相位构建、传播计算、图像分析、参数扫描可以在同一个环境里完成而且每一步都能直观看到中间结果。尤其是当你需要快速验证“超表面排布方式是否能把高斯光转成涡旋光”“不同拓扑荷叠加后相位长什么样”“传播一段距离后相位演化是否稳定”这类问题时MATLAB写个脚本通常几分钟就能出图比先在CAD里建模再丢进FDTD仿真要快几个量级。当然FDTD那种全波仿真是另一回事它是用来做单元级验证的MATLAB和FDTD的正确关系是“先系统级快速验证再单元级精确仿真”。这也就是标题里“phasechange”的准确价值相位变化不是一个抽象概念而是可以用几十行代码定量计算、可视化、分析的具体量。2. 相位构建的数值细节网格、坐标、采样精度2.1 网格划分是相位仿真的地基相位仿真最容易被忽略的就是网格设置。网格的物理尺寸和采样点数直接决定后续所有计算的有效性。推荐这样设置参数lambda 632.8e-9; % He-Ne激光波长 L 4e-3; % 仿真区域边长 4mm N 1024; % 网格点数 dx L / N; % 空间采样间隔 x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x);这里有几个原则值得说明。仿真区域边长L要远大于光束束腰否则光场会被边界截断采样点数N取2的幂1024、2048那种这样不管是用FFT还是角谱传播都方便。采样间隔dx必须满足奈奎斯特采样条件特别是超表面相位梯度很陡的时候如果dx过大高频相位变化会被直接抹平。我见过最典型的问题就是有人把L设得过大比如10mm但N只有256结果一个涡旋相位圆环上就十几个像素画出来的相位图跟马赛克一样中心暗核根本看不出来。经验值是束腰约1mm时仿真区域4mm配1024点是比较稳的组合最小特征尺寸大约4微米足够分辨绝大多数亚波长结构相位梯度的宏观效果。2.2 构建标准的拉盖尔-高斯涡旋光束很多教程会说“涡旋光就是exp(ilθ)”这只是相位因子。光束还得有强度包络否则直接用平面波乘exp(ilθ)也能得到涡旋相位但强度分布太均匀不贴近实际实验。最常用的是拉盖尔-高斯模式LG模式其复振幅化简到二维截面可以写成w0 1e-3; % 束腰半径 r sqrt(X.^2 Y.^2); theta atan2(Y, X); l 2; % 拓扑荷 E_lg (sqrt(2)*r/w0).^abs(l) ... .* exp(-r.^2/w0^2) ... .* exp(1i*l*theta);这个公式里的(sqrt(2)*r/w0).^abs(l)决定了高阶涡旋的强度环会随着拓扑荷增大而向外扩展而exp(-r.^2/w0^2)负责把整体能量约束在束腰范围内。实际用的时候还可以加上曲率因子和Gouy相位但对相位变化分析来说上面这个简化版本已经足够。为什么要强调用LG模式而不是纯exp(ilθ)因为纯相位项在r趋近零的地方不存在强度为零的约束数值上可能因为相位不连续导致暗核不够暗。LG模式在r0处自然为零符合真实涡旋光的物理特性。2.3 超表面相位屏的三种常见写法超表面的相位设计按功能可以分成几类在MATLAB里都需要写成二维相位矩阵。第一种是螺旋相位用于把高斯光转换成涡旋光或改变涡旋光拓扑荷l_meta 3; phase_spiral l_meta * atan2(Y, X);第二种是梯度相位让光束发生偏转对应超表面的“广义斯涅尔定律”常见于分束器设计period 500e-6; % 超表面单元周期对应的宏观相位周期 phase_gradient 2*pi/period * X;第三种是聚焦相位相当于超表面透镜用于汇聚光场f 10e-3; % 焦距 phase_lens -k * (X.^2 Y.^2) / (2*f);实际项目中超表面上的相位分布往往是几种基本相位的复合。比如既要偏转又要涡旋就写成θ_meta l_meta·θ 2π/Λ·x。这种线性叠加在物理上是合理的因为超表面每个单元只能给出一个确定的相位延迟这个相位延迟可以同时满足多个功能需求只要空间采样足够。要注意的是写成相位矩阵之后真正作用在光场上时要用exp(1i * phase)构造复振幅调制项。这一步别写成phase .* E那是错误的相位必须放在指数里。3. 把相位变化动态演出来传播与可视化3.1 相位叠加的核心操作入射涡旋光经过超表面后的出射光场在薄屏近似下就是两串复振幅逐点相乘E_out E_lg .* exp(1i * phase_meta);这行代码是整篇文章里最核心的东西。它的物理含义是超表面上每个位置对入射光的相位做一次局部修改透射光的等相位面因此被重新塑形。如果phase_meta本身是螺旋相位就会改变入射涡旋光的拓扑荷如果phase_meta是梯度相位就会给整个光束叠加一个横向波矢偏折如果phase_meta是透镜相位就会改变波前曲率。相位叠加之后最好在频域里看一眼角谱这样可以快速判断出射光是否混入了不需要的高频分量。尤其是梯度相位叠加涡旋相位时频谱中心会往一侧偏移偏折角度跟梯度周期成反比。3.2 用角谱法传播到任意距离只看零距离处的相位分布是不够的实际系统里探测器或者下一级光学元件都离超表面有一段距离。这里推荐用角谱法Angular Spectrum Method做光场传播因为它物理意义清晰MATLAB实现也方便。角谱法的思路是把空间光场分解成一系列平面波每个平面波传播一段距离后在频域乘上一个传播相位因子再逆变换回空间域。z 100e-3; % 传播距离 10cm k 2*pi / lambda; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); H exp(1i * k * z .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); E_far ifft2(ifftshift(fftshift(fft2(E_out)) .* H));这里面有几个容易出错的地方单独说一下。第一频域坐标的构造必须与空间域匹配。fftshift(fft2(E_out))之后零频在矩阵中心频域坐标也必须从负值到正值排列-N/2 : N/2-1而且频率间隔是1/L其中L是空间域仿真区域的边长。第二sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)里面如果出现负数表示倏逝波可以直接令其为零否则数值上会指数发散。虽然绝大多数传播距离下这个项不会碰到但当你把z设得极小或者网格设得很大时一定要加一个判断mask (1 - (lambda*FX).^2 - (lambda*FY).^2) 0; H zeros(N); H(mask) exp(1i * k * z .* sqrt(1 - (lambda*FX(mask)).^2 - (lambda*FY(mask)).^2));第三角谱法的适用距离很宽近场远场都行但注意它要求传播路径上是均匀介质不能穿过透镜或其他折射率突变界面。要模拟经过透镜后再传播就得拆成三段透镜前传播、透镜相位作用、透镜后传播。3.3 相位图、强度图、截面图组合出效果相位图是判断涡旋光质量最直观的方式。MATLAB里可以用imagesc显示相位figure; imagesc(x*1e3, x*1e3, angle(E_far)); axis xy; axis square; colormap(hsv); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(远场相位分布);这里用hsv色图而不是默认parula因为hsv在±π处是同一颜色首尾衔接特别适合显示周期性相位。如果用parula你会看到一条明显的颜色断层那不是物理问题是色图不匹配造成的视觉假象。强度图用imagesc(abs(E_far).^2)显示为了看清暗核通常会把强度做归一化再取对数或者直接叠加一个对比度增强I abs(E_far).^2; I I / max(I(:)); imagesc(x*1e3, x*1e3, log10(I 1e-4));相位变化不光要看二维分布还可以提取沿某个圆的相位值来验证拓扑荷。比如在以光轴为中心、半径r0处采样一圈相位unwrap之后如果线性变化2πl就说明拓扑荷确实是l。这个验证方法实验和仿真都通用。r0 0.5e-3; % 采样子圆半径 th_sample linspace(0, 2*pi, 360); x_sample r0*cos(th_sample); y_sample r0*sin(th_sample); phase_sample interp2(X, Y, angle(E_far), x_sample, y_sample, spline); phase_unwrap unwrap(phase_sample);如果你看到unwrap之后相位曲线不是单调上升或下降的斜线而是来回起伏说明这一圈相位不纯净可能是叠加了别的像差或者涡旋中心没对准网格原点。4. 实操过程一整套涡旋光超表面相位变化仿真4.1 参数初始化我习惯把所有可能需要调整的量集中放在脚本最前面方便扫描参数。这里给出一份可直接跑通的基础配置。clear; clc; close all; % 光场参数 lambda 632.8e-9; % 波长 (m) w0 1e-3; % 束腰半径 (m) l_inc 1; % 入射涡旋拓扑荷 % 超表面参数 l_meta 2; % 超表面叠加的拓扑荷 period_g 1.2e-3; % 梯度相位周期 (m)可选 f_meta 20e-3; % 超表面透镜焦距 (m)可选 % 仿真参数 L 4e-3; % 仿真区域边长 (m) N 1024; % 网格点数 z_prop 50e-3; % 传播距离 (m) dx L / N; x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); r sqrt(X.^2 Y.^2); theta atan2(Y, X); k 2*pi / lambda;把拓扑荷、焦距、周期都写在顶部的好处是后面想扫一个参数曲线时直接套一层for循环就行不用改核心代码。4.2 构建光场与超表面相位入射光场用LG模式超表面相位做一个“涡旋偏转”的复合设计这在很多光束操纵场景里是实用组合。% 入射涡旋光场 E_inc (sqrt(2)*r/w0).^abs(l_inc) .* exp(-r.^2/w0^2) .* exp(1i * l_inc * theta); % 超表面复合相位 phase_meta l_meta * theta 2*pi/period_g * X; % 出射光场 E_out E_inc .* exp(1i * phase_meta);复合相位的意思是这个超表面同时做两件事把入射的l1涡旋变成l3涡旋并让整个光束朝x方向偏转一个角度。偏转角理论上等于arcsin(λ/period_g)代入上面的参数大约是0.03度偏转量很小但在远场图上可以看出来。如果只关心纯涡旋变化把梯度项去掉就行。透镜聚焦项也一样需要就加不需要就删整个脚本是模块化的。4.3 传播与相位提取传播用角谱法按第三节的代码封装成一个函数会更方便。function E_prop asm_propagate(E_in, dx, lambda, z) N size(E_in, 1); L N * dx; fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx, fx); k 2*pi / lambda; H exp(1i * k * z .* sqrt(max(0, 1 - (lambda*FX).^2 - (lambda*FY).^2))); E_prop ifft2(ifftshift(fftshift(fft2(E_in)) .* H)); end调用方式很简单E_far asm_propagate(E_out, dx, lambda, z_prop);然后用前面提到的方法分别画入射相位、超表面相位、出射相位、传播后相位四张图排成2×2一眼就能看出相位变化链条。4.4 一组完整结果判读画出图之后建议按下面顺序检查结果。第一看入射光强度图中心应该有暗核暗核尺寸大约在波长量级周围是亮环。第二看出射光的相位分布如果超表面只加螺旋相位那么相位图绕中心一圈应该有l_incl_meta个2π跳变。如果加了梯度项相位图上会叠加一系列平行斜条纹这是偏转相位的表现。第三看传播后强度图涡旋中心依然保持暗核说明相位奇点稳定。如果暗核被填充大概率是角谱法里没做倏逝波滤波或者网格采样不够导致高频信息混叠。第四看提取的圆周相位unwrap之后斜率对应的拓扑荷数应该等于l_incl_meta。这是量化验证比肉眼看图靠谱得多。5. 常见问题与排查技巧实录5.1 相位图看起来一团乱全是细条纹出现这种情况通常是相位值的2π跳变太密导致imagesc的色图断层密集。在显示相位之前可以用wrapToPi对相位做缠绕处理或者直接用angle(exp(1i*phase))把相位压到[-π, π]区间。如果是显示传播后的相位由于光束整体有一个很强的平面波因子相位随空间快速振荡这种情况不适合直接看总相位应该先减去参考平面波相位比如k*z或者k*x*sinθ只看相对相位变化。5.2 中心暗核消失涡旋变成实心亮斑从仿真角度说原因八成是网格中心点没有精确落在坐标原点上。x (-N/2 : N/2-1) * dx在N为偶数时坐标点并不包含零导致涡旋中心落在四个像素的交叉点上相位奇点在离散网格中无法正确表达。解决办法是让N为奇数或者把坐标稍微平移半个像素x (-(N-1)/2 : (N-1)/2) * dx;这样中心像素就是0暗核不会被抹平。另一个可能原因是入射光场里加了过大的束腰w0导致暗核半径小于采样间隔中心看起来就是亮的。可以把w0调小或者把仿真区域L调大让暗核占据5个像素以上。5.3 角谱传播后能量明显减少角谱法本身是能量守恒的在无损耗条件下能量减少通常是因为边界截断。光场传播时会在空间上逐渐扩展如果仿真区域L不够大光场边缘会碰到边界FFT的周期性会把泄漏过来的场叠加到另一边造成数值误差。解决方法是把仿真区域设为光束尺寸的3到5倍或者像实际仿真里常用的做法一样在初始光场上乘一个超高斯窗口来抑制边缘衍射。还有一种情况是传播距离太长远场光斑直径超过L这种时候应该减小z_prop或者改用傅里叶变换法Fraunhofer衍射直接计算远场分布。5.4 拓扑荷提取不准unwrap后斜率不是整数这个问题的根源通常是采样本圆穿过涡旋中心导致相位跳变不连续。采样半径不能选太大也不能选太小太大会进入强度很低的区域相位噪声大太小会切割到暗核区域相位奇点附近相位本身就不连续。经验做法是把采样半径设为亮环的峰值位置大概在LG模式的半径处或者选束腰的0.7到1.0倍。另外interp2默认的线性插值在相位跳变边缘会出问题建议用spline插值并且先对复振幅插值再取相位角而不是直接对相位矩阵插值E_complex_sample interp2(X, Y, E_far, x_sample, y_sample, spline); phase_sample angle(E_complex_sample);这样处理之后unwrap的稳定性会好很多。5.5 复合相位扫描时结果不理想如果你在做参数扫描比如连续改变l_meta要注意绕回问题。相位从π到-π的跳变在数值上表现为一个很大的负向跳变unwrap可以解决但在做乘法叠加时可能会引入意想不到的伪影。建议所有相位操作都用复指数形式进行也就是把相位一直放在exp(1i*φ)里面等到最后需要显示相位时才取angle()。5.6 仿真配置问题速查表现象可能原因排查与解决相位图条纹过密相位周期性太强或显示方式不对使用wrapToPi压缩到[-π,π]或显示相对相位中心暗核消失N为偶数导致中心不在网格点改用奇数N或加半像素平移传播后能量减少边界截断或z过长增大L或改Fraunhofer远场法拓扑荷提取非整数采样半径不合适或插值方法不当在强度峰值半径处采样对复振幅插值再取相位低频背景太强入射高斯分量未完全转化为涡旋分量检查超表面相位幅值是否精确为2π整数倍周期图中有很强网格状条纹FFT采样混叠提高N或增大L检查频域坐标是否与空间域匹配6. 超表面相位设计的几个实用心得6.1 相位量化不一定非要8级有的同学一开始设计超表面单元就会陷入“8台阶量化”的思维定式因为传统的二元光学元件确实需要量化台阶。但这里要分清楚台阶量化是加工工艺限制不是物理必需。在MATLAB初步验证时完完全全可以先用连续相位做仿真确认系统级效果没问题再量化。如果量化后效果变差再考虑用一些优化算法来补偿量化误差。我实际做过一个对比连续相位、4台阶、8台阶、16台阶当单元的相位调控范围小于2π时比如只调控0.6π4台阶和8台阶的结果几乎没有差别。量化的主要损失发生在相位梯度很陡的区域那里相邻像素的相位跳变大量化台阶把本来就稀疏的采样信息又砍了一刀。6.2 孔径效应对相位变化的影响超表面不可能无限大但它有边界。这种有限孔径会给涡旋光引入衍射效应反映在远场就是强度环展宽、暗核边缘弥散。仿真的时候想模拟真实的有限孔径超表面只要在超表面相位上乘一个圆形孔径函数aperture_r 1.5e-3; aperture_mask double(r aperture_r); E_out E_inc .* exp(1i * phase_meta) .* aperture_mask;孔径边缘会导致高频衍射如果传播距离较长你会发现强度分布会在孔径边缘出现菲涅尔环纹这是正常现象。想要压低这个效应可以把硬边改成软边超高斯渐变但物理上实验加工出来的超表面基本都是硬边所以这个衍射效应其实也限制着实际系统的性能。6.3 别忘了偏振自由度涡旋光与超表面的结合很多时候还会涉及偏振。最典型的就是几何相位Pancharatnam-Berry相位超表面它的相位调控依赖单元的旋转角入射圆偏振光经过后转换为反向圆偏振光同时携带2倍旋转角的相位延迟。如果你研究的超表面是基于这类机制那么MATLAB仿真里就不能只用标量相位要考虑琼斯矩阵。做偏振相关仿真时推荐用琼斯矢量形式E_in_xy [E_x; E_y]; % 2xN矩阵两个偏振分量 J_meta [cos(alpha), -sin(alpha); sin(alpha), cos(alpha)]; % 单元旋转计算过程会比标量形式稍微复杂但思路一致每个像素位置用单元旋转角α生成一个局域琼斯矩阵作用在入射偏振态上然后分别传播两个正交分量。这样既能分析强度也能分析偏振转换效率。6.4 从远场反推超表面相位有一个前提很多人想做逆向设计给定想要的远场涡旋光反推超表面的相位分布。这个概念没错但有一个前提——超表面只能提供一个2π范围内的相位且这个相位对波前的作用是局部的。简单说它只能改“相位面”的形状不能凭空改变光场的振幅分布。所以逆推的时候如果目标光场的振幅分布和入射光差异很大单靠纯相位超表面做不好这时候要么加入振幅调控这会牺牲效率要么接受部分能量转换不理想。这个限制在做涡旋光生成时尤其明显高斯光转涡旋光效率天然受制于模式匹配。6.5 MATLAB仿真流程图在脑中的正确顺序把整套仿真浓缩成一句话光源构建、超表面相位叠加、角谱传播、相位/强度分析。这个顺序千万别乱。如果你先传播再叠加超表面相位等于把超表面放在传播路径末端物理场景完全不同。7. 从相位仿真到实际器件的衔接建议仿真不等于实际两者之间还有一大截路要走。我建议做完MATLAB系统级仿真之后再按这个顺序往下推进首先做单元级全波仿真CST或Lumerical FDTD验证你设计的超表面单元确实能实现所需相位的2π覆盖。单元结构的具体尺寸和材料参数在这个阶段确定仿真周期一般为亚波长微米量级。其次把单元仿真得到的相位-尺寸/旋转角关系表导出映射回MATLAB。这个映射很关键因为单元仿真出来的相位响应通常不是线性的直接按理想相位排布单元会导致实际相位误差产生明显的散射噪声。然后在MATLAB里做一个“真实相位”的二次验证用映射后的实际相位矩阵替换理想的phase_meta再次执行传播仿真。这一步能暴露很多单元设计问题比如相位覆盖不足、相邻单元耦合导致相位偏差等。最后才考虑加工和光学测量。到实验阶段你手上已经有一套“理想仿真→真实参数仿真”的对照数据测试结果出来也知道该怎么分析。这个流程很稳能避免那种“软件仿真看着完美加工出来完全不能用”的尴尬情况。我个人在实际操作中的体会是MATLAB最容易被低估的价值不是计算本身而是它让你在不到半小时内把整套物理逻辑过了一遍。相位构建、叠加、传播、分析每一步都在逼着你回答“这一步的物理意义是什么”等这些概念在脑子里彻底清晰了再上手任何专业商用软件都很快。分享一个小技巧收尾在你第一次跑通仿真后把入射拓扑荷改成0、把超表面改成纯螺旋相位对比一下出射光和入射光的相位分布你会在图上直观看到“一个平面波如何被超表面卷成涡旋”的整个过程。这个对比图特别适合用来检查自己的仿真链路是否正确也非常适合写论文时放给审稿人看——一张图解释清楚相位变化的核心机理。本文还有配套的精品资源点击获取