【光学】高斯光束传播模拟器 matlab代码

 ✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知，求助可私信。

🔥 内容介绍

光学领域中，高斯光束因其在激光器、光纤通信和光学测量等方面的重要应用而备受关注。精确模拟高斯光束的传播特性对于理解和预测光学系统性能至关重要。本文将深入探讨高斯光束传播模拟器的理论基础、常用的数值算法以及其在不同领域的应用前景，并分析其局限性及未来发展方向。

一、 高斯光束的理论基础

高斯光束是一种理想化的光束模型，其电场分布可以用复振幅表示，其横向强度分布服从高斯函数。描述高斯光束的关键参数包括束腰半径 ω₀，波长 λ，以及光束传播方向上的位置 z。通过这些参数，可以计算光束在任意位置 z 的束腰半径 ω(z)，曲率半径 R(z)，以及光束的Gouy相移。

具体来说，高斯光束的电场分布可以表示为：

E(r, z) = E₀ * ω₀ / ω(z) * exp[-r²/ω(z)²] * exp[-j(kz + kr²/2R(z) + ζ(z))]

其中，E₀ 为束腰处的电场振幅，r 为径向坐标，k 为波数 (k = 2π/λ)，ζ(z) 为Gouy相移。 ω(z) 和 R(z) 分别由以下公式给出：

ω(z) = ω₀ √[1 + (z/z₀)²]

R(z) = z[1 + (z₀/z)²]

z₀ = πω₀²/λ 为瑞利长度，表示光束从束腰到其面积加倍的传播距离。

上述公式准确描述了理想高斯光束的传播特性。然而，实际光束往往会受到衍射、散射、吸收等因素的影响，偏离理想高斯模型。因此，需要采用数值模拟方法来更精确地预测光束的传播行为。

二、 高斯光束传播模拟器的算法实现

高斯光束传播模拟器主要采用数值方法来求解光束传播方程，例如角谱法、分步傅里叶变换法 (FFT) 和光线追迹法等。

(1) 角谱法: 该方法基于光束的角谱表示，将光束的传播过程分解为一系列平面波的叠加。通过计算每个平面波的传播相位，可以得到光束在传播后任意位置的电场分布。角谱法具有较高的计算精度，但对于大尺寸光束或长传播距离的模拟，计算量较大。

(2) 分步傅里叶变换法 (FFT): 该方法利用傅里叶变换将光束的传播过程分解为一系列傅里叶变换和相位调制步骤。通过迭代计算，可以模拟光束在非均匀介质中的传播。FFT方法计算效率高，适合模拟大尺寸光束或长传播距离的传播情况。然而，其精度受限于采样间隔和FFT点数。

(3) 光线追迹法: 该方法基于几何光学原理，将光束表示为一系列光线，通过计算每条光线的传播路径，模拟光束的传播过程。光线追迹法计算简单，但精度较低，不适用于描述光束衍射等波动光学效应。

选择哪种算法取决于具体的应用需求。对于精度要求高的模拟，角谱法是首选；对于计算效率要求高的模拟，FFT方法更具优势；对于简单的几何光学模拟，光线追迹法则较为适用。

三、 高斯光束传播模拟器的应用前景

高斯光束传播模拟器在诸多领域具有广泛的应用前景：

(1) 激光器设计与优化: 模拟不同激光器腔型内高斯光束的传播特性，优化激光器设计参数，提高激光输出功率和光束质量。

(2) 光纤通信系统设计: 模拟光信号在光纤中的传输，分析色散、非线性效应等对光信号的影响，优化光纤通信系统性能。

(3) 光学显微镜设计与图像处理: 模拟光束在显微镜系统中的传播，分析像差、分辨率等对成像质量的影响，优化显微镜设计和图像处理算法。

(4) 自由空间光通信: 模拟光束在自由空间中的传播，分析大气湍流、光束漂移等对光通信的影响，设计更鲁棒的光通信系统。

(5) 光学传感技术: 模拟光束与待测物质的相互作用，分析光学传感器的灵敏度和精度，设计新型光学传感器。

四、 局限性与未来发展方向

尽管高斯光束传播模拟器具有广泛的应用价值，但其仍然存在一些局限性：

(1) 理想化模型: 高斯光束模型是理想化的，实际光束往往具有非高斯成分，需要更复杂的模型来描述。

(2) 计算资源: 对于复杂的模拟，计算资源消耗较大，需要更高效的算法和更强大的计算平台。

(3) 模型参数的准确性: 模拟结果的准确性依赖于模型参数的准确性，需要更精确的测量手段来获取模型参数。

未来，高斯光束传播模拟器的发展方向主要包括：

(1) 更精确的模型: 发展更精确的光束模型，例如考虑非高斯成分、光束畸变等因素。

(2) 更有效的算法: 开发更高效的数值算法，提高计算速度和精度。

(3) 更强大的计算平台: 利用高性能计算技术，例如GPU加速计算，提高模拟效率。

(4) 与实验结果的结合: 将模拟结果与实验结果进行比较，验证模型的准确性，并不断改进模型。

📣 部分代码

classdef Beam < handle

    properties

        lambda      % wavelength of the beam

        M2      % M_square factor

        beam_segments = double.empty(0, 6)      % an array storing the beam parameters of each beam segments separated by lenses

                                                % Example: beam_segments(i, :) = [beam waist diameter, waist position, divergence full angle, Rayleigh length, lens position, lens focal length]

        table = OpticTable.empty        % optic table where the beam is added

        profile = struct('z', [], 'radius', [])     % structure array storing the beam profile (i.e. beam waist radius evolution along the propagation direction z) 

        color = cell.empty(0, 2)        % color used to generate the beam profile plot

    end

    

    methods

        function self = Beam(lambda, waist, z0, M2)

            % constructor of a object of Beam class

            if nargin < 4

                self.M2 = 1;        % default M2 factor, 1 for ideal Gaussian beam

            else

                self.M2 = M2;

            end

            

            self.lambda = lambda;

            init_divg = 4 * self.M2 * self.lambda / pi / waist;

            init_zR = waist / init_divg;

            self.beam_segments(1, :) = [waist, z0, init_divg, init_zR, 0, 0];

            self.beam_color();

        end

        

        function update(self, lens_list)

            % After adding lenses, update the beam parameters of each beam segments based on the lens_lst

            self.beam_segments(2:end, :) = [];    

            if ~isempty(lens_list)

                for i = 1 : size(lens_list, 1)

                    w1 = self.beam_segments(i, 1);

                    z1 = self.beam_segments(i, 2);

                    zR1 = self.beam_segments(i, 4);

                    d1 = lens_list(i, 1) - z1;      % spacing between the beam waist and the position of the next lens

                    [w2, zR2, d2] = gaussian_focusing(w1, zR1, d1, lens_list(i, 2));        % Call the gaussian_focusing function to calculate the beam parameters after the lens

                    z2 = d2 + lens_list(i, 1);

                    divg2 = w2 / zR2;

                    self.beam_segments(i+1, :) = [w2, z2, divg2, zR2, lens_list(i, 1), lens_list(i, 2)];        % store the calculated beam parameters of the next segment into the beam_segments array

                end

            end

        end

        function profile_gen(self)

            % generate the beam propagation profile between the origin and table_end of the optic table.

            if size(self.beam_segments, 1) == 1     % when no lenses are in the beam

                z = linspace(self.table.origin, self.table.table_end, 1000);

                w_p = 0.5 * self.equation(z, 1);

            else        % when at least one lens is in the beam line

                z = linspace(self.table.origin, self.beam_segments(2, 5), 1000);

                w_p = 0.5 * self.equation(z, 1);

                

                if size(self.beam_segments, 1) > 2      % when at least two lenses are in the beam line

                    for i = 2 : size(self.beam_segments, 1) - 1

                        z_seg = linspace(self.beam_segments(i, 5), self.beam_segments(i+1, 5), 1000);

                        w_p_seg = 0.5 * self.equation(z_seg, i);

                        z = [z z_seg];

                        w_p = [w_p w_p_seg];

                    end 

                end

                z_seg = linspace(self.beam_segments(end, 5), self.table.table_end, 1000);

                w_p_seg = 0.5 * self.equation(z_seg, size(self.beam_segments, 1));

                z = [z z_seg];

                w_p = [w_p w_p_seg];

            end

            self.profile.radius = w_p;      % store the radius profile into the profile.radius array

            self.profile.z = z;     % store the propagation direction sampling into the profile.z array 

        end

        function add_to(self, optic_table)

            % add the beam to the optic_table

            if isempty(optic_table.gaussian_beam)

                optic_table.gaussian_beam = self;

                self.table = optic_table;

            else

                disp('There is already a beam on the table!') 

            end

        end

        function beam_color(self)

            % determine the beam profile color based on its wavelength

            if self.lambda > 620e-9

                self.color = {'#F44336', [244/255 67/255 54/255]};

            elseif self.lambda > 590e-9

                self.color = {'#FB8C00', [251/255 140/255 0]};

            elseif self.lambda > 570e-9

                self.color = {'#FFEB3B', [1 235/255 59/255]};

            elseif self.lambda > 495e-9

                self.color = {'#B2FF59', [178/255 255/255 89/255]};

            elseif self.lambda > 476e-9

                self.color = {'#00ACC1', [0 172/255 193/255]};

            elseif self.lambda > 450e-9

                self.color = {'#2196F3', [33/255 150/255 243/255]};

            else

                self.color = {'#AB47BC', [171/255 71/255 188/255]};

            end

        end

        function change_to(self, new_lambda, new_M2, new_z0, new_waist)

            % change the beam initial parameters

            init_waist = self.beam_segments(1, 1);

            init_z0 = self.beam_segments(1, 2);

            if nargin > 1

                self.lambda = new_lambda;

                if nargin > 2

                    self.M2 = new_M2;

                    if nargin > 3

                        init_z0 = new_z0;

                        if nargin > 4

                            init_waist = new_waist;

                        end

                    end

                end

            end

            init_divg = 4 * self.M2 * self.lambda / pi / init_waist;

            init_zR = init_waist / init_divg;

            self.beam_segments(1, :) = [init_waist, init_z0, init_divg, init_zR, 0, 0];

            self.beam_color();

            if ~isempty(self.table)         

⛳️ 运行结果

🔗 参考文献

🎈 部分理论引用网络文献，若有侵权联系博主删除

👇 关注我领取海量matlab电子书和数学建模资料

🎁  私信完整代码和数据获取及论文数模仿真定制

擅长领域：

🌈 各类智能优化算法改进及应用

生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱调度、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划（2E-VRP）、充电车辆路径规划（EVRP）、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题、港口调度、港口岸桥调度、停机位分配、机场航班调度、泄漏源定位

🌈 机器学习和深度学习时序、回归、分类、聚类和降维

2.1 bp时序、回归预测和分类

2.2 ENS声神经网络时序、回归预测和分类

2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类

2.4 CNN|TCN|GCN卷积神经网络系列时序、回归预测和分类

2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类

2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类

2.7 ELMAN递归神经网络时序、回归\预测和分类

2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类

2.9 RBF径向基神经网络时序、回归预测和分类

2.10 DBN深度置信网络时序、回归预测和分类

2.11 FNN模糊神经网络时序、回归预测

2.12 RF随机森林时序、回归预测和分类

2.13 BLS宽度学习时序、回归预测和分类

2.14 PNN脉冲神经网络分类

2.15 模糊小波神经网络预测和分类

2.16 时序、回归预测和分类

2.17 时序、回归预测预测和分类

2.18 XGBOOST集成学习时序、回归预测预测和分类

2.19 Transform各类组合时序、回归预测预测和分类

方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

🌈图像处理方面

图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

🌈 路径规划方面

旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻

🌈 无人机应用方面

无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划

🌈 通信方面

传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化、水声通信、通信上传下载分配

🌈 信号处理方面

信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化、心电信号、DOA估计、编码译码、变分模态分解、管道泄漏、滤波器、数字信号处理+传输+分析+去噪、数字信号调制、误码率、信号估计、DTMF、信号检测

🌈电力系统方面

微电网优化、无功优化、配电网重构、储能配置、有序充电、MPPT优化、家庭用电

🌈 元胞自动机方面

交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀

🌈 雷达方面

卡尔曼滤波跟踪、航迹关联、航迹融合、SOC估计、阵列优化、NLOS识别

🌈 车间调度

零等待流水车间调度问题NWFSP 、 置换流水车间调度问题PFSP、 混合流水车间调度问题HFSP 、零空闲流水车间调度问题NIFSP、分布式置换流水车间调度问题 DPFSP、阻塞流水车间调度问题BFSP

👇