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

Matlab辅助激光光学分析与应用作者：刘良清Email：[e-mail]单位：武汉凌云光电科技有限公司毕业院校：华中科技大学激光技术与工程学院学历：硕士研究生研究方向：自适应光学、非线性光学、激光光学、固体激光器件2008年3月第一版第一章光的波动性和衍射1.1Maxwell方程组和电磁波十八世纪中叶，JamesMaxwell将已知的各种电磁作用关系用一组方程组合起来，形成了一个方程组：0ερ∇=·E(源于库伦定律的高斯定律)(1.1)0∇=·B(源于毕奥-萨瓦尔定律的高斯定律)(1.2)0t∂∇×+=∂BE(法拉第定律)(1.3)00Jt∂∇×−ε=μ∂BE(Maxwell修正的安培定律)(1.4)式中，和分别代表了电场和磁场分量。电荷密度EBρ描述路径空间单位体积内的电荷量分布；电流描述电荷的移动(单位电荷乘以速度)。J0ε表示真空介电常数，其值为。表示真空磁导率常数，其值为（或者）。122208.85410/CNm−ε=×·0μ70410/TmA−μ=π×·2k/gmC·在安培定律中引入了一个关键参数之后，Maxwell意识到，方程组构成了一个完美的电磁现象自洽理论。此外，方程组预言了电磁波的存在，并以光速传播。在Maxwell时代之前就已经有人对光速进行了测量，因此一个显而易见的结果(当时还难以令人置信)便是，光是一种高频振荡表现，类似并超越了支配电流和电荷的影响因素。而在此之前，光学还仍然作为一种独立于电学和磁学的主体进行讨论的。这里，我们不再对电磁学的基本知识进行详细的讨论，因为它们在普通物理课程中都有讲述，并且有大量的文献和书籍对其进行了细致的分析。但我们要简要的从波动方程出发，求解旁轴近似下的Maxwell方程组，得到激光传输与变换的基本方程，以方便我们后续的讨论和应用。为了体现Matlab在可视化方面的优势，我们先以一个简单的例子作为本书的开篇，以达到抛砖引玉的效果。在电动力学中，我们会遇到真空电磁场波动方程的旁轴近似解，众所周知，其解为具有高斯分布的电场复振幅：[{}220022(,)expexp()()2()()wkrrrzjjkzzwzRzwz⎡⎤Ψ=Ψ−−−−φ⎢⎥π⎣⎦](1.5)式中，为光波传播常数。、、2/=πλk()wz()Rz()φz是与光束有关的传播参数。分别表示为：2020()1zwzww⎛⎞λ=+⎜π⎝⎠⎟(1.6)1220()1⎡⎤⎛⎞π⎢⎥=+⎜⎟λ⎢⎥⎝⎠⎣⎦wRzzz(1.7)210()tan,−⎛⎞πφ==⎜⎟λ⎝⎠RRwzzZZ(1.8)光束远场发散角为：00()limzwzzw→∞λθ==π(1.9)或者220()()()wzRzwz⎛⎞⎛λθ=+⎜⎟⎜π⎝⎠⎝⎞⎟⎠(1.10)我们可以用几行简单matlab程序就可以画出具有高斯分布的电场强度，如图1.1所示，图形美观，方便对光强的分布有一个感性的视觉认识。程序代码为：clear;clc;w0=0.5;r=linspace(0,3*w0,200);eta=linspace(0,2*pi,200);[rho,theta]=meshgrid(r,eta);[x,y]=pol2cart(theta,rho);Iopt=exp(-2*rho.ˆ2/w0ˆ2);surf(x,y,Iopt);shadinginterp;xlabel(’位置/mm’);ylabel(’位置/mm’);zlabel(’相对强度/a.u.’);title(’高斯强度分布’);axis([-3*w0,3*w0,-3*w0,3*w0,0,1]);colorbar;colormap(’hot’);boxon;gridoff;图1.1高斯光强分布另外，我们还可以画出高斯光束在自由传输过程中的强度变化，如图1.2所示，程序代码如下：clear;clc;lambda=1.064e-3;w0=0.5;ZR=pi*w0ˆ2/lambda;z=linspace(-2*ZR,2*ZR,200);y=linspace(-4*w0,4*w0,200);[py,pz]=meshgrid(y,z);2wz=w0*sqrt(1+(lambda*pz/pi/w0ˆ2).ˆ2);Iopt=w0ˆ2./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统一了电磁理论以后，他马上意识到，波动可能是该方程组的解的形式。事实上，他希望找到一组满足波动形式的方程组，以辅助他完成找到真正的波动方程。既然已经知道了光是以波动方式传播的，基尔霍夫首先注意到了001/εμ正好给出了精确的光速(之前就已经被测量过)，并且法拉第和克尔已经观测到强磁场和强电场会影响光在晶体中的传播。对初始接触Maxwell方程组的人来说，并不能一眼就看出它的解具有波动形式。但是经过适当的数学操作，我们就可以将它变为波动方程的形式。83.0010/c=×ms我们来推到电场E的波动方程，磁场B的波动方程的推到过程是类似的。我们将方程(1.3)进行卷积，可得：()()0t∂∇×∇×+∇×=∂EB(1.11)该方程可以由矢量微分恒等式简化：()()2•E∇×∇×=∇∇−∇EE(1.12)卷积可由∇×B(1.4)式代换，由此得到：()2000•Jtt∂∂⎛⎞∇∇−∇+εμ+μ=⎜∂∂⎝⎠EEE0⎟(1.13)3再由(1.1)式代入上式，经过整理就可得到：2200020Jtt⎛⎞∂ρ∂∇−εμ=∇+μ⎜⎟∂ε∂⎝⎠EE(1.14)需要指出的是，上式中没有考虑到介质的极化。弱考虑到介质的极化和实际一般光学问题中自由电荷为零的条件，上式修正为：(22200002201•freeJPPttt∂∂∂∇−εμ=μ+μ−∇∇∂∂∂εEE)(1.15)式中，P为极化强度矢量。这样我们得到了一般的电场传播方程，该方程在非线性光学中有很重要的地位。当光在真空中传播时，式(1.15)中的右边所有项均为零，方程简化为：220020t∂∇−εμ=∂EE(1.16)这样我们就得到了电场传播的波动方程形式。当然在有些实际问题中，式(1.15)中右边的项并不是都为零，至少会有一项不为零，这与介质的性质有关。1.3衍射考虑一个振动频率为的光场，其复振幅可以表述为ω()jtre−ωE，则它也必须满足波动方程：()()222220jtjtnererct−ω−ω∂∇−∂EE=(1.17)由于电场振幅的含时部分是显式给出的，则方程(1.17)可以简化为：()()220rkr∇−EE=(1.18)式中是波矢量的大小。/knc≡ω(1.18)式就是所谓的赫姆霍兹方程。如果我们忽略波动的矢量特性，而只考虑它的振幅(这里不再详细讨论其过程)，那么在标量近似下，就得到了标量赫姆霍兹方程：()()220ErkEr∇−=(1.19)然后，我们考虑一束沿z轴传播的光束，它的电场复振幅写成(),,jkzExyze的形式。我们将它代入标量赫姆霍兹方程(1.19)式，得到：2222222jkzEEEEjkexyzz⎛⎞∂∂∂∂+++=⎜∂∂∂∂⎝⎠0⎟(1.20)在旁轴近似下，有222EEkzz∂∂∂∂。即是说，我们假设了电场的复振幅沿z轴传播方向是缓慢变化的，与平面波类似。但是我们允许振幅沿z轴在远大于波长量级的范围上有明显的变化。这样就得到了旁轴波方程：22222jkExyz⎛⎞∂∂∂0++≅⎜∂∂∂⎝⎠⎟(1.21)求解方程(1.21)式，得到：()()()()222,,,,0kjxxyyzjExyzExyedxdyz⎡⎤′′−+−⎢⎥⎣⎦Σ′′′≅−λ∫∫′(1.22)4于是电场的表达式为：()()()()()222,,,,,,0jkzxxyyjkzzExyzExyzejExyedxdyz⎡⎤′′−+−⎢⎥+⎢⎥⎣⎦Σ=′′′≅−λ∫∫′(1.23)值得一提的是，基尔霍夫早在1887年就提出了著名的菲涅耳-基尔霍夫衍射公式：()()1cos(,),,,,02aperturerjkrjeExyzExydxdy+⎡⎤′′′′=−∫∫rz(1.24)式(1.23)和(1.24)在分母时具有一致性，并在指数上：)⎢⎥λ⎣⎦rz≅()(222zxxyyrz′′−+−(1.25)同时，式(1.23)是(1.24)式在满足)≅+()(22zxxyy′′−+−条件下的菲尼尔旁轴近似。另外，如果进一步满足远场条件()222kxyz′′+，就得到夫琅禾费衍射近似：()()222xyjkzzjez⎛⎞+′′++⎜⎟⎜⎟⎝⎠Σ,,,,0xxyyjkzExyzExyedxdy′′′′≅−λ∫∫(1.26)1.3.1小孔衍射对称的小孔，这时，孔径上的场分布可以写为：假设光场透过一个圆柱()(),,0,0ExyE′′′=ρ(1.27)这样积分公式中，得到简化衍射积分式：，二维衍射积分可以简化为一维衍射积分。将(1.27)式代入到菲尼尔衍射()()()222cos22,,jkzjkjkzjEz⎛⎞ρ′πρρρ+⎜⎟′00zzedEedez−θ−θ⎜⎟⎝⎠ρ≅Σ′′′′−ρρρθλ∫∫(1.28)对角度的积分项，我们可以借助下面的公式完成：()2cos002jkzkdeJz′πρρ′−θ−θ′ρρ⎛⎞′θ=π⎜⎟⎝⎠∫(1.29)式中，称为零阶Bessel函数。这样，(1.28)式可以简化为：0J()()222jkzjkjkzz⎛⎞ρρ+⎜⎟⎜⎟Σ′πρρ⎛⎞⎜⎟⎝⎠(1.30)式(1.30)中的积分项也称为220,,0zzEzedEeJ⎝⎠′′′ρ≅−ρρρλ∫()22,0jkzEeρ′ρ2的汉克尔变换。在夫琅禾费衍射近似下，2jkze项ρ等于1，积分项变为的汉克尔变换。于是夫琅禾费柱对称圆孔衍射方程为：(),0E′ρ()()22,,jkzRjkzzz⎛ρ+⎜⎟⎜⎟200zEedEJ⎞⎝⎠′πρρ⎛⎞ρ≅⎝⎠31)′′′−ρρρ⎜⎟λ∫(1.5(),0E′ρ当然，虽然经过了一系列简化，然而是复振幅通常都是不确定的，即使知道了强度分布，相位分布也可能是比较难预测的。也可以通过辅助手段测量强度分布和相位分布。讨论圆这里，我们以平面波入射为例，孔夫琅禾费衍射问题。这时(),0E′ρ可用常数代替，不妨设为1。利用Bessel函数的递推关系，我们可以得到解析的圆孔夫琅禾费衍射公式：()22122jkzjkzkRJE⎛⎞⎛⎞ρρ++⎜⎟⎜⎟⎜⎟⎛⎞⎜⎟′πρρπ⎛⎞ρ2220,/zzRjkjzzedJeRzzzkRz⎜⎟⎝⎠⎝⎠ρ⎝⎠′′≅−ρρ=−⎜⎟λλρ⎝⎠∫(1.32)即使是解析式，我们还是不能直观地感受到衍射斑的样式。下面我们利用Matlab给出夫琅禾费圆孔衍射的强度分布。程序代码如下：R=0.1;ce(0,2*1.22*lambda/2/R*z,201);01);grid(r,eta);rt(theta,rho);*k).ˆ2/lambdaˆ2;,max(r),0,max(Ie(:))]);off;1,:),’k’);lambda=1.064e-3;k=2*pi/lambda;z=1.0e3;r=linspaeta=linspace(0,2*pi,2[rho,theta]=mesh[x,y]=pol2caBess=besselj(1,rho*R*k/z);Ie=4*piˆ2*Rˆ2*Bess.ˆ2./(rh