高斯光束及其MATLAB仿真

高斯光束的定义

​ 在光学中，高斯光束（英语：Gaussian beam）是横向电场以及辐照度分布近似满足高斯函数的电磁波光束。许多激光都近似满足高斯光束的条件，在这种情况中，激光在光谐振腔中以$TE{M}_{00}$波模（横向基模）传播。当它在满足近衍射极限的镜片中发生折射时，高斯光束会变换成另一种不同参数的高斯光束，因此，高斯光束是激光光学中一种方便、广泛应用的模型。
​ 描述高斯光束的数学函数是亥姆霍兹方程的一个近轴近似解（属于小角近似的一种）。这个解具有高斯函数的形式，代表了光束中电场分量的复振幅。尽管电磁波的传播包括电场和磁场两部分，研究其中任一个场，就足以描述波在传播时的性质。
​ 亥姆霍兹方程的近轴近似解可能不止一个。笛卡尔坐标系下求解可得一类称为厄米-高斯模的解，在柱坐标中求解则得到一类称为拉盖尔-高斯模的解。

基模高斯光束

一般表达式：
$$\begin{gathered} E(x,y,z)=\frac{w_0}{w}\exp (-\frac{x^2+y^2}{w^2}) \\ \exp \left\{i\left[\frac{k(x^2+y^2)}{2R}-\psi \right]\right\}\exp (ikz) \end{gathered}$$
简化表达式为：
$$E(x,y,z)=\frac{1}{1+iz/z_R}\exp \left[-\frac{(x^2+y^2)/w_0^2}{1+iz/z_R}\right]\exp (ikz)$$
$w_0$是中心束腰半径；$k=2\pi /\lambda$是波数；$\lambda$是光束的波长；$z_R=1/2k{w}_0$表示瑞利距离；$w=w_0\sqrt{1+(z/z_R{)}^2}$表示光束传播到$z$处的束宽，$R=z_R(\frac{z}{z_R}+\frac{z_R}{z})$表示等相位面的曲率半径，$\Psi$表示相位因子。

拉盖尔-高斯光束

一般表达式：
$$\begin{gathered} E_p^l(x,y,z)=\frac{w_0}{w}{\left(\sqrt{2}\frac{r}{w}\right)}^l{L}_p^l\left(-\frac{r^2}{w^2}\right)\exp \left(-\frac{r^2}{w^2}\right) \\ \exp \left\{i\left[k\frac{r^2}{2R}-(2p+l+1)\psi \right]\right\} \\ \exp (il\phi \left)\exp \right(ikz) \end{gathered}$$
简化表达式：
$$\begin{gathered} E(x,y,z)={\left(\sqrt{2}\frac{r}{w}\right)}^l{L}_p^l\left(2\frac{r^2}{w^2}\right){\left[\frac{1-iz/z_R}{\sqrt{1+(z/z_R{)}^2}}\right]}^{2p+l} \\ \frac{1}{1+iz/z_R}\exp \left[-\frac{(x^2+y^2)/w_0^2}{1+iz/z_R}\right]\exp (ikz) \end{gathered}$$

厄密-高斯光束

一般表达式：
$$\begin{gathered} E_{mn}(x,y,z)=\frac{w_0}{w}{H}_m\left(\sqrt{2}\frac{x}{w}\right)H_n\left(\sqrt{2}\frac{y}{w}\right)\exp \left(-\frac{r^2}{w^2}\right) \\ \exp \left\{i\left[k\frac{r^2}{2R}-(2p+l+1)\psi \right]\right\}\exp (ikz) \end{gathered}$$
简化表达式：
$$\begin{gathered} E_{mn}(x,y,z)=H_m\left(\sqrt{2}\frac{x}{w}\right)H_n\left(\sqrt{2}\frac{y}{w}\right){\left[\frac{1-iz/z_R}{\sqrt{1+(z/z_R{)}^2}}\right]}^{2p+l} \\ \frac{1}{1+iz/z_R}\exp \left[-\frac{(x^2+y^2)/w_0^2}{1+iz/z_R}\right]\exp (ikz) \end{gathered}$$

MATLAB仿真程序

采用MATLAB对上述高斯光束进行仿真模拟，部分子函数未编写，可以直接从其他途径搜索得到，这里只是展示主函数和最终的结果。

%============================================
clc;
close all
clear;

%% =========基本参数的输入===================
lambda0=6.328e-6;                 %基本光束的输入波长
k0=2*pi/lambda0;                   %波数
w0=1.5*lambda0;                    %高斯光束的中心束腰半径
n_i=1.0;                                   %背景空间折射率
k=k0*n_i;                                 %背景空间的波数
ZR=k*w0^2/2;                         %瑞利距离
z=0;                                          %光束传播距离
w=w0*sqrt(1+(z/ZR)^2);           %传播到z处的束宽
R=ZR*(z/ZR+ZR/z);                 %等相位面曲率半径
Phi=atan(z/ZR);                       %相位因子
m=1;
n=2;                                        %高阶高斯光束的阶数
p=1;                                         %拉盖尔高斯光束的阶数
l=1;                                          %拉盖尔高斯光束的拓扑荷数


%% =========坐标参数的设置=====================
Dx=10*lambda0;                        %x轴方向距离
Dy=10*lambda0;                        %y轴方向距离
step=0.01*lambda0;                   %步长
Numx=Dx/step+1;
Numy=Dy/step+1;                      %分割数
x=linspace(-Dx/2,Dx/2,Numx);  %x方向网格
y=linspace(-Dy/2,Dy/2,Numy);   %y方向网格
[x1,y1]=meshgrid(x,y);                 %行成网格数据
[phi,r]=cart2pol(x1,y1);                %转化成极坐标


%% ========不同光束的表达式观察平面===============
% 此处采用一般表达式进行仿真计算
%基模高斯光束表达式
E_Gauss=(1/(1+1i*z./ZR)).*exp(-(r.^2./w0^2)./(1+1i*z./ZR)).*exp(1i*k.*z);
I_Gauss=E_Gauss.*conj(E_Gauss);           %计算基模高斯光束的光强
normalizationprocessI_Gauss=I_Gauss/max(max(I_Gauss));   %光强的归一化处理

%厄米高斯光束
E_Hermite=Hermite(m,sqrt(2)*x1/w).*Hermite(n,sqrt(2)*y1/w).*((1-1i.*z/ZR)./...
                   sqrt(1+(z/ZR)^2)).^(m+n).*E_Gauss;
I_Hermite=E_Hermite.*conj(E_Hermite);     %计算厄密高斯光束的光强
normalizationprocessI_Hermite=I_Hermite/max(max(I_Hermite)); %光强的归一化处理

%拉盖尔高斯光束
E_Lague=power(sqrt(2)*r./w,l).*laguerre(p,l,2*r.^2/w.^2).*power((1-1i.*z/ZR)/...
                     sqrt(1+(z/ZR).^2),2*p+l).*E_Gauss.*exp(1i*l*phi);
I_Lague=E_Lague.*conj(E_Lague);           %计算拉盖尔高斯光束的光强
normalizationprocessI_Lague=I_Lague/max(max(I_Lague));      %光强的归一化处理


%% =========观察平面处画图设置==========================
%高斯光束图形设置
figure()                       
plot(x,normalizationprocessI_Gauss((size(x1,1))/2,:),'linewidth',1.5); %一维图像
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('Normalization Process Density','fontname','times new roman','fontsize',20); %y方向标注

figure()                                        
pcolor(x1,y1,normalizationprocessI_Gauss);                             %二维图像           
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
caxis([0,1]);
shading interp;colormap jet;colorbar; 


figure()
surf(x1,y1,normalizationprocessI_Gauss);                               %三维图像
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
caxis([0,1]);
shading interp;colormap jet;colorbar; 

%厄密高斯光束图形设置
figure()                       
plot(x,normalizationprocessI_Hermite((size(x1,1))/2,:),'linewidth',1.5); %一维图像
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('Normalization Process Density','fontname','times new roman','fontsize',20); %y方向标注

figure()                                        
pcolor(x1,y1,normalizationprocessI_Hermite);                             %二维图像           
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
caxis([0,1]);
shading interp;colormap jet;  colorbar; 

figure()
pcolor(x1,y1,angle(E_Hermite));
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
caxis([0,1]);
shading interp;  colormap jet;  colorbar;

figure()
surf(x1,y1,normalizationprocessI_Hermite);                               %三维图像
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
caxis([0,1]);
shading interp;colormap jet;colorbar; 

%拉盖尔高斯光束图形设置
figure()                       
plot(x,normalizationprocessI_Lague((size(x1,1))/2,:),'linewidth',1.5); %一维图像
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('Normalization Process Density','fontname','times new roman','fontsize',20); %y方向标注

figure()                                        
pcolor(x1,y1,normalizationprocessI_Lague);                             %二维图像           
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
caxis([0,1]);
shading interp;colormap jet;colorbar; 

figure()
pcolor(x1,y1,angle(E_Lague));
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
shading interp;colormap jet;

figure()
surf(x1,y1,normalizationprocessI_Lague);                               %三维图像
set(gca,'fontname','times new roman','fontsize',15);     %坐标轴设置
xlabel('\itx(\it\lambda)','fontname','times new roman','fontsize',20); %x方向标注
ylabel('\ity(\it\lambda)','fontname','times new roman','fontsize',20); %y方向标注
caxis([0,1]);
shading flat;colormap jet;colorbar; 
theAxes=axis;                                        %保存坐标变量


%% =======绘制动态图像=========================
%保存11组数据,用来表征传播距离变化时光强的变化
z=0:w0:10*w0;                                     %将传播距离进行离散化 
for i=1:size(z,2)
    %基模高斯光束表达式
       E_Gauss=(1/(1+1i*z(i)./ZR)).*exp(-(r.^2./w0^2)./(1+1i*z(i)./ZR)).*exp(1i*k.*z(i));
       I_Gauss=E_Gauss.*conj(E_Gauss);           %计算基模高斯光束的光强
       normalizationprocessI_Gauss(:,:,i)=I_Gauss/max(max(I_Gauss));   %光强的归一化处理
       %厄米高斯光束
        E_Hermite=Hermite(m,sqrt(2)*x1/w).*Hermite(n,sqrt(2)*y1/w).*((1-1i.*z(i)/ZR)./...
                   sqrt(1+(z(i)/ZR)^2)).^(m+n).*E_Gauss;
        I_Hermite=E_Hermite.*conj(E_Hermite);     %计算厄密高斯光束的光强
        normalizationprocessI_Hermite(:,:,i)=I_Hermite/max(max(I_Hermite)); %光强的归一化处理

        %拉盖尔高斯光束
        E_Lague=power(sqrt(2)*r./w,l).*laguerre(p,l,2*r.^2/w.^2).*power((1-1i.*z(i)/ZR)/...
                     sqrt(1+(z(i)/ZR).^2),2*p+l).*E_Gauss.*exp(1i*l*phi);
        I_Lague=E_Lague.*conj(E_Lague);           %计算拉盖尔高斯光束的光强
        normalizationprocessI_Lague(:,:,i)=I_Lague/max(max(I_Lague));      %光强的归一化处理
end