Matlab辅助激光光学分析与应用

本文格式为Word版，下载可任意编辑——Matlab辅助激光光学分析与应用

Matlab激光光学

Matlab辅助激光光学分析与应用

作者：刘良清Email：llq-hust@http://.77单位：武汉凌云光电科技有限公司

毕业院校：华中科技大学激光技术与工程学院学历：硕士研究生研究方向：自适应光学、非线性光学、激光光学、固体激光

器件

2023年3月第一版

Matlab激光光学

第一章光的波动性和衍射

1.1Maxwell方程组和电磁波

十八世纪中叶，JamesMaxwell将已知的各种电磁作用关系用一组方程组合起来，形成了一个方程组：

ρ

(源于库伦定律的高斯定律)(1.1)ε0

B=0(源于毕奥-萨瓦尔定律的高斯定律)(1.2)

B

E+=0(法拉第定律)(1.3)

tEB

ε0=J(Maxwell修正的安培定律)(1.4)

μ0t

式中，E和B分别代表了电场和磁场分量。电荷密度ρ描述路径空间单位体积内的电荷量分布；电流J描述电荷的移动(单位电荷乘以速度)。ε0表示真空介电常数，其值为

E=

7

ε0=8.8541012C2/Nm2。其值为μ0=4π10Tm/A（或μ0表示真空磁导率常数，

者kgm/C）。

在安培定律中引入了一个关键参数之后，Maxwell意识到，方程组构成了一个完美的电磁现象自洽理论。此外，方程组预言了电磁波的存在，并以光速传播。在Maxwell时代之前就已经有人对光速进行了测量，因此一个显而易见的结果(当时还难以令人置信)便是，光是一种高频振荡表现，类似并超越了支配电流和电荷的影响因素。而在此之前，光学还依旧作为一种独立于电学和磁学的主体进行探讨的。

这里，我们不再对电磁学的基本知识进行详细的探讨，由于它们在普通物理课程中都有陈述，并且有大量的文献和书籍对其进行了细致的分析。但我们要简要的从波动方程出发，求解旁轴近似下的Maxwell方程组，得到激光传输与变换的基本方程，以便利我们后续的探讨和应用。

为了表达Matlab在可视化方面的优势，我们先以一个简单的例子作为本书的开篇，以达到抛砖引玉的效果。在电动力学中，我们会遇到真空电磁场波动方程的旁轴近似解，众所周知，其解为具有高斯分布的电场复振幅：

2

kr2r2

jΨ(r,z)=Ψ2exp{j[kzφ(z)]}(1.5)

2R(z)w(z)

式中，k=2π/λ为光波传播常数。w(z)、R(z)、φ(z)是与光束有关的传播参数。分别

表示为：

w(z)=w

Matlab激光光学

wz=w0*sqrt(1+(lambda*pz/pi/w02).2);Iopt=w02./wz.2.*exp(-2*py.2./wz.2);surf(pz,py,Iopt);shadinginterp;xlabel(’位置/mm’);ylabel(’位置/mm’);zlabel(’相对强度/a.u.’);title(’高斯强度分布的传输’);colorbar;

colormap(’hot’);boxon;

gridoff;

图1.2高斯光束自由传输强度变化

以上我们以简单的例子展示了Matlab在可视化方面的强大功能，但本文不再对Matlab的基本功能和语法常识进行介绍，我们认为本书的读者已经具备了基本的Matlab编程技巧。或者说，我们所做的只是将我们的实际运用跟读者进行交流探讨，促进大家共同进步。当然，我们会在一些比较关键的地方指出编程过程中需要注意的问题。

1.2波动方程

当Maxwell统一了电磁理论以后，他马上意识到，波动可能是该方程组的解的形式。事实上，他希望找到一组满足波动形式的方程组，以辅助他完成找到真正的波动方程。既然

已经知道了光是以波动方式传播的，基尔霍夫首先注意到了1/正好给出了确切的光速c=3.00108m/s(之前就已经被测量过)，并且法拉第和克尔已经观测到强磁场和强电场会影响光在晶体中的传播。对初始接触Maxwell方程组的人来说，并不能一眼就看出它的解具有波动形式。但是经过适当的数学操作，我们就可以将它变为波动方程的形式。

我们来推到电场E的波动方程，磁场B的波动方程的推到过程是类似的。我们将方程(1.3)进行卷积，可得：

(E)+

(B)=0(1.11)t

该方程可以由矢量微分恒等式简化：

卷积B可由(1.4)式代换，由此得到：

(E)=(E)2E(1.12)

(E)2E+

Eεμ+μ0J=0(1.13)00tt

Matlab激光光学

再由(1.1)式代入上式，经过整理就可得到：

2

ρ2EJ

(1.14)Eε0μ02=+μ0

tεt0

需要指出的是，上式中没有考虑到介质的极化。弱考虑到介质的极化和实际一般光学问题中

自由电荷为零的条件，上式修正为：

Jfree2E2P1

Eε0μ02=μ0+μ02(P)(1.15)

tttε0

2

式中，P为极化强度矢量。这样我们得到了一般的电场传播方程，该方程在非线性光学中有很重要的地位。当光在真空中传播时，式(1.15)中的右边所有项均为零，方程简化为：

2

2E

Eε0μ02=0(1.16)

t

这样我们就得到了电场传播的波动方程形式。当然在有些实际问题中，式(1.15)中右边的项并不是都为零，至少会有一项不为零，这与介质的性质有关。

1.3衍射

考虑一个振动频率为ω的光场，其复振幅可以表述为E(r)e动方程：

jωt

，则它也必需满足波

E(r)e

2

jωt

n22ejωt

2E(r)=0(1.17)2ct

由于电场振幅的含时部分是显式给出的，则方程(1.17)可以简化为：

式中k≡nω/c是波矢量的大小。(1.18)式就是所谓的赫姆霍兹方程。假使我们忽略波动的矢量特性，而只考虑它的振幅(这里不再详细探讨其过程)，那么在标量近似下，就得到了标量赫姆霍兹方程：

2E(r)k2E(r)=0(1.18)

2E(r)k2E(r)=0(1.19)

(x,y,z)ejkz的形式。我们然后，我们考虑一束沿z轴传播的光束，它的电场复振幅写成E

将它代入标量赫姆霍兹方程(1.19)式，得到：

2E2Ejkz2EE

+2e=0(1.20)2+2+2jk

xyzzE2E

。即是说，我们假设了电场的复振幅沿z轴传播方向是在旁轴近似下，有2k

zz2

缓慢变化的，与平面波类似。但是我们允许振幅沿z轴在远大于波长量级的范围上有明显的变化。这样就得到了旁轴波方程：

求解方程(1.21)式，得到：

22

++2jkE0(1.21)22yzx

(x,y,z)jE

λz

(x′,y′,0)e∫∫E

Σ

j

k

(xx′)2+(yy′)22z

dx′dy′(1.22)

Matlab激光光学

于是电场的表达式为：

(x,y,z)=E(x,y,z)ejkzE

j

λz

(x′,y′,0)e∫∫E

Σ

(xx′)2+(yy′)2

+zjk

2z

(1.23)

dx′dy′

值得一提的是，基尔霍夫早在1887年就提出了著名的菲涅耳-基尔霍夫衍射公式：

jejkr

E(x,y,z)=E(x′,y′,0)∫∫rλaperture

2

1+cos(r,z)′′

dxdy(1.24)2

2

式(1.23)和(1.24)在分母rz时具有一致性，并在指数上：

2z

22

同时，式(1.23)是(1.24)式在满足z(xx′)+(yy′)条件下的菲尼尔旁轴近似。另

外，假使进一步满足远场条件

xx′)+(yy′)(rz+

(1.25)

k(x′2+y′2)2

z，就得到夫琅禾费衍射近似：

j

E(x,y,z)e

λz

x2+y2jk2z+z

∫∫E(x′,y′,0)e

Σ

jk

xx′+yy′

z

′dy′(1.26)

1.3.1小孔衍射

假设光场透过一个圆柱对称的小孔，这时，孔径上的场分布可以写为：

得到简化衍射积分式：

这样，二维衍射积分可以简化为一维衍射积分。将(1.27)式代入到菲尼尔衍射积分公式中，

ρ2jk+z2z

ρ22πρρ′

jkcos(θθ′)jk

z2z

E(x′,y′,0)=E(ρ′,0)(1.27)

E(ρ,z)

j

eλz

2π

∫ρ′dρ′E(ρ′,0)e

Σ

∫dθ′e

(1.28)

对角度的积分项，我们可以借助下面的公式完成：

∫dθ′e

ρρ′

(θθ′)jkz

kρρ′=2πJ0(1.29)

z

ρ22z

式中，J0称为零阶Bessel函数。这样，(1.28)式可以简化为：

2πj

E(ρ,z)e

λz

ρ2

jkz2z

∫ρ′dρ′E(ρ′,0)e

Σ

jk

kρρ′

J0(1.30)z

jkρ2

2z

式(1.30)中的积分项也称为E(ρ′,0)e

jk

ρ2z

2

的汉克尔变换。在夫琅禾费衍射近似下，e

项

等于1，积分项变为E(ρ′,0)的汉克尔变换。于是夫琅禾费柱对称圆孔衍射方程为：

E(ρ,z)

2πjeλz

ρ2jkz2z

kρρ′′′′ρρρdE0J,()0(1.31)∫zR

Matlab激光光学

1.4波前畸变

在前面的计算与仿真中，好多时候都用到了平面波假设，然而实际光束会经过各种光学元件，即使是大气传输也会导致光束的波前发生畸变。波前畸变导致光束质量下降，光束可聚焦能力降低，在激光系统中会严重影响系统的稳定性和可靠性。

1.4.1Zernike多项式

Zernike多项式是定义在单位圆上正交多项式：

mm

NR(r)cosmθ(m≥0)nnm

(1.40)Zn(r,θ)=mm

NnRn(r)sinmθ(m0)

式中，n和m分别为径向自由度和角向频率，m的变化范围是:-n,-n+2,…,n-2,n。其中，径向多项式定义为：

m

m

Rn(r)=

(nm)/2

s=0

(1)(ns)!rn2s

s!n+m/2s!nm/2s!

s

(1.41)

式中Nn为归一化因子，表述为：

m

=Nn

(1.42)

δm,0的条件为，当m=0时，δm,0=1；反之，δm,0=0。这样一来，我们可以用Zernike

多项式来拟合定义在一个归一化单位圆内的畸变表面W(r,θ)：

m

W(r,θ)=∑∑cn,mZn(r/rmax,θ)(1.43)

n=0mN

式中，N表示用于拟合描述的最大阶数，cn,m为Zernike系数，它由下式决定之：

cn,m=∫

rmax

∫

2π

m

W(r,θ)Zn(r/rmax,θ)rdrdθ(1.44)

下面我们看看Zernike多项式前几阶的图形结构。如图1.10所示，为前六阶Zernike模式

结构，它们分别代表了像差理论中的倾斜、离焦和像散等波前畸变。程序代码如下：

functionvarargout=zernikepol(num,hang,lie,str)%使用方法

%num为依照定义顺序给出的Zernike阶数,可以只输入该参数%hang和lie分别定义径向和角向网格划分数%str为指定是圆域’circ’或是方形域’square’%例如：zernikepol(3)

%[x,y,Z]=zernikepol(3);可将网格坐标输出ifnargin4

error(’输入参数过多！’);end

ifnargin4str=’circ’;ifnargin3lie=200;ifnargin2hang=200;

Matlab激光光学

fordf=1:mr2form=1:ne0