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

M atlab辅助激光光学分析与应用 2008 年 3 月 第一版 Matlab 辅助激光光学分析与应用 II MMaattllaabb辅助激光光学分析与应用.I 第一章 光的波动性和衍射.1 1.1 Maxwell方程组和电磁波.1 1.2 波动方程.3 1.3 衍射.4 1.3.1 小孔衍射.5 1.3.2 双缝衍射.11 1.4 波前畸变.13 1.4.1 Zernike 多项式.13 1.4.2 畸变光束的衍射.15 1.5 光束通过光学元件的变换.18 1.5.1 平行光束通过透镜的聚焦.18 1.5.2 高斯光束通过透镜的聚焦.23 1.5.3 自聚焦透镜.25 1.6 高斯光束.27 1.6.1 高阶高斯光束.27 1.6.2 高斯光束的传输变换.31 第二章 激光谐振腔.35 2.1 激光谐振腔的本征方程.35 2.2 无源腔的 Fox-Li 迭代方法.35 2.3 无源腔的矩阵特征向量方法.40 2.3.1 平行平面腔.40 2.3.2 双凹腔.43 2.4 基模谐振腔.49 Matlab 辅助激光光学分析与应用 1 第一章 光的波动性和衍射 1.1 Maxwell方程组和电磁波 十八世纪中叶，James Maxwell将已知的各种电磁作用关系用一组方程组合起来，形成了一个方程组：0 E (源于库伦定律的高斯定律)(1.1)0 B (源于毕奥-萨瓦尔定律的高斯定律)(1.2)0BEt (法拉第定律)(1.3)00BEJt (Maxwell 修正的安培定律)(1.4)式中，E和B分别代表了电场和磁场分量。电荷密度描述路径空间单位体积内的电荷量分布；电流J描述电荷的移动(单位电荷乘以速度)。0表示真空介电常数，其值为122208.854 10/CN m。0表示真空磁导率常数，其值为70410/T mA（或者2k/g mC）。在安培定律中引入了一个关键参数之后，Maxwell 意识到，方程组构成了一个完美的电磁现象自洽理论。此外，方程组预言了电磁波的存在，并以光速传播。在 Maxwell 时代之前就已经有人对光速进行了测量，因此一个显而易见的结果(当时还难以令人置信)便是，光是一种高频振荡表现，类似并超越了支配电流和电荷的影响因素。而在此之前，光学还仍然作为一种独立于电学和磁学的主体进行讨论的。这里，我们不再对电磁学的基本知识进行详细的讨论，因为它们在普通物理课程中都有讲述，并且有大量的文献和书籍对其进行了细致的分析。但我们要简要的从波动方程出发，求解旁轴近似下的 Maxwell 方程组，得到激光传输与变换的基本方程，以方便我们后续的讨论和应用。为了体现 Matlab 在可视化方面的优势，我们先以一个简单的例子作为本书的开篇，以达到抛砖引玉的效果。在电动力学中，我们会遇到真空电磁场波动方程的旁轴近似解，众所周知，其解为具有高斯分布的电场复振幅：220022(,)expexp()()2()()wkrrr zjj kzzw zR zwz(1.5)式中，2/k为光波传播常数。()w z、()R z、()z是与光束有关的传播参数。分别表示为：Matlab 辅助激光光学分析与应用 2 2020()1zw zww(1.6)220()1wR zzz(1.7)210()tan,RRwzzZZ(1.8)光束远场发散角为：00()limzw zzw(1.9)或者 220()()()w zR zw z(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/w02);surf(x,y,Iopt);shading interp;xlabel(位置/mm);ylabel(位置/mm);zlabel(相对强度/a.u.);title(高斯强度分布);axis(-3*w0,3*w0,-3*w0,3*w0,0,1);colorbar;colormap(hot);box on;grid off;图 1.1 高斯光强分布 另外，我们还可以画出高斯光束在自由传输过程中的强度变化，如图 1.2 所示，程序代码如下：clear;Matlab 辅助激光光学分析与应用 3 clc;lambda=1.064e-3;w0=0.5;ZR=pi*w02/lambda;z=linspace(-2*ZR,2*ZR,200);y=linspace(-4*w0,4*w0,200);py,pz=meshgrid(y,z);wz=w0*sqrt(1+(lambda*pz/pi/w02).2);Iopt=w02./wz.2.*exp(-2*py.2./wz.2);surf(pz,py,Iopt);shading interp;xlabel(位置/mm);ylabel(位置/mm);zlabel(相对强度/a.u.);title(高斯强度分布的传输);colorbar;colormap(hot);box on;grid off;图 1.2 高斯光束自由传输强度变化 以上我们以简单的例子展示了 Matlab 在可视化方面的强大功能，但本文不再对 Matlab的基本功能和语法常识进行介绍，我们认为本书的读者已经具备了基本的 Matlab 编程技巧。或者说，我们所做的只是将我们的实际运用跟读者进行交流讨论，促进大家共同进步。当然，我们会在一些比较关键的地方指出编程过程中需要注意的问题。1.2 波动方程 当 Maxwell 统一了电磁理论以后，他马上意识到，波动可能是该方程组的解的形式。事实上，他希望找到一组满足波动形式的方程组，以辅助他完成找到真正的波动方程。既然已经知道了光是以波动方式传播的，基尔霍夫首先注意到了001/正好给出了精确的光速83.00 10/cms(之前就已经被测量过)，并且法拉第和克尔已经观测到强磁场和强电场会影响光在晶体中的传播。对初始接触 Maxwell 方程组的人来说，并不能一眼就看出它的解具有波动形式。但是经过适当的数学操作，我们就可以将它变为波动方程的形式。我们来推到电场 E 的波动方程，磁场 B 的波动方程的推到过程是类似的。我们将方程(1.3)进行卷积，可得：0EBt(1.11)Matlab 辅助激光光学分析与应用 4 该方程可以由矢量微分恒等式简化：2EEE(1.12)卷积B可由(1.4)式代换，由此得到：20000EEEJtt(1.13)再由(1.1)式代入上式，经过整理就可得到：2200020EEJtt(1.14)需要指出的是，上式中没有考虑到介质的极化。若考虑到介质的极化和实际一般光学问题中自由电荷为零的条件，上式修正为：22200002201EEfreeJPPttt(1.15)式中，P为极化强度矢量。这样我们得到了一般的电场传播方程，该方程在非线性光学中有很重要的地位。当光在真空中传播时，式(1.15)中的右边所有项均为零，方程简化为：220020EEt(1.16)这样我们就得到了电场传播的波动方程形式。当然在有些实际问题中，式(1.15)中右边的项并不是都为零，至少会有一项不为零，这与介质的性质有关。1.3 衍射 考虑一个振动频率为的光场，其复振幅可以表述为Ej tr e，则它也必须满足波动方程：222220EEj tj tner erct(1.17)由于(是假设)电场振幅的含时部分是显式给出的，则方程(1.17)可以简化为：220EErkr(1.18)式中/knc是波矢量的大小。(1.18)式就是所谓的赫姆霍兹方程。如果我们忽略波动的矢量特性，而只考虑它的振幅(这里不再详细讨论其过程)，那么在标量近似下，就得到了标量赫姆霍兹方程：220E rk E r(1.19)然后，我们考虑一束沿 z 轴传播的光束，它的电场复振幅写成,jkzE x y z e的形式。我们将它代入标量赫姆霍兹方程(1.19)式，得到：22222220jkzEEEEjkexyzz(1.20)Matlab 辅助激光光学分析与应用 5 在旁轴近似下，有222EEkzz。即是说，我们假设了电场的复振幅沿 z 轴传播方向是缓慢变化的，与平面波类似。但是我们允许振幅沿 z 轴在远大于波长量级的范围上有明显的变化。这样就得到了旁轴波方程：222220jkExyz(1.21)求解方程(1.21)式，得到：222,0kjxxyyzjE x y zE x yedx dyz(1.22)于是电场的表达式为：222,0jkzxxyyjkzzE x y zE x y z ejE x yedx dyz(1.23)值得一提的是，基尔霍夫早在 1887 年就提出了著名的菲涅耳-基尔霍夫衍射公式：1cos(,),02r zjkraperturejeE x y zE x ydx dyr(1.24)式(1.23)和(1.24)在分母rz时具有一致性，并在指数上：222xxyyrzz(1.25)同时，式(1.23)是(1.24)式在满足22zxxyy条件下的菲尼尔旁轴近似。另外，如果进一步满足远场条件222k xyz，就得到夫琅禾费衍射近似：222,0 xyxxyyjkzjkzzjE x y zeE x yedx dyz(1.26)1.3.1 小孔衍射 假设光场透过一个圆柱对称的小孔，这时，孔径上的场分布可以写为：,0,0E x yE(1.27)这样，二维衍射积分可以简化为一维衍射积分。将(1.27)式代入到菲尼尔衍射积分公式中，得到简化衍射积分式：222cos220,0jkzjkjkzzzjEzedEed ez(1.28)对角度的积分项，我们可以借助下面的公式完成：Matlab 辅助激光光学分析与应用 6 2cos002jkzkd eJz(1.29)式中，0J称为零阶 Bessel 函数。这样，(1.28)式可以简化为：222202,0jkzjkzzjkEzedEeJzz(1.30)式(1.30)中的积分项也称为22,0jkzEe的汉克尔变换。在夫琅禾费衍射近似下，22jkze项等于 1，积分项变为,0E的汉克尔变换。于是夫琅禾费柱对称圆孔衍射方程为：2202,0jkzzRjkEzedEJzz(1.31)虽然经过了一系列简化，然而,0E是复振幅通常都是不确定的，即使知道了强度分布，相位分布也可能是比较难预测的。当然，也可以通过辅助手段测量强度分布和相位分布。这里，我们以平面波入射为例，讨论圆孔夫琅禾费衍射问题。这时,0E可用常数代替，不妨设为 1。利用 Bessel 函数的递推关系，我们可以得到解析的圆孔夫琅禾费衍射公式：221222022,/jkzjkzzzRk RJjkjzEzedJeRzzzk Rz(1.32)-15-10-505101500.10.20.30.40.50.60.70.80.91x 10 图 1.3 夫琅禾费圆孔衍射(孔径 0.2mm，距离 1m 远)Matlab 辅助激光光学分析与应用 7 即使是解析式，我们还是不能直观地感受到衍射斑的样式。下面我们利用 Matlab 给出夫琅禾费圆孔衍射的强度分布。程序代码如下：R=0.1;lambda=1.064e-3;k=2*pi/lambda;z=1.0e3;r=linspace(0,2*1.22*lambda/2/R*z,201);eta=linspace(0,2*pi,201);rho,theta=meshgrid(r,eta);x,y=pol2cart(theta,rho);Bess=besselj(1,rho*R*k/z);Ie=4*pi2*R2*Bess.2./(rho*k).2/lambda2;surf(x,y,Ie);axis(-max(r),max(r),-max(r),max(r),0,max(Ie(:);shading interp;box on;grid off;figure;plot(x(1,:),Ie(1,:),k,x(101,:),Ie(101,:),k);程序中使用了这样一个参数1.22/D，因为平面波假设具有最小的衍射角，这个参数就称为衍射极限角，所以在衍射区域里，只取二倍衍射极限范围就可以大致看出小孔的衍射特性。并且，对于解析表达，Matlab 提供了 Bessel 函数工具箱，可以直接调用。当然，为了能够从图上看到跟实际观测相近的衍射环，可以将相对强度分布取四次开根号，如图 1.3下图所示。-15-10-505101500.10.20.30.40.50.60.70.80.91x 10-3位 置 /mm相对强度 /a.u.图 1.4 数值积分法计算的夫琅禾费圆孔衍射 上面，我们使用的解析表达式绘图。接下来，我们采用数值积分方法直接求解(1.26)式夫琅禾费圆孔衍射分布。这需要将目标平面划分网格，然后用网格格点上的复振幅代替领域内的平均振幅分布求相面上的振幅分布。为了简化计算量，我们不再计算二维分布，只计算一维分布，计算结果如图 1.4 所示。程序代码如下：R=0.1;lambda=1.064e-3;k=2*pi/lambda;z=1.0e3;r=linspace(0,2*1.22*lambda/2/R*z,201);eta=linspace(0,2*pi,201);rho,theta=meshgrid(r,eta);x,y=pol2cart(theta,rho);r0=linspace(0,R,201);eta0=linspace(0,2*pi,201);rho0,theta0=meshgrid(r0,eta0);x0,y0=pol2cart(theta0,rho0);deta=R/200*2*pi/200;E2=zeros(201,1);for gk=1:201 Matlab 辅助激光光学分析与应用 8 for m=1:200 for n=1:201 E2(gk)=E2(gk)-j/lambda/z*exp(x(1,gk)2+y(1,gk)2)/z/2+z)*j*k)*exp(j*k*(x(1,gk)*x0(m,n)+y(1,gk)*y0(m,n)/z)*deta*rho0(m,n);end end end Ie=conj(E2).*E2;plot(r,Ie,k,-r,Ie,k);比较图 1.3 和图 1.4，可以看出，数值积分法计算出的衍射斑与理论解析结果是一致的。另外，值得注意的是，从图上看到，一级衍射峰值比零级衍射峰值要低很多，其他级次就更低了。那为什么我们在实验中仍然可以清晰地看到许多级次的衍射环呢？那是因为人眼对光子的感光度比较高，而实验所用的激光都具有很高的光子数密度，因此还是可以看到很多的衍射级次。对于旁轴远场条件不够满足的情况，夫琅禾费近似不成立，即使是菲涅耳近似也可能不成立，这时必须使用更为严格的菲涅耳-基尔霍夫标量衍射公式：1cos,02jkrjE x y zE x yedsr(1.33)式中，222rxxyyz。我们仍然以平面波的圆孔衍射为例，对(1.33)式数值计算不同相面上的衍射图样，如图 1.5 所示。程序代码如下：%exmp1_3_1 R=0.1;lambda=1.064e-3;k=2*pi/lambda;for z=0.5:5:25;r=linspace(0,2*1.22*lambda/2/R*z,201);eta=linspace(0,2*pi,201);rho,theta=meshgrid(r,eta);x,y=pol2cart(theta,rho);r0=linspace(0,R,201);eta0=linspace(0,2*pi,201);rho0,theta0=meshgrid(r0,eta0);x0,y0=pol2cart(theta0,rho0);deta=R/200*2*pi/200;E2=zeros(201,1);for gk=1:201 for m=1:200 for n=1:201 Rrho=sqrt(x(1,gk)-x0(m,n)2+(y(1,gk)-y0(m,n)2+z2);Rtheta=z/Rrho;E2(gk)=E2(gk)-j/lambda/2*exp(Rrho*j*k)*(1+Rtheta)/Rrho*deta*rho0(m,n);end end end Ie=conj(E2).*E2;Ie=Ie/max(Ie);plot3(z*ones(size(r),r,Ie,k,z*ones(size(r),-r,Ie,k);hold on;end 0510152025-0.4-0.200.20.400.51propagating derection /mmMatlab 辅助激光光学分析与应用 9 图 1.5 不同传输距离处的圆孔衍射样式(孔径 0.2mm)可见，在远场条件不满的情况下，衍射样式有较大的差异。远场条件由菲涅耳数决定：2aFL(1.34)式中，a 为对光束起实际限制作用元件的横向限度(半径)，L 为传输距离。当 F 满足约为 1的数量级时，菲涅耳衍射近似成立，当 F 满足远小于 1 的数量级时，夫琅禾费衍射近似成立。(a)(b)图 1.6 矩孔衍射图样 D=0.2mm(a)z=20,(b)z=500 当衍射孔为矩形孔时，同样利用(1.33)式，可以数值求解不同相面上的衍射样式。如图 1.6 所示。其中左图为衍射样式的相对强度分布，右图为了突出与实际观测的近似，对强度分布取四次开根号。程序代码如下：%exmp1_3_2 矩孔衍射 clear;tic;R=0.1;lambda=1.064e-3;k=2*pi/lambda;z=500;xmax=2*1.22*lambda/2/R*z;x=linspace(-xmax,xmax,61);y=x;x,y=meshgrid(x,y);x0=linspace(-R,R,61);y0=x0;x0,y0=meshgrid(x0,y0);deta=(2*R/60)2;E2=zeros(61,61);for bk=1:61 for gk=1:61 for m=1:60 for n=1:60 Matlab 辅助激光光学分析与应用 10 Rrho=sqrt(x(gk,bk)-x0(m,n)2+(y(gk,bk)-y0(m,n)2+z2);Rtheta=z/Rrho;E2(gk,bk)=E2(gk,bk)-j/lambda/2*exp(Rrho*j*k)*(1+Rtheta)/Rrho*deta;end end end end Ie=conj(E2).*E2;surf(x,y,Ie);shading interp;axis(-xmax,xmax,-xmax,xmax);figure;Ie=conj(E2).*E2;surf(x,y,Ie.(1/4);shading interp;axis(-xmax,xmax,-xmax,xmax);toc;另外，矩孔衍射也可以求得平面波衍射夫琅禾费近似解析表达式：222222,sin/2sin/2/2/2xyxxyyjkzjkzzxyjkzzjAE x y zeedx dyzkxazkybzjAezkxzkyz(1.35)则远场衍射的强度分布为：22222222sin/2sin/2/2/2FkxazkybzAIzkxzkyz(1.36)式中，a 和 b 分别为矩孔的边长。它的衍射图样如图 1.7 所示，单缝衍射与此有相类似的分布，只不过去除一个自由度，这里就不再赘述。程序代码如下：%矩孔衍射的解析计算 exmp1_3_3 R=0.1;lambda=1.064e-3;k=2*pi/lambda;z=1.0e3;xmax=8*1.22*lambda/2/R*z;x=linspace(-xmax,xmax,200);y=x;x,y=meshgrid(x,y);IF=sin(k*x*R/z).2.*sin(k*y*R/z).2./(k*x/2/z).2./(k*y/2/z).2/lambda2/z2;surf(x,y,IF.(1/2);colormap(hot)axis equal shading interp;Matlab 辅助激光光学分析与应用 11 00.010.020.030.04-60-40-200204060 图 1.7 夫琅禾费矩孔衍射振幅绝对值图样 1.3.2 双缝衍射 对于双缝夫琅禾费衍射，根据相干叠加原理，只需要将(1.35)稍作调整并叠加：22/22/22sin/2/2,/2/2sin/2/2/2/2xdjkzzxdjkzzk xdazjAE x y zezk xdazk xdazjAezk xdaz(1.37)式中，a 为缝宽，d 为缝的中心间距。计算结果如图 1.8 所示。程序代码如下：%双缝衍射解析计算 exmp1_3_4 clear;R=0.1;d=4*R;lambda=1.064e-3;k=2*pi/lambda;z=1.0e3;xmax=3*1.22*lambda/2/R*z;x=linspace(-xmax,xmax,200);y=x;x,y=meshgrid(x,y);IF1=j*exp(j*k*(x-d/2).2/2/z+z).*sin(k*(x-d/2)*R/z)./(k*(x-d/2)/2/z)/lambda/z;IF2=j*exp(j*k*(x+d/2).2/2/z+z).*sin(k*(x+d/2)*R/z)./(k*(x+d/2)/2/z)/lambda/z;IF=IF1+IF2;surf(y,x,abs(IF);colormap(hot)shading interp;figure;plot(abs(IF(100,:).2,x(100,:);Matlab 辅助激光光学分析与应用 12 00.020.040.060.080.10.120.14-20-15-10-505101520 图 1.8 夫琅禾费双缝衍射(孔径=0.2mm，缝间距 d=0.4mm)如果不满足夫琅禾费远场条件或者入射光场不是平面波，则不能直接求取目标平面上的衍射场的解析表达式。这时，也只能使用菲涅耳-基尔霍夫标量衍射公式(1.33)式。例如，假设入射到两狭缝的光场振幅均为超高斯分布：00expnsxEAw(1.38)当指数项中的 n2 时，称为超高斯分布，且 n 越大越接近平顶分布。同时，双缝衍射可以减少一个自由度，积分公式变为：1cos,02jkrDjE x zE xedxr(1.39)式中22rxxz，cos/z r。设缝宽 a=0.2mm，sw=0.08mm，n=4,缝的间距d=0.4mm，入射场的传播轴与缝的中心平行。计算过节如图 1.9 所示，程序代码如下：%exmp1_3_5 双缝衍射数值解 clear;tic;R=0.1;d=0.2;ws=0.08;lambda=1.064e-3;k=2*pi/lambda;z=1.0e3;xmax=3*1.22*lambda/2/R*z;x=linspace(-xmax,xmax,61);y=x;x,y=meshgrid(x,y);x0=linspace(-R-d,R+d,61);y0=x0;x0,y0=meshgrid(x0,y0);E0=exp(-(x0-d).4/ws4)+exp(-(x0+d).4/ws4);deta=(2*R/60)2;E2=zeros(61,61);for bk=1:61 for gk=1:61 for m=1:60 for n=1:60 Rrho=sqrt(x(gk,bk)-x0(m,n)2+z2);Rtheta=z/Rrho;E2(gk,bk)=E2(gk,bk)-j/lambda/2*E0(m,n)*exp(Rrho*j*k)*(1+Rtheta)/Rrho*deta;end end Matlab 辅助激光光学分析与应用 13 end end Ie=conj(E2).*E2;surf(x,y,Ie);shading interp;axis(-xmax,xmax,-xmax,xmax);figure;Ie=conj(E2).*E2;surf(x,y,Ie.(1/2);shading interp;axis(-xmax,xmax,-xmax,xmax);toc;图 1.9 超高斯光束双缝衍射，左图为相对强度，右图为振幅绝对值 (缝宽 D=0.2mm 缝间距 d=0.4mm,距离 z=1m)虽然可以减少一个自由度，但上面的程序仍然使用了二维数值积分，这是为了画出二维强度分布图。如果只是关注一维强度分布的话，则程序可以再简化，以减少计算量，节约程序运行时间。原则上，我们可以模拟任意复振幅分布入射的双缝衍射，甚至其他几何孔径的衍射，只要确定知道了复振幅分布的具体数值。例子很多，这里就不再一一演示了。1.4 波前畸变 在前面的计算与仿真中，很多时候都用到了平面波假设，然而实际光束会经过各种光学元件，即使是大气传输也会导致光束的波前发生畸变。波前畸变导致光束质量下降，光束可聚焦能力降低，在激光系统中会严重影响系统的稳定性和可靠性。1.4.1 Zernike 多项式 Zernike 多项式是定义在单位圆上正交多项式：cos0,sin0mmnnmnmmnnNRrmmZrNRrmm(1.40)式中，n 和 m 分别为径向自由度和角向频率，m 的变化范围是:-n,-n+2,n-2,n。其中，径向多项式定义为：/2201!/2!/2!snmnsmnsns rRrsnmsnms(1.41)Matlab 辅助激光光学分析与应用 14 式中mnN为归一化因子，表述为：,0211mnmnN(1.42),0m的条件为，当 m=0 时，,0m=1；反之，,0m=0。这样一来，我们可以用 Zernike 多项式来拟合定义在一个归一化单位圆内的畸变表面,Wr：,max0,/,Nmn mnnmWrcZr r(1.43)式中，N 表示用于拟合描述的最大阶数，,n mc为 Zernike 系数，它由下式决定之：max2,max00,/,rmn mncWrZr rrdrd(1.44)下面我们看看 Zernike 多项式前几阶的图形结构。如图 1.10 所示，为前六阶 Zernike 模式结构，它们分别代表了像差理论中的倾斜、离焦和像散等波前畸变。程序代码如下：function varargout=zernikepol(num,hang,lie,str)%num 为按照定义顺序给出的 Zernike 阶数,可以只输入该参数%hang 和 lie 分别定义径向和角向网格划分数，str 为指定是圆域circ或是方形域square%示例：zernikepol(3)，x,y,Z=zernikepol(3);可将网格坐标输出 if nargin4 error(输入参数过多！);end if nargin4 str=circ;if nargin3 lie=200;if nargin2 hang=200;lie=hang;if nargin1 error(参数不足！);end end end end if num0)：);end if strcmp(str,circ)rad=linspace(0,1,hang);ang=linspace(0,2*pi,lie);rho,theta=meshgrid(rad,ang);x,y=pol2cart(theta,rho);elseif strcmp(str,square)lon=sqrt(2)/2;x,y=meshgrid(linspace(-lon,lon,hang),linspace(-lon,lon,lie);theta,rho=cart2pol(x,y);end cont=0;flag=0;n=1;while contnum for m=-n:2:n cont=cont+1;if cont=num flag=1;break;end end if flag=1 Matlab 辅助激光光学分析与应用 15 break;end n=n+1;end if m=0 deltam=1;else deltam=0;end Norm=sqrt(2*(n+1)/(1+deltam);Rnm_rho=zeros(size(rho);for s=0:(n-abs(m)/2 Rnm_rho=Rnm_rho+(-1)s.*prod(1:(n-s)*rho.(n-2*s)/(prod(1:s)*prod(1:(n+abs(m)/2-s)*prod(1:(n-abs(m)/2-s);end if m10 hmax=10;else hmax=Dr;end h0=linspace(-hmax,hmax,mh);mz=1000;z0=20;%初始光线与透镜平面的距离 y=zeros(size(z0);theta1=asin(hmax/R);theta2=asin(n*hmax/R);theta=theta2-theta1;f=hmax/tan(theta);%透镜的近似焦距 z=linspace(0,f+z0+f/3,mz);figure;for gh=1:mh theta1=asin(h0(gh)/R);theta2=asin(n*h0(gh)/R);theta=theta2-theta1;for gz=1:mz L=sqrt(R2-h0(gh)2)-(R-d);if z(gz)lambda/10 mr0=mr0+1;end ne0=mr0;Matlab 辅助激光光学分析与应用 21 rmax=5*f*phy;r=linspace(0,rmax,mr2);eta=linspace(0,2*pi,ne2);rho,theta=meshgrid(r,eta);x,y=pol2cart(theta,rho);r0=linspace(0,R0,mr0);eta0=linspace(0,2*pi*(ne0-1),ne0-1);rho0,theta0=meshgrid(r0,eta0);x0,y0=pol2cart(theta0,rho0);deta=R0/(mr0-1)*2*pi/(ne0-1);E2=zeros(size(x);for gk=1:ne2 for df=1:mr2 Rrho=sqrt(x(gk,df)-x0).2+(y(gk,df)-y0).2+z2);Rtheta=z./Rrho;opd=exp(j*k*(n-1)*(sqrt(RL2-rho0.2)-(RL-d)+d);Ep=-j/lambda/2*exp(Rrho*j*k).*(1+Rtheta)./Rrho*deta.*rho0.*opd;E2(gk,df)=sum(Ep(:);end end Ie=conj(E2).*E2;%Ie=Ie/max(Ie(:);figure;surf(x,y,Ie);shading interp;axis(-rmax,rmax,-rmax,rmax)grid off;box on;toc;-0.1-0.08-0.06-0.04-0.0200.020.040.060.080.105001000150020002500300035004000 图 1.15 平面波经过透镜的聚焦(入射光斑半径 1mm 平凸透镜凸面曲率半径 25mm，中心厚度 3mm)有必要指出的是，随着焦距的缩短，对初始光场的网格划分要求会不同，这是为了保证衍射积分的精度，尽量减少计算的误差。网格划分的要求为，任意两个格点对衍射场点的光程差应该不大于十分之一波长。否者，相位起伏的跳变会造成计算结果出现较大的偏差，甚至得出错误的计算结果。当然，这首先也要保证初始光场两个格点之间的相位不能有较大的跳变。而衍射场的网格划分则没有具体的要求，只要能保证可以较为连续地进行曲线拟合就行，这个过程是由画图函数自己完成的。但是衍射场的坐标范围需要仔细地设计，保证可以看到完整的光斑大小即可。而聚焦光斑的大小是由入射光束的发散角和透镜的焦距决定的，通常可以由下式估计之：Df(1.51)式中 D 为聚焦光斑的直径，f 为透镜的几何焦距，为入射光束的远场发散角(全角)。由于一般情况下，入射光束均为高斯光束，所以其发散角或者通过测量得到，或者由入射光束的束腰大小估计理想高斯光束的远场发散角。这个问题，我们后面还会有所讨论，这里不再赘述。如图 1.16 所示，当透镜的焦距减小时(透镜凸面曲率半径为 15mm)，聚焦光斑也变小，这时光斑半径大约为 0.012mm。同时，我们可以看见，聚焦光斑的强度分布与小孔远场衍射有一致的极大和极小点，说明透镜具有远场变换性质。Matlab 辅助激光光学分析与应用 22 -0.06-0.04-0.0200.020.040.06010002000300040005000600070008000900010000 图 1.16 平面波经过透镜的聚焦(入射光斑半径 1mm 平凸透镜凸面曲率半径 15mm，中心厚度 3mm)当平面波的波前存在畸变，或者由于透镜本身存在像差的时候，聚焦光斑的形状也将随之发生改变：222 1,1cos,2jknRxydjkx yjkrjE x y zeeedsr(1.52)式中，,x y为波前附加的相位畸变项。这里，只举两个 Zernike 波前畸变的例子，图 1.17所示分别为波前存像散和彗差时，透镜几何焦面上的衍射光强分布。其结果与畸变光束的圆孔远场衍射是一致的。程序代码如下：%exmp1_5_3 畸变平面波透镜焦面衍射数值计算 clear;clc;tic;n=1.5062;%1064nm 波长折射率，k9 玻璃 d=3;%透镜中心厚度 RL=0.025e3;%透镜凸面曲率半径 f=RL/(n-1);%透镜的焦距 R0=1;%入射光束半径 lambda=1.064e-3;k=2*pi/lambda;phy=lambda/pi/R0;z=f;mr2=51;ne2=61;mr0=81;while sqrt(R02+z2)-sqrt(R02*(1-1/mr0)2+z2)lambda/10 mr0=mr0+1;end ne0=mr0;rmax=15*f*phy;r=linspace(0,rmax,mr2);eta=linspace(0,2*pi,ne2);rho,theta=meshgrid(r,eta);x,y=pol2cart(theta,rho);r0=linspace(0,R0,mr0);eta0=linspace(0,2*pi*(ne0-1),ne0-1);rho0,theta0=meshgrid(r0,eta0);x0,y0=pol2cart(theta0,rho0);deta=R0/(mr0-1)*2*pi/(ne0-1);E2=zeros(size(x);E1=exp(j*k*zernikepol(3,mr0,ne0)*2.0e-4);for gk=1:ne2 for df=1:mr2 Rrho=sqrt(x(gk,df)-x0).2+(y(gk,df)-y0).2+z2);Rtheta=z./Rrho;opd=exp(j*k*(n-1)*(sqrt(RL2-rho0.2)-(RL-d)+d);Matlab 辅助激光光学分析与应用 23 Ep=-j/lambda/2*exp(Rrho*j*k).*(1+Rtheta)./Rrho*deta.*rho0.*opd.*E1(1:ne0-1,:);E2(gk,df)=sum(Ep(:);end end Ie=conj(E2).*E2;figure;surf(x,y,Ie);shading interp;axis(-rmax,rmax,-rmax,rmax)grid off;box on;toc;(a)像散光束聚焦 (b)彗差光束聚焦 图 1.17 畸变平面波的聚焦(入射光斑半径 1mm 平凸透镜凸面曲率半径 25mm，中心厚度 3mm)1.5.2 高斯光束通过透镜的聚焦 高斯光束与平面波不同的是，高斯光束是衍射自洽的，不会像平面波