2026/8/13 12:59:14
FEATURED · 精选文章 
1. 项目概述从理论公式到可视化光束高斯光束这大概是光学和激光领域里最基础也最重要的概念之一了。但凡接触过激光原理、光纤通信或者光学设计都绕不开它。但说实话光看教科书上那一堆关于束腰、瑞利长度、发散角的公式总感觉隔着一层纱不够直观。我记得自己刚开始学的时候就特别想“看见”这个光束到底长什么样它在空间中是怎么传播和演变的。这就是为什么我觉得用MATLAB来做高斯光束的仿真特别有价值。它不是一个复杂的科研项目而是一个将抽象理论具象化的绝佳工具。通过几行代码你就能亲手“创造”出一束光观察它的强度分布看着它从束腰处最细的状态慢慢发散开来所有公式里的参数都变成了屏幕上可以调节的滑块和可以直观看到的图像变化。这个过程不仅能帮你牢牢记住那些关键概念更能培养一种“计算光学”的思维——用编程和数值计算来解决和验证光学问题。这篇内容我就想和你一起从最基础的高斯光束公式出发一步步在MATLAB里把它仿真出来。我们会涵盖一维、二维的强度分布三维的空间传播以及一些常用的分析技巧。无论你是正在学习相关课程的学生还是需要快速验证光学设计的工程师这个简单直接的仿真过程都能给你带来实实在在的帮助。2. 高斯光束理论基础与MATLAB实现思路在动手写代码之前我们得先统一一下“语言”明确我们要仿真的对象到底是什么。高斯光束顾名思义其横截面上的光强分布遵循高斯函数。它的数学描述是解析的这为我们的仿真提供了清晰的路径。2.1 核心公式拆解我们通常从最基本的傍轴近似下的标量波方程解说起。一个沿z轴传播的高斯光束其复振幅可以表示为U(x, y, z) A0 * (w0 / w(z)) * exp(-(x^2 y^2) / w(z)^2) * exp(-i*k*z - i*k*(x^2y^2)/(2*R(z)) i*ζ(z))这个公式看起来有点唬人但我们把它拆开看每个部分都有明确的物理意义振幅部分A0 * (w0 / w(z)) * exp(-(x^2 y^2) / w(z)^2)A0是常数代表振幅。w0是束腰半径即光束最细处的半径光强下降到中心值的1/e^2处的半径。w(z)是z位置处的光束半径它随传播距离变化w(z) w0 * sqrt(1 (z / zR)^2)。这里引入了瑞利长度zR π * w0^2 / λ它是衡量光束准直程度的重要参数可以理解为光束大致保持平行的传播距离。相位部分exp(-i*k*z - i*k*(x^2y^2)/(2*R(z)) i*ζ(z))exp(-i*k*z)是平面波的相位延迟。exp(-i*k*(x^2y^2)/(2*R(z)))是波前曲率带来的相位因子R(z)是波前曲率半径R(z) z * (1 (zR / z)^2)。在束腰处(z0)R(z)为无穷大波前是平面远离束腰波前变为球面。exp(i*ζ(z))是古依相位ζ(z) atan(z / zR)这是一个额外的相位移动在光束通过焦点时会发生π的相位跃变。我们最常关心的是光强分布I(x, y, z) |U(x, y, z)|^2。对于许多仿真我们可以暂时忽略相位部分专注于强度仿真公式简化为I(x, y, z) I0 * (w0 / w(z))^2 * exp(-2*(x^2 y^2) / w(z)^2)其中I0是束腰中心处的峰值光强。2.2 MATLAB仿真策略设计有了公式如何在MATLAB中构建它我的思路是分层次、模块化地进行参数化定义将所有物理参数λ,w0,zR等设为变量。这样我们只需要改变几个参数就能立刻看到不同波长、不同束腰大小的光束行为仿真就活了。网格化空间使用meshgrid函数创建代表x, y, z坐标的数值网格。这是数值计算的基础我们的光束将在这个网格上被“计算”出来。向量化计算利用MATLAB强大的矩阵运算能力避免使用低效的循环。直接对整个坐标网格应用上述公式一次性计算出整个空间的光场或光强分布。这是MATLAB仿真的性能关键。可视化呈现这是点睛之笔。我们将使用imagesc或pcolor来显示二维截面图。surf或mesh来绘制三维强度分布。plot来绘制一维光强剖面用于精确测量束宽。subplot将多角度视图组合在一起形成全面的分析面板。实操心得在开始写代码前最好在纸上或注释里把公式和计算步骤列清楚。尤其是w(z)和R(z)这些z的函数确保你在计算每个z平面时都用对了对应的值。一个常见的错误是在计算三维体数据时错误地广播了数组维度导致结果不对。3. 从二维截面到三维传播完整仿真实现理论清晰了策略定好了现在打开MATLAB我们开始动手实现。我会从最简单的静态二维截面开始逐步扩展到完整的三维空间传播仿真。3.1 基础二维束腰光强分布仿真我们从仿真的起点——束腰位置z0的二维光强分布开始。这是最直观的一步。%% 1. 基础参数设置 clear; close all; clc; % 清空环境好习惯 lambda 632.8e-9; % 波长单位米这里用常见的He-Ne激光波长 w0 1e-3; % 束腰半径1毫米 I0 1; % 中心峰值光强归一化为1 % 设置观察平面范围和采样点数 x linspace(-3*w0, 3*w0, 500); % x轴范围-3倍束腰到3倍束腰 y linspace(-3*w0, 3*w0, 500); % y轴范围 [X, Y] meshgrid(x, y); % 生成二维网格坐标 %% 2. 计算束腰处光强分布 (z0, 此时 w(z)w0) % 二维高斯分布公式 I_xy I0 * exp(-2*(X.^2 Y.^2) / w0^2); %% 3. 可视化 figure(Position, [100, 100, 1200, 400]); % 设置一个宽幅图窗 % 子图1二维伪彩色图 subplot(1, 3, 1); imagesc(x*1e3, y*1e3, I_xy); % 坐标转换为毫米显示 axis image; % 保持纵横比相等 xlabel(x (mm)); ylabel(y (mm)); title(束腰处光强分布 (二维视图)); colorbar; colormap(hot); % 使用‘hot’色图模拟光强 % 子图2三维曲面图 subplot(1, 3, 2); surf(X*1e3, Y*1e3, I_xy, EdgeColor, none); xlabel(x (mm)); ylabel(y (mm)); zlabel(相对光强); title(束腰处光强分布 (三维视图)); colormap(hot); view(30, 30); % 调整视角 % 子图3通过中心的一维剖面 subplot(1, 3, 3); plot(x*1e3, I_xy(ceil(end/2), :), LineWidth, 2); % 取中间一行数据 xlabel(x (mm)); ylabel(相对光强); title(束腰处水平中心光强剖面); grid on; hold on; % 标记1/e^2强度点对应束腰半径w0 plot([-w0, w0]*1e3, [exp(-2), exp(-2)], r--, LineWidth, 1.5); legend(光强分布, 1/e^2 阈值, Location, best);这段代码运行后你会得到一个包含三个视图的图形。二维图让你看到经典的高斯光斑三维图让你感受强度的起伏而一维剖面图则是定量分析的关键。你可以从剖面图上精确读出光强下降到中心值1/e^2约0.135的位置那就是±w0。注意事项linspace生成的坐标点数量决定了图像的分辨率和计算量。500个点对于演示是清晰的但如果你的计算区域很大或者需要非常精细可以增加到1000或更多。同时注意meshgrid生成的X和Y是矩阵这确保了后续的点乘(.^2)和矩阵运算能正确进行。3.2 三维空间传播仿真静态截面看完了我们让光束“动”起来看看它怎么传播。这意味着我们要计算一系列不同z位置上的二维光强分布然后堆叠或选取视图来观察。%% 1. 扩展参数定义传播轴 z linspace(-3*zR, 3*zR, 150); % 传播距离范围-3倍瑞利长度到3倍瑞利长度 % 计算瑞利长度 zR pi * w0^2 / lambda; fprintf(束腰半径 w0 %.2f mm, 瑞利长度 zR %.2f mm\n, w0*1e3, zR*1e3); % 创建三维空间网格注意全三维网格数据量巨大这里采用循环计算每个z面 % 我们依然使用之前的x-y网格 I_3D zeros(length(y), length(x), length(z)); % 预分配内存提升效率 %% 2. 循环计算每个z平面的光强分布 for k 1:length(z) z_current z(k); % 计算当前z处的光束半径w(z)和波前曲率半径R(z)此处仅用于强度R(z)未用 w_z w0 * sqrt(1 (z_current / zR)^2); % 计算当前z平面的光强分布 I_3D(:, :, k) I0 * (w0 / w_z)^2 .* exp(-2*(X.^2 Y.^2) / w_z^2); end %% 3. 三维传播可视化 figure(Position, [100, 100, 1400, 500]); % 子图1X-Z平面y0的切片图 subplot(1, 3, 1); % 提取y0中心线上的数据 [~, idx_y0] min(abs(y - 0)); % 找到y坐标最接近0的索引 I_xz squeeze(I_3D(idx_y0, :, :)); % 挤压维度得到(x, z)矩阵 imagesc(z*1e3, x*1e3, I_xz); xlabel(传播距离 z (mm)); ylabel(x (mm)); title(X-Z平面 (y0) 光强分布); axis xy; % 确保y轴方向正确 colorbar; colormap(hot); % 标记束腰位置(z0)和瑞利长度位置 hold on; plot([0,0], [x(1), x(end)]*1e3, w--, LineWidth, 1); plot([-zR, zR]*1e3, [0,0], c--, LineWidth, 1); legend(, 束腰位置 z0, 瑞利长度 ±zR, Location, best); % 子图2不同z位置的光束半径变化曲线 subplot(1, 3, 2); w_z_array w0 * sqrt(1 (z / zR).^2); % 理论值 plot(z*1e3, w_z_array*1e3, b-, LineWidth, 2); xlabel(传播距离 z (mm)); ylabel(光束半径 w(z) (mm)); title(光束半径随传播距离变化); grid on; hold on; % 标记束腰和渐近线 plot([z(1), z(end)]*1e3, [w0, w0]*1e3, r--, LineWidth, 1); plot([z(1), z(end)]*1e3, [w0/zR * abs(z(end)), w0/zR * abs(z(end))]*1e3, g--, LineWidth, 1); legend(w(z) 理论曲线, 束腰半径 w0, 渐近线 (发散角), Location, best); % 子图3选取三个特征z位置显示其二维光斑 subplot(1, 3, 3); z_plots [-zR, 0, zR]; % 选择 z -zR, 0, zR 三个面 plot_titles {z -zR (束腰上游), z 0 (束腰处), z zR (束腰下游)}; for p 1:3 [~, idx_z] min(abs(z - z_plots(p))); % 找到最接近的索引 I_slice I_3D(:, :, idx_z); % 创建一个临时坐标为了在子图中显示 ax subplot(3, 3, 6p); % 定位到右列的下三个子图位置 imagesc(ax, x*1e3, y*1e3, I_slice); axis(ax, image); xlabel(ax, x (mm)); ylabel(ax, y (mm)); title(ax, plot_titles{p}); colormap(ax, hot); if p 3 colorbar(eastoutside); end end这段代码生成了更丰富的视图。X-Z平面切片图清晰地展示了光束从汇聚到束腰再到发散的全过程你能看到光束宽度如何平滑变化。光束半径变化曲线则定量验证了w(z)的双曲线公式。最后对比三个特征位置的光斑直观展示了“束腰处最细离束腰越远光斑越大”的规律。实操心得计算三维数据I_3D时预分配数组zeros至关重要。如果是在循环内动态扩展数组当z的点数很多时MATLAB会反复申请内存导致速度极慢甚至内存不足。另外对于非常大的三维网格直接存储所有数据可能内存吃不消。这时可以采用“按需计算”的策略即只在需要画图的时候计算当前切片的数据而不是一次性计算全部。4. 进阶分析与常见问题调试基础仿真跑通后我们可以玩点更深入的同时也会遇到一些典型问题。这部分是区分“会写代码”和“理解仿真”的关键。4.1 相位信息仿真与波前可视化之前我们只关注了强度。要完整描述光场相位必不可少尤其是在涉及干涉、衍射或光束整形时。%% 1. 计算包含相位的复振幅场 (在束腰处z0) z_target 0; % 选择束腰位置 k 2 * pi / lambda; % 波数 % 注意在束腰处R(z) - inf因此波前曲率相位项为0。古依相位 ζ(z)0。 U_xy sqrt(I0) * exp(-(X.^2 Y.^2) / w0^2); % 复振幅此时相位均匀 %% 2. 计算离开束腰后的复振幅场 (例如 z zR) z_target zR; w_z w0 * sqrt(1 (z_target / zR)^2); R_z z_target * (1 (zR / z_target)^2); % 波前曲率半径 zeta_z atan(z_target / zR); % 古依相位 % 完整的复振幅计算 U_xy_z sqrt(I0) * (w0 / w_z) .* exp(-(X.^2 Y.^2) / w_z^2) .* ... exp(-1i * k * z_target) .* ... exp(-1i * k * (X.^2 Y.^2) / (2 * R_z)) .* ... exp(1i * zeta_z); %% 3. 可视化相位 phase_z angle(U_xy_z); % 提取相位范围 [-π, π] figure; subplot(1,2,1); imagesc(x*1e3, y*1e3, phase_z); axis image; colorbar; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(z %.1f mm 处的相位分布, z_target*1e3)); colormap(hsv); % HSV色图非常适合表示周期性相位 subplot(1,2,2); % 绘制通过中心的相位剖面 plot(x*1e3, phase_z(ceil(end/2), :), LineWidth, 2); xlabel(x (mm)); ylabel(相位 (弧度)); title(相位中心剖面); grid on; % 可以看到相位随x^2变化的抛物线形状这正是球面波前的特征。4.2 常见仿真问题与排查技巧在实际操作中你可能会遇到下面这些问题。这里是我的排查清单问题现象可能原因排查与解决方法图形显示一片空白或全黑/全白1. 数据值超出显示范围或全部为0/NaN。2.imagesc的坐标轴设置错误数据画到看不见的地方。3. 色图范围不合适。1. 在imagesc后使用caxis([0, 1])手动设置颜色轴范围。或用max(I_xy(:))检查数据最大值。2. 检查imagesc(x, y, I)中的x, y向量是否与矩阵I的维度匹配size(I,2)length(x),size(I,1)length(y)。3. 尝试colormap(gray)或imagesc(..., CDataMapping, scaled)。三维曲面图非常粗糙或畸变网格点太少 (linspace点数不足)。增加linspace的采样点数例如从100增加到500。对于曲面图也可以使用shading interp命令进行平滑插值。计算速度极慢1. 使用了多重嵌套for循环。2. 未预分配大数组。3. 网格分辨率过高。1.务必使用向量化操作。将exp,sqrt等函数直接作用于整个矩阵X,Y。2. 对于像I_3D这样的结果数组使用zeros预分配。3. 在保证可视效果的前提下适当降低采样点数。光束形状不对称或奇怪1. x和y的坐标范围或采样点数不一致。2. 公式中输入了错误的坐标变量。1. 确保linspace(-range, range, N)中的range和N在x和y方向一致除非你特意需要各向异性的光束。2. 检查公式中是X.^2 Y.^2而不是x.^2 y.^2。前者是矩阵后者是向量。“发散角”与理论值对不上1. 在图中测量发散角的方法不对。2.w0或λ的单位不一致。1. 理论发散角远场θ λ / (π * w0)。在仿真中可通过拟合远场w(z) ≈ θ *相位图出现剧烈跳变锯齿相位包裹问题。真实相位变化超过[-π, π]范围angle函数将其折叠到此区间。这是正常现象表示相位变化剧烈。若要分析连续的相位变化可以使用unwrap函数对一维数据效果好但对于二维相位解包裹需要更复杂的算法这属于另一个专题。可视化时用hsv色图即可清晰观察相位周期。避坑技巧在编写每个关键公式后用一两行代码做个快速验证。比如计算完w_z后fprintf(在z%.2f处w_z%.2e\n, z_target, w_z)看看输出是否合理。对于光束半径在束腰处应该等于w0在zzR处应该等于sqrt(2)*w0。这种即时的小验证能帮你快速定位是公式写错了还是参数单位错了。5. 仿真扩展与应用实例一个扎实的基础仿真框架建立后你可以像搭积木一样扩展出很多有趣且实用的应用。5.1 模拟光束通过透镜的变换这是非常经典的应用。一个薄透镜可以改变高斯光束的波前曲率从而变换其束腰位置和大小。%% 模拟透镜对高斯光束的变换 f 0.1; % 透镜焦距100mm z_from_lens 0.2; % 透镜到输入光束束腰的距离200mm % 假设入射光束参数在透镜位置处 lambda 632.8e-9; w0_in 0.5e-3; % 入射光束束腰半径 zR_in pi * w0_in^2 / lambda; % 计算入射光束在透镜位置的光斑大小w_in和波前曲率半径R_in w_in w0_in * sqrt(1 (z_from_lens / zR_in)^2); R_in z_from_lens * (1 (zR_in / z_from_lens)^2); % 应用ABCD定律薄透镜变换 % 对于薄透镜传输矩阵为 [1, 0; -1/f, 1] % 复数q参数满足 1/q2 1/q1 - 1/f % 其中 q z i*zR 或者用更常用的形式1/q 1/R - i*λ/(π*w^2) % 计算入射光束的q参数 q1_inv 1/R_in - 1i*lambda/(pi*w_in^2); % 经过透镜后的q参数 q2_inv q1_inv - 1/f; % 从q2参数提取出射光束的R_out和w_out q2 1 / q2_inv; R_out 1 / real(1/q2); w_out sqrt(-lambda / (pi * imag(1/q2))); % 进一步可以计算出射光束的新束腰位置相对于透镜和束腰大小 z_out R_out / (1 (lambda*R_out/(pi*w_out^2))^2); w0_out w_out / sqrt(1 (pi*w_out^2/(lambda*R_out))^2); fprintf(经过焦距%.3f m的透镜后\n, f); fprintf( 出射光束在透镜处的光斑半径: %.3f mm\n, w_out*1e3); fprintf( 出射光束在透镜处的波前曲率半径: %.3f m\n, R_out); fprintf( 新的束腰位置相对于透镜: %.3f m\n, z_out); fprintf( 新的束腰半径: %.3f mm\n, w0_out*1e3);这个例子展示了如何将高斯光束的q参数与ABCD矩阵理论结合用代码实现复杂光学系统的仿真。你可以改变焦距f和入射距离z_from_lens观察输出光束参数如何变化这对于光学系统设计中的光路计算非常有帮助。5.2 生成动态传播GIF动画静态图片虽然好但动态图更能直观展示传播过程。我们可以用循环生成一系列图片然后合成GIF。%% 生成光束传播动态GIF filename gaussian_beam_propagation.gif; % 沿用之前的三维数据 I_3D 或者重新计算一段更短的z轴用于动画 z_anim linspace(-2*zR, 2*zR, 50); figure(Position, [100, 100, 800, 600], Color, white); for n 1:length(z_anim) % 计算或提取第n帧的数据 (这里假设重新计算) w_z w0 * sqrt(1 (z_anim(n) / zR)^2); I_frame I0 * (w0 / w_z)^2 .* exp(-2*(X.^2 Y.^2) / w_z^2); % 绘图 imagesc(x*1e3, y*1e3, I_frame); axis image; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(高斯光束传播 | z %.2f mm, z_anim(n)*1e3)); colormap(hot); colorbar; caxis([0, 1]); % 固定颜色轴让动画稳定 drawnow; % 捕获帧并写入GIF frame getframe(gcf); im frame2im(frame); [imind, cm] rgb2ind(im, 256); if n 1 imwrite(imind, cm, filename, gif, Loopcount, inf, DelayTime, 0.1); else imwrite(imind, cm, filename, gif, WriteMode, append, DelayTime, 0.1); end end fprintf(GIF动画已保存为: %s\n, filename);这个动画能非常生动地展示光束从汇聚到发散束腰处最细的整个过程。生成的GIF可以用于报告、演示或教学材料中。5.3 像散与椭圆高斯光束仿真现实中的激光器可能输出非圆对称的即椭圆高斯光束。它的两个正交方向x和y有不同的束腰和发散角。%% 仿真椭圆高斯光束 w0x 1.0e-3; % x方向束腰半径 w0y 0.5e-3; % y方向束腰半径 zRx pi * w0x^2 / lambda; zRy pi * w0y^2 / lambda; % 计算特定z位置例如z0的光强 z_astig 0; wx_z w0x * sqrt(1 (z_astig / zRx)^2); wy_z w0y * sqrt(1 (z_astig / zRy)^2); % 光强分布为两个方向高斯函数的乘积 I_ellipse I0 * exp(-2*(X.^2)/wx_z^2 - 2*(Y.^2)/wy_z^2); figure; imagesc(x*1e3, y*1e3, I_ellipse); axis image; xlabel(x (mm)); ylabel(y (mm)); title(椭圆高斯光束 (束腰处)); colorbar; colormap(hot); % 可以看到光斑在x和y方向宽度不同呈椭圆形。通过调整w0x和w0y的比例你可以仿真出从细线到扁椭圆的各种光束形状。更进一步你还可以让两个方向的束腰位置不在同一z点这就是像散光束的典型特征仿真方法类似只是分别计算x和y方向的光束半径w(z)时使用不同的z偏移量。这些扩展实例只是抛砖引玉。基于这个核心框架你还可以去仿真高阶模如厄米-高斯或拉盖尔-高斯光束模拟光束在湍流大气中的传播或者与傅里叶光学结合计算光束的衍射。关键是理解物理模型然后用MATLAB这个强大的计算工具把它实现出来。
MATLAB仿真高斯光束：从理论公式到可视化传播
1. 项目概述从理论公式到可视化光束高斯光束这大概是光学和激光领域里最基础也最重要的概念之一了。但凡接触过激光原理、光纤通信或者光学设计都绕不开它。但说实话光看教科书上那一堆关于束腰、瑞利长度、发散角的公式总感觉隔着一层纱不够直观。我记得自己刚开始学的时候就特别想“看见”这个光束到底长什么样它在空间中是怎么传播和演变的。这就是为什么我觉得用MATLAB来做高斯光束的仿真特别有价值。它不是一个复杂的科研项目而是一个将抽象理论具象化的绝佳工具。通过几行代码你就能亲手“创造”出一束光观察它的强度分布看着它从束腰处最细的状态慢慢发散开来所有公式里的参数都变成了屏幕上可以调节的滑块和可以直观看到的图像变化。这个过程不仅能帮你牢牢记住那些关键概念更能培养一种“计算光学”的思维——用编程和数值计算来解决和验证光学问题。这篇内容我就想和你一起从最基础的高斯光束公式出发一步步在MATLAB里把它仿真出来。我们会涵盖一维、二维的强度分布三维的空间传播以及一些常用的分析技巧。无论你是正在学习相关课程的学生还是需要快速验证光学设计的工程师这个简单直接的仿真过程都能给你带来实实在在的帮助。2. 高斯光束理论基础与MATLAB实现思路在动手写代码之前我们得先统一一下“语言”明确我们要仿真的对象到底是什么。高斯光束顾名思义其横截面上的光强分布遵循高斯函数。它的数学描述是解析的这为我们的仿真提供了清晰的路径。2.1 核心公式拆解我们通常从最基本的傍轴近似下的标量波方程解说起。一个沿z轴传播的高斯光束其复振幅可以表示为U(x, y, z) A0 * (w0 / w(z)) * exp(-(x^2 y^2) / w(z)^2) * exp(-i*k*z - i*k*(x^2y^2)/(2*R(z)) i*ζ(z))这个公式看起来有点唬人但我们把它拆开看每个部分都有明确的物理意义振幅部分A0 * (w0 / w(z)) * exp(-(x^2 y^2) / w(z)^2)A0是常数代表振幅。w0是束腰半径即光束最细处的半径光强下降到中心值的1/e^2处的半径。w(z)是z位置处的光束半径它随传播距离变化w(z) w0 * sqrt(1 (z / zR)^2)。这里引入了瑞利长度zR π * w0^2 / λ它是衡量光束准直程度的重要参数可以理解为光束大致保持平行的传播距离。相位部分exp(-i*k*z - i*k*(x^2y^2)/(2*R(z)) i*ζ(z))exp(-i*k*z)是平面波的相位延迟。exp(-i*k*(x^2y^2)/(2*R(z)))是波前曲率带来的相位因子R(z)是波前曲率半径R(z) z * (1 (zR / z)^2)。在束腰处(z0)R(z)为无穷大波前是平面远离束腰波前变为球面。exp(i*ζ(z))是古依相位ζ(z) atan(z / zR)这是一个额外的相位移动在光束通过焦点时会发生π的相位跃变。我们最常关心的是光强分布I(x, y, z) |U(x, y, z)|^2。对于许多仿真我们可以暂时忽略相位部分专注于强度仿真公式简化为I(x, y, z) I0 * (w0 / w(z))^2 * exp(-2*(x^2 y^2) / w(z)^2)其中I0是束腰中心处的峰值光强。2.2 MATLAB仿真策略设计有了公式如何在MATLAB中构建它我的思路是分层次、模块化地进行参数化定义将所有物理参数λ,w0,zR等设为变量。这样我们只需要改变几个参数就能立刻看到不同波长、不同束腰大小的光束行为仿真就活了。网格化空间使用meshgrid函数创建代表x, y, z坐标的数值网格。这是数值计算的基础我们的光束将在这个网格上被“计算”出来。向量化计算利用MATLAB强大的矩阵运算能力避免使用低效的循环。直接对整个坐标网格应用上述公式一次性计算出整个空间的光场或光强分布。这是MATLAB仿真的性能关键。可视化呈现这是点睛之笔。我们将使用imagesc或pcolor来显示二维截面图。surf或mesh来绘制三维强度分布。plot来绘制一维光强剖面用于精确测量束宽。subplot将多角度视图组合在一起形成全面的分析面板。实操心得在开始写代码前最好在纸上或注释里把公式和计算步骤列清楚。尤其是w(z)和R(z)这些z的函数确保你在计算每个z平面时都用对了对应的值。一个常见的错误是在计算三维体数据时错误地广播了数组维度导致结果不对。3. 从二维截面到三维传播完整仿真实现理论清晰了策略定好了现在打开MATLAB我们开始动手实现。我会从最简单的静态二维截面开始逐步扩展到完整的三维空间传播仿真。3.1 基础二维束腰光强分布仿真我们从仿真的起点——束腰位置z0的二维光强分布开始。这是最直观的一步。%% 1. 基础参数设置 clear; close all; clc; % 清空环境好习惯 lambda 632.8e-9; % 波长单位米这里用常见的He-Ne激光波长 w0 1e-3; % 束腰半径1毫米 I0 1; % 中心峰值光强归一化为1 % 设置观察平面范围和采样点数 x linspace(-3*w0, 3*w0, 500); % x轴范围-3倍束腰到3倍束腰 y linspace(-3*w0, 3*w0, 500); % y轴范围 [X, Y] meshgrid(x, y); % 生成二维网格坐标 %% 2. 计算束腰处光强分布 (z0, 此时 w(z)w0) % 二维高斯分布公式 I_xy I0 * exp(-2*(X.^2 Y.^2) / w0^2); %% 3. 可视化 figure(Position, [100, 100, 1200, 400]); % 设置一个宽幅图窗 % 子图1二维伪彩色图 subplot(1, 3, 1); imagesc(x*1e3, y*1e3, I_xy); % 坐标转换为毫米显示 axis image; % 保持纵横比相等 xlabel(x (mm)); ylabel(y (mm)); title(束腰处光强分布 (二维视图)); colorbar; colormap(hot); % 使用‘hot’色图模拟光强 % 子图2三维曲面图 subplot(1, 3, 2); surf(X*1e3, Y*1e3, I_xy, EdgeColor, none); xlabel(x (mm)); ylabel(y (mm)); zlabel(相对光强); title(束腰处光强分布 (三维视图)); colormap(hot); view(30, 30); % 调整视角 % 子图3通过中心的一维剖面 subplot(1, 3, 3); plot(x*1e3, I_xy(ceil(end/2), :), LineWidth, 2); % 取中间一行数据 xlabel(x (mm)); ylabel(相对光强); title(束腰处水平中心光强剖面); grid on; hold on; % 标记1/e^2强度点对应束腰半径w0 plot([-w0, w0]*1e3, [exp(-2), exp(-2)], r--, LineWidth, 1.5); legend(光强分布, 1/e^2 阈值, Location, best);这段代码运行后你会得到一个包含三个视图的图形。二维图让你看到经典的高斯光斑三维图让你感受强度的起伏而一维剖面图则是定量分析的关键。你可以从剖面图上精确读出光强下降到中心值1/e^2约0.135的位置那就是±w0。注意事项linspace生成的坐标点数量决定了图像的分辨率和计算量。500个点对于演示是清晰的但如果你的计算区域很大或者需要非常精细可以增加到1000或更多。同时注意meshgrid生成的X和Y是矩阵这确保了后续的点乘(.^2)和矩阵运算能正确进行。3.2 三维空间传播仿真静态截面看完了我们让光束“动”起来看看它怎么传播。这意味着我们要计算一系列不同z位置上的二维光强分布然后堆叠或选取视图来观察。%% 1. 扩展参数定义传播轴 z linspace(-3*zR, 3*zR, 150); % 传播距离范围-3倍瑞利长度到3倍瑞利长度 % 计算瑞利长度 zR pi * w0^2 / lambda; fprintf(束腰半径 w0 %.2f mm, 瑞利长度 zR %.2f mm\n, w0*1e3, zR*1e3); % 创建三维空间网格注意全三维网格数据量巨大这里采用循环计算每个z面 % 我们依然使用之前的x-y网格 I_3D zeros(length(y), length(x), length(z)); % 预分配内存提升效率 %% 2. 循环计算每个z平面的光强分布 for k 1:length(z) z_current z(k); % 计算当前z处的光束半径w(z)和波前曲率半径R(z)此处仅用于强度R(z)未用 w_z w0 * sqrt(1 (z_current / zR)^2); % 计算当前z平面的光强分布 I_3D(:, :, k) I0 * (w0 / w_z)^2 .* exp(-2*(X.^2 Y.^2) / w_z^2); end %% 3. 三维传播可视化 figure(Position, [100, 100, 1400, 500]); % 子图1X-Z平面y0的切片图 subplot(1, 3, 1); % 提取y0中心线上的数据 [~, idx_y0] min(abs(y - 0)); % 找到y坐标最接近0的索引 I_xz squeeze(I_3D(idx_y0, :, :)); % 挤压维度得到(x, z)矩阵 imagesc(z*1e3, x*1e3, I_xz); xlabel(传播距离 z (mm)); ylabel(x (mm)); title(X-Z平面 (y0) 光强分布); axis xy; % 确保y轴方向正确 colorbar; colormap(hot); % 标记束腰位置(z0)和瑞利长度位置 hold on; plot([0,0], [x(1), x(end)]*1e3, w--, LineWidth, 1); plot([-zR, zR]*1e3, [0,0], c--, LineWidth, 1); legend(, 束腰位置 z0, 瑞利长度 ±zR, Location, best); % 子图2不同z位置的光束半径变化曲线 subplot(1, 3, 2); w_z_array w0 * sqrt(1 (z / zR).^2); % 理论值 plot(z*1e3, w_z_array*1e3, b-, LineWidth, 2); xlabel(传播距离 z (mm)); ylabel(光束半径 w(z) (mm)); title(光束半径随传播距离变化); grid on; hold on; % 标记束腰和渐近线 plot([z(1), z(end)]*1e3, [w0, w0]*1e3, r--, LineWidth, 1); plot([z(1), z(end)]*1e3, [w0/zR * abs(z(end)), w0/zR * abs(z(end))]*1e3, g--, LineWidth, 1); legend(w(z) 理论曲线, 束腰半径 w0, 渐近线 (发散角), Location, best); % 子图3选取三个特征z位置显示其二维光斑 subplot(1, 3, 3); z_plots [-zR, 0, zR]; % 选择 z -zR, 0, zR 三个面 plot_titles {z -zR (束腰上游), z 0 (束腰处), z zR (束腰下游)}; for p 1:3 [~, idx_z] min(abs(z - z_plots(p))); % 找到最接近的索引 I_slice I_3D(:, :, idx_z); % 创建一个临时坐标为了在子图中显示 ax subplot(3, 3, 6p); % 定位到右列的下三个子图位置 imagesc(ax, x*1e3, y*1e3, I_slice); axis(ax, image); xlabel(ax, x (mm)); ylabel(ax, y (mm)); title(ax, plot_titles{p}); colormap(ax, hot); if p 3 colorbar(eastoutside); end end这段代码生成了更丰富的视图。X-Z平面切片图清晰地展示了光束从汇聚到束腰再到发散的全过程你能看到光束宽度如何平滑变化。光束半径变化曲线则定量验证了w(z)的双曲线公式。最后对比三个特征位置的光斑直观展示了“束腰处最细离束腰越远光斑越大”的规律。实操心得计算三维数据I_3D时预分配数组zeros至关重要。如果是在循环内动态扩展数组当z的点数很多时MATLAB会反复申请内存导致速度极慢甚至内存不足。另外对于非常大的三维网格直接存储所有数据可能内存吃不消。这时可以采用“按需计算”的策略即只在需要画图的时候计算当前切片的数据而不是一次性计算全部。4. 进阶分析与常见问题调试基础仿真跑通后我们可以玩点更深入的同时也会遇到一些典型问题。这部分是区分“会写代码”和“理解仿真”的关键。4.1 相位信息仿真与波前可视化之前我们只关注了强度。要完整描述光场相位必不可少尤其是在涉及干涉、衍射或光束整形时。%% 1. 计算包含相位的复振幅场 (在束腰处z0) z_target 0; % 选择束腰位置 k 2 * pi / lambda; % 波数 % 注意在束腰处R(z) - inf因此波前曲率相位项为0。古依相位 ζ(z)0。 U_xy sqrt(I0) * exp(-(X.^2 Y.^2) / w0^2); % 复振幅此时相位均匀 %% 2. 计算离开束腰后的复振幅场 (例如 z zR) z_target zR; w_z w0 * sqrt(1 (z_target / zR)^2); R_z z_target * (1 (zR / z_target)^2); % 波前曲率半径 zeta_z atan(z_target / zR); % 古依相位 % 完整的复振幅计算 U_xy_z sqrt(I0) * (w0 / w_z) .* exp(-(X.^2 Y.^2) / w_z^2) .* ... exp(-1i * k * z_target) .* ... exp(-1i * k * (X.^2 Y.^2) / (2 * R_z)) .* ... exp(1i * zeta_z); %% 3. 可视化相位 phase_z angle(U_xy_z); % 提取相位范围 [-π, π] figure; subplot(1,2,1); imagesc(x*1e3, y*1e3, phase_z); axis image; colorbar; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(z %.1f mm 处的相位分布, z_target*1e3)); colormap(hsv); % HSV色图非常适合表示周期性相位 subplot(1,2,2); % 绘制通过中心的相位剖面 plot(x*1e3, phase_z(ceil(end/2), :), LineWidth, 2); xlabel(x (mm)); ylabel(相位 (弧度)); title(相位中心剖面); grid on; % 可以看到相位随x^2变化的抛物线形状这正是球面波前的特征。4.2 常见仿真问题与排查技巧在实际操作中你可能会遇到下面这些问题。这里是我的排查清单问题现象可能原因排查与解决方法图形显示一片空白或全黑/全白1. 数据值超出显示范围或全部为0/NaN。2.imagesc的坐标轴设置错误数据画到看不见的地方。3. 色图范围不合适。1. 在imagesc后使用caxis([0, 1])手动设置颜色轴范围。或用max(I_xy(:))检查数据最大值。2. 检查imagesc(x, y, I)中的x, y向量是否与矩阵I的维度匹配size(I,2)length(x),size(I,1)length(y)。3. 尝试colormap(gray)或imagesc(..., CDataMapping, scaled)。三维曲面图非常粗糙或畸变网格点太少 (linspace点数不足)。增加linspace的采样点数例如从100增加到500。对于曲面图也可以使用shading interp命令进行平滑插值。计算速度极慢1. 使用了多重嵌套for循环。2. 未预分配大数组。3. 网格分辨率过高。1.务必使用向量化操作。将exp,sqrt等函数直接作用于整个矩阵X,Y。2. 对于像I_3D这样的结果数组使用zeros预分配。3. 在保证可视效果的前提下适当降低采样点数。光束形状不对称或奇怪1. x和y的坐标范围或采样点数不一致。2. 公式中输入了错误的坐标变量。1. 确保linspace(-range, range, N)中的range和N在x和y方向一致除非你特意需要各向异性的光束。2. 检查公式中是X.^2 Y.^2而不是x.^2 y.^2。前者是矩阵后者是向量。“发散角”与理论值对不上1. 在图中测量发散角的方法不对。2.w0或λ的单位不一致。1. 理论发散角远场θ λ / (π * w0)。在仿真中可通过拟合远场w(z) ≈ θ *相位图出现剧烈跳变锯齿相位包裹问题。真实相位变化超过[-π, π]范围angle函数将其折叠到此区间。这是正常现象表示相位变化剧烈。若要分析连续的相位变化可以使用unwrap函数对一维数据效果好但对于二维相位解包裹需要更复杂的算法这属于另一个专题。可视化时用hsv色图即可清晰观察相位周期。避坑技巧在编写每个关键公式后用一两行代码做个快速验证。比如计算完w_z后fprintf(在z%.2f处w_z%.2e\n, z_target, w_z)看看输出是否合理。对于光束半径在束腰处应该等于w0在zzR处应该等于sqrt(2)*w0。这种即时的小验证能帮你快速定位是公式写错了还是参数单位错了。5. 仿真扩展与应用实例一个扎实的基础仿真框架建立后你可以像搭积木一样扩展出很多有趣且实用的应用。5.1 模拟光束通过透镜的变换这是非常经典的应用。一个薄透镜可以改变高斯光束的波前曲率从而变换其束腰位置和大小。%% 模拟透镜对高斯光束的变换 f 0.1; % 透镜焦距100mm z_from_lens 0.2; % 透镜到输入光束束腰的距离200mm % 假设入射光束参数在透镜位置处 lambda 632.8e-9; w0_in 0.5e-3; % 入射光束束腰半径 zR_in pi * w0_in^2 / lambda; % 计算入射光束在透镜位置的光斑大小w_in和波前曲率半径R_in w_in w0_in * sqrt(1 (z_from_lens / zR_in)^2); R_in z_from_lens * (1 (zR_in / z_from_lens)^2); % 应用ABCD定律薄透镜变换 % 对于薄透镜传输矩阵为 [1, 0; -1/f, 1] % 复数q参数满足 1/q2 1/q1 - 1/f % 其中 q z i*zR 或者用更常用的形式1/q 1/R - i*λ/(π*w^2) % 计算入射光束的q参数 q1_inv 1/R_in - 1i*lambda/(pi*w_in^2); % 经过透镜后的q参数 q2_inv q1_inv - 1/f; % 从q2参数提取出射光束的R_out和w_out q2 1 / q2_inv; R_out 1 / real(1/q2); w_out sqrt(-lambda / (pi * imag(1/q2))); % 进一步可以计算出射光束的新束腰位置相对于透镜和束腰大小 z_out R_out / (1 (lambda*R_out/(pi*w_out^2))^2); w0_out w_out / sqrt(1 (pi*w_out^2/(lambda*R_out))^2); fprintf(经过焦距%.3f m的透镜后\n, f); fprintf( 出射光束在透镜处的光斑半径: %.3f mm\n, w_out*1e3); fprintf( 出射光束在透镜处的波前曲率半径: %.3f m\n, R_out); fprintf( 新的束腰位置相对于透镜: %.3f m\n, z_out); fprintf( 新的束腰半径: %.3f mm\n, w0_out*1e3);这个例子展示了如何将高斯光束的q参数与ABCD矩阵理论结合用代码实现复杂光学系统的仿真。你可以改变焦距f和入射距离z_from_lens观察输出光束参数如何变化这对于光学系统设计中的光路计算非常有帮助。5.2 生成动态传播GIF动画静态图片虽然好但动态图更能直观展示传播过程。我们可以用循环生成一系列图片然后合成GIF。%% 生成光束传播动态GIF filename gaussian_beam_propagation.gif; % 沿用之前的三维数据 I_3D 或者重新计算一段更短的z轴用于动画 z_anim linspace(-2*zR, 2*zR, 50); figure(Position, [100, 100, 800, 600], Color, white); for n 1:length(z_anim) % 计算或提取第n帧的数据 (这里假设重新计算) w_z w0 * sqrt(1 (z_anim(n) / zR)^2); I_frame I0 * (w0 / w_z)^2 .* exp(-2*(X.^2 Y.^2) / w_z^2); % 绘图 imagesc(x*1e3, y*1e3, I_frame); axis image; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(高斯光束传播 | z %.2f mm, z_anim(n)*1e3)); colormap(hot); colorbar; caxis([0, 1]); % 固定颜色轴让动画稳定 drawnow; % 捕获帧并写入GIF frame getframe(gcf); im frame2im(frame); [imind, cm] rgb2ind(im, 256); if n 1 imwrite(imind, cm, filename, gif, Loopcount, inf, DelayTime, 0.1); else imwrite(imind, cm, filename, gif, WriteMode, append, DelayTime, 0.1); end end fprintf(GIF动画已保存为: %s\n, filename);这个动画能非常生动地展示光束从汇聚到发散束腰处最细的整个过程。生成的GIF可以用于报告、演示或教学材料中。5.3 像散与椭圆高斯光束仿真现实中的激光器可能输出非圆对称的即椭圆高斯光束。它的两个正交方向x和y有不同的束腰和发散角。%% 仿真椭圆高斯光束 w0x 1.0e-3; % x方向束腰半径 w0y 0.5e-3; % y方向束腰半径 zRx pi * w0x^2 / lambda; zRy pi * w0y^2 / lambda; % 计算特定z位置例如z0的光强 z_astig 0; wx_z w0x * sqrt(1 (z_astig / zRx)^2); wy_z w0y * sqrt(1 (z_astig / zRy)^2); % 光强分布为两个方向高斯函数的乘积 I_ellipse I0 * exp(-2*(X.^2)/wx_z^2 - 2*(Y.^2)/wy_z^2); figure; imagesc(x*1e3, y*1e3, I_ellipse); axis image; xlabel(x (mm)); ylabel(y (mm)); title(椭圆高斯光束 (束腰处)); colorbar; colormap(hot); % 可以看到光斑在x和y方向宽度不同呈椭圆形。通过调整w0x和w0y的比例你可以仿真出从细线到扁椭圆的各种光束形状。更进一步你还可以让两个方向的束腰位置不在同一z点这就是像散光束的典型特征仿真方法类似只是分别计算x和y方向的光束半径w(z)时使用不同的z偏移量。这些扩展实例只是抛砖引玉。基于这个核心框架你还可以去仿真高阶模如厄米-高斯或拉盖尔-高斯光束模拟光束在湍流大气中的传播或者与傅里叶光学结合计算光束的衍射。关键是理解物理模型然后用MATLAB这个强大的计算工具把它实现出来。 
