尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

用MATLAB绘制高斯光束三维光强分布与峰值定位

用MATLAB绘制高斯光束三维光强分布与峰值定位 简介面向光电、物理及电子信息类专业学生这份Matlab代码用于模拟高斯光束在空间中的三维光强分布与峰值分布可支撑课程设计、期末大作业或毕业设计中的光学仿真环节。资源包共4个文件其中2个为.m源程序分别对应主要仿真模型与辅助计算另外2张PNG图片展示了运行得到的三维光强与峰值分布效果图压缩包整体仅65KB轻量便于下载和复现。目前已有118人学习使用。代码采用参数化编程光束腰、波长、传播距离等关键参数均可方便调整注释明细、思路清晰即使Matlab基础一般的读者也能读懂修改在Matlab2014至2024a版本中均可直接运行省去环境配置的麻烦。通过该资源读者既能直观理解高斯光束的强度演化规律又可获得一套可复用的仿真脚本方便在此基础上拓展其他光束或光学系统模型。1. 为什么高斯光束的三维光强分布值得自己用 MATLAB 写一遍高斯光束的三维光强分布是激光物理课程设计里的高频题目但网上很多 MATLAB 代码只给一个二维截面或者把传播方向直接当成第三个坐标轴画出来的曲面像一段被拉亮的圆柱完全看不出瑞利距离内的发散过程。真正能交差的版本至少包含两个独立脚本gaosiguangshu00.m 画三维光强分布gaosiguangshu2.m 提取沿传播方向的峰值分布。这篇文章把这两个脚本的常见写法拆开从参数设定、网格生成到峰值定位逐段解释。适合正在做课程设计、期末大作业的同学也适合想用 MATLAB 快速验证高斯光束参数变化的人。2. 高斯光束参数化建模束腰、瑞利距离与强度公式模拟之前先要把强度表达式写对。很多示例代码把高斯分布写成I exp(-2*(x^2y^2)/w0^2)这个式子只描述某一个横截面上的光强放进三维图里所有传播位置的光斑半径都是同一个w0实际物理过程完全被抹掉了。基模高斯光束的强度应该同时依赖径向坐标r和传播坐标z。2.1 基模高斯光束的强度表达式忽略相位项后归一化强度可以写成I(r,z) I0 * (w0 / w(z))^2 * exp(-2 * r^2 / w(z)^2)其中w(z) w0 * sqrt(1 (z / zR)^2)zR pi * w0^2 / lambda。这个式子里的(w0/w(z))^2是峰值强度的衰减因子exp(-2r^2/w(z)^2)是每个截面上的高斯横向分布。zR叫瑞利距离是区分近场和远场的关键尺度z zR时光斑半径基本不变z zR后光斑近似按发散角线性扩大。做代码前先列一张参数表避免后面改参数时不知道动哪里。物理量符号代码变量案例取值单位波长λlambda1.064e-6m束腰半径w0w00.5e-3m瑞利距离zRzR由 lambda、w0 算出m初始峰值强度I0I01相对量传播方向采样zz-5zR 到 5zRm径向采样rr0 到 3w0m表里的波长选了 1064 nm这是 Nd:YAG 激光常见波长束腰取 0.5 mm对应普通实验光路。zR算出来约 0.738 m所以 z 轴取 ±3.7 m 才能覆盖从近场到远场的完整变化。2.2 参数化代码骨架把上面公式翻译成 MATLAB我习惯先建立两个一维坐标向量再用meshgrid扩展成矩阵最后用向量化运算一次算完。代码开头集中放参数后续要改只需要动前几行。% 文件 gaosiguangshu00.m lambda 1.064e-6; % 波长单位m w0 0.5e-3; % 束腰半径单位m I0 1; % 初始峰值强度归一化 zR pi * w0^2 / lambda; % 瑞利距离单位m z linspace(-5*zR, 5*zR, 400); % 传播方向采样点 r linspace(0, 3*w0, 200); % 径向采样点取到 3 倍束腰 [Z, R] meshgrid(z, r); % Z 和 R 同尺寸矩阵 w w0 * sqrt(1 (Z / zR).^2); % 每个 z 位置的光斑半径 I I0 * (w0 ./ w).^2 .* exp(-2 * (R ./ w).^2); % 三维强度矩阵代码逻辑说明meshgrid(z, r)生成两个矩阵Z的每一行是同一个 z 坐标R的每一列是同一个 r 坐标这样R./w才能逐元素计算径向归一化坐标。w0./w和(R./w).^2里的点运算符是矩阵逐位运算不能用/和^替代否则 MATLAB 会按线性代数规则处理结果完全错误。生成后可以用whos I查看矩阵尺寸预期是200x400对应length(r)行、length(z)列。这个维度顺序在后续用max提取峰值分布时很关键提前记住能省很多调试时间。2.3 常见参数误用与收敛性检查最容易犯的错是把w固定成w0这样画出来的三维强度曲面高度处处相同只是横向高斯包络被拉伸。正确做法里w必须随 Z 矩阵变化。另一个常被忽略的点是 z 轴向范围。如果zR很小比如束腰 0.05 mm、波长 1064 nm 时zR ≈ 7.4 mm坐标范围就要改成毫米量级。我的做法是先算zR再用zR作为尺度单位去设置linspace这样无论参数怎么改图像形状始终保持一致。若想输出发散角远场半角可以用theta lambda / (pi * w0)估算这个值在 1064 nm、0.5 mm 束腰下约 0.68 mrad可以在注释里写上方便和仿真图对照。3. 用 meshgrid 和 surf 把三维光强分布画出来强度矩阵算好之后绘图本身不复杂但坐标安排要清楚。这里的“三维”不是医学成像那种体渲染而是把传播方向z当作横轴、径向r当作纵轴、强度I当作高度得到一张三维曲面。这是 MATLAB 模拟高斯光束时最直观的表示方式课程设计和论文插图都够用。3.1 三维图的坐标安排如果沿用第二段的矩阵I的第一维是r第二维是z。surf(x, y, Z)要求Z是矩阵x和y可以是与 Z 同样尺寸的网格矩阵。我一般直接传z、r和I并让surf自动对齐网格。需要提醒的是网上有代码直接用surf(I)这样横纵坐标变成矩阵下标没有物理单位后面想标“传播距离 mm”就很麻烦。正确做法是显式传入坐标向量或网格矩阵。3.2 gaosiguangshu00.m 的绘图段沿用前面算好的I、z、r绘图段可以写成这样% 继续使用 gaosiguangshu00.m 中的变量 figure(Color, white); surf(z * 1e3, r * 1e3, I, EdgeColor, none); xlabel(传播距离 z (mm)); ylabel(径向 r (mm)); zlabel(归一化强度 I(r,z)); colormap(jet); colorbar; view([45, 30]);代码里把z和r都乘以1e3单位从米换算成毫米坐标轴标签更直观。EdgeColor,none去掉网格线曲面不会显得杂乱。view([45,30])把观察方向设为方位角 45 度、仰角 30 度能同时看见束腰收窄和峰值强度下降。jet是 MATLAB 经典配色如果学校报告要求黑白打印可以换colormap(parula)或直接去掉颜色条。3.3 叠加峰值轨迹和光斑半径包络只看曲面还不够把轴向峰值轨迹叠加上去能一眼看出强度衰减和光斑展宽的位置关系hold on; plot3(z * 1e3, zeros(size(z)), I0 * (w0 ./ w).^2, r--, LineWidth, 1.5); hold off; legend(三维光强分布, 轴向峰值强度);这段代码里的zeros(size(z))表示径向位置在 0 处也就是光轴第三个参数是每个 z 位置的峰值强度公式已经在第二章推过。红色虚线从三维曲面的“脊线”穿过去物理意义和后续峰值分布曲线一致只是这里把三维和轨迹放在同一张图里。绘图时还有一个细节网格点不要一开始就拉满。先用z取 150 点、r取 100 点确认角度和坐标范围满意后再把采样点提上去。surf对 2000×2000 矩阵的旋转交互会明显卡顿课程设计阶段没必要追求那种密度。4. 沿传播轴提取峰值分布并定位束腰峰值分布就是每个传播位置横截面上的最大光强随 z 的变化曲线。高斯光束任意横截面都是中心最强所以可以用解析式直接写出来I_peak(z) I0 * (w0 / w(z))^2。然而在真实项目里仿真数据往往来自数值计算或实验采集没有现成公式必须从矩阵里搜索最大值。gaosiguangshu2.m 做的就是这件事。4.1 从强度矩阵中搜索每个截面的最大值前面讲过I的尺寸是[length(r), length(z)]第一维是径向第二维是传播方向。用max沿第一维取最大就能得到一组随 z 变化的峰值强度。% 文件 gaosiguangshu2.m % 输入 I 来自 gaosiguangshu00.m peakI max(I, [], 1); % 沿径向求最大结果 1 x length(z) [peakMax, idx] max(peakI); % 全局最大强度及其索引 z_est z(idx); % 离散网格估计的束腰位置 figure(Color, white); plot(z * 1e3, peakI, b-, LineWidth, 1.2); xlabel(传播距离 z (mm)); ylabel(峰值强度); title(沿传播方向的峰值分布); grid on;max(I, [], 1)的第二个参数必须是 1表示沿第一维、也就是径向取最大值。如果写成max(I)MATLAB 会直接对全部元素求最大返回单个数值脚本就废了。peakMax理论上应该接近 1因为束腰处峰值强度就是I0如果明显偏离先回检查一下w0和lambda的单位是否统一。4.2 峰值位置的网格误差z(idx)直接作为束腰估计精度受linspace步长限制。z 范围取到 ±5zR采样 400 点时步长约 0.025 zR直接用max找到的束腰最多偏一个步长。对课程设计来说误差可接受但论文里如果需要精确束腰位置通常要做抛物线拟合。取峰值索引idx附近三个点用二次多项式拟合顶点位置就是更精确的束腰% 在峰值点附近做抛物线拟合修正离散误差 k idx; if k 1 k numel(z) p polyfit(z(k-1:k1)., peakI(k-1:k1)., 2); z_fit -p(2) / (2 * p(1)); % 抛物线顶点横坐标 endpolyfit输入必须是列向量所以z(k-1:k1)后面加了转置。抛物线顶点公式是-b/(2a)其中p(1)是二次项系数、p(2)是一次项系数。这里只取峰值附近三个点因为远离峰顶的曲线不再接近抛物线取太多点反而引入偏差。不同网格密度下的定位误差可以参考下面的数量级z 采样点数步长 / zR直接用 max 误差二次拟合误差1010.1约 0.05 zR约 0.002 zR4010.025约 0.013 zR约 0.0002 zR10010.01约 0.005 zR约 0.0001 zR表格说明一个结论网格越密直接max和拟合修正都越准但拟合方法在低采样率下提升更明显。如果峰值在采样区间边界拟合条件k 1 k numel(z)不成立这时需要扩大 z 范围重新计算。4.3 峰值分布与解析式对比gaosiguangshu2.m 除了画数据曲线我还会把解析式叠加上去用来检查数值计算是否正确I_analytic I0 * (w0 ./ w(1, :)).^2; % w 的第 1 行对应 r0 处实际上w所有行都相同因为光斑半径不依赖 r所以取w(1,:)即可。将I_analytic和peakI画在同一张图里如果两条线完全重合说明强度矩阵计算无误。实际数值结果里peakI应该等于中心位置r0的那一行强度两者是同一回事。5. 把课程设计代码改成自己版本的三步收尾技巧前面两个脚本已经能出图但交作业或写报告时还想做得更像样一点。最后说三个常用技巧直接改代码即可。5.1 用 subplot 把三维图和峰值分布拼在一起课程设计要“一张图看清所有结论”用subplot(1,2,1)和subplot(1,2,2)把两张图放同一窗口figure(Color, white); subplot(1,2,1); surf(z * 1e3, r * 1e3, I, EdgeColor, none); view([45, 30]); colorbar; xlabel(z (mm)); ylabel(r (mm)); subplot(1,2,2); plot(z * 1e3, peakI, b-, LineWidth, 1.2); xlabel(z (mm)); ylabel(峰值强度); grid on;两张图共用一个 z 坐标便于对比光斑发散和强度衰减。subplot中间不加结论文本直接交给读者看图比在标题里写解释性文字更专业。5.2 用 exportgraphics 导出高清位图MATLAB 2019a 及更早版本常用printR2020a 以后推荐这样导出exportgraphics(gcf, GaussianBeam_result.png, Resolution, 300);300 dpi 足够打印在课程设计报告里。如果要提交矢量图把后缀改成.pdf或.svg插入 Word 时不会出现锯齿。5.3 快速验证参数变化是否合理最后把束腰半径从0.5e-3改成0.2e-3重新运行两个脚本三维曲面会明显变得“瘦高”峰值分布曲线下降得更快束腰位置仍然保持在 z0 附近。这是因为束腰越小瑞利距离越短光束在更短的传播距离内发散掉。如果改成束腰变大峰值分布曲线则更平缓。这个验证能确认参数化编程没有写死任何中间量也适合在答辩时现场演示。本文还有配套的精品资源点击获取
返回列表