Matlab辅助激光光学分析与应用（ 第三版 word文件）

二零一四年 第三版Matlab辅助激光光学分析与应用作 者：刘良清 单 位：武汉凌云光电科技毕业院校：华中科技大学激光技术与工程学院学 历：硕士研究生研究方向：自适应光学、非线性光学、激光光学、固体激光器件、激光工业应用2014年五月 第三版前 言（序）每个人都会有自己的理想和梦想，只不过我们常常走着走着就把它们弄丢了。也许你和我一样，小时候也曾信誓旦旦地说过自己的将来的理想是要当科学家什么的。那时候不是因为我们理想很远大，其实我们根本就不知道理想是什么。我们不知道真正科学家到底是什么，也不知道除了科学家之外还有什么可以叫做理想！当然，我现在没能成为科学家，将来也没有太大的可能会成为科学家。我既不是什么教授，也不是什么专家，顶多算是一个普通的技术员，普通得不能再普通了。记得研二的时候，导师还曾担心过我能不能正常毕业。整个大学过得似乎都有些浑浑噩噩，没什么明确的目标，只是盼望着大学毕业了能否找到一个温饱不愁的工作。现在我知道，这世界上根本就没有哪一样东西是真正会让你不愁的。很早，我就有想写一本书的想法。开始想写小说，结果构思来、构思去，发现好的情节似乎都被作家们写完了。突然有一个新奇的思路的时候，才发现自己文笔真的不怎地，开篇几段就写不下去了。罢了，作家之梦就这么完蛋了！研究生的最后一年，为了毕业论文和文章，慢慢自己学会了一点MATLAB编程，处理一些数据和图形。但是，后来发现之前做过的一些研究，突然有一天再去看的时候，竟然怎么也看不懂。于是，我就想，干脆把自己研究过的或者是纯属娱乐研究过的东西都总结在一起，把整个思路都写下来，这样可以方便自己以后参考。既然要写，那就得像模像样地写，有多少就写多少。我不知道什么时候会写完，照这个想法，也许一辈子也写不完。因为，总会有新的东西值得我们去研究。那么就这样吧，这就是我的第一本书了。书中肯定会有各种问题和错误，很多结论只代表我个人观点，如有雷同请找专家。尽信书不是好书！目 录MATLAB辅助激光光学分析与应用IMATLAB辅助激光光学分析与应用II作 者：刘良清II第一章 光的波动性和衍射11.1 Maxwell方程组和电磁波11.2 波动方程41.3 衍射51.3.1 小孔衍射61.3.2 双缝衍射141.4 波前畸变171.4.1 Zernike多项式171.4.2 畸变光束的衍射201.5 光束通过光学元件的变换241.5.1 平行光束通过透镜的聚焦24.2 高斯光束通过透镜的聚焦301.5.3 自聚焦透镜321.6 高斯光束361.6.1 高阶高斯光束361.6.2 高斯光束的传输变换411.6.3 高斯光束的干涉45第二章 激光谐振腔602.1 激光谐振腔的本征方程602.2 无源腔的Fox-Li迭代方法602.3 无源腔的矩阵特征向量方法662.3.1 平行平面腔672.3.2 双凹腔712.4 基模谐振腔782.4.1 高斯光束的ABCD定律与谐振腔802.4.2 基模谐振腔的稳定性912.4.3 基模谐振腔设计实例972.4.4 基模谐振腔设计的逆向方法107第三章 激光工程设计1183.1 速率方程的求解1183.1.1 主动调Q速率方程的求解1183.1.2 被动调Q速率方程的求解1253.1.3 脉冲泵浦速率方程的求解1323.2 非线性激光光学1363.2.1 非线性晶体1363.2.2 非线性晶体相位匹配161第四章 MATLAB数据与图像处理1724.1 Matlab图像处理应用1724.1.1 自制特殊图形172.2 图像数据处理1804.2 Matlab数据处理应用1974.2.1 批量数据处理1974.3 Matlab方程求解应用211含参数非线性方程组的求解211第一章 光的波动性和衍射1.1 Maxwell方程组和电磁波十八世纪中叶，James Maxwell将已知的各种电磁作用关系用一组方程组合起来，形成了一个方程组： (源于库伦定律的高斯定律)(1.1) (源于毕奥-萨瓦尔定律的高斯定律)(1.2) (法拉第定律)(1.3) (Maxwell修正的安培定律)(1.4)式中，和分别代表了电场和磁场分量。电荷密度描述路径空间单位体积内的电荷量分布；电流描述电荷的移动(单位电荷乘以速度)。表示真空介电常数，其值为。表示真空磁导率常数，其值为（或者）。在安培定律中引入了一个关键参数之后，Maxwell意识到，方程组构成了一个完美的电磁现象自洽理论。此外，方程组预言了电磁波的存在，并以光速传播。在Maxwell时代之前就已经有人对光速进行了测量，因此一个显而易见的结果(当时还难以令人置信)便是，光是一种高频振荡表现，类似并超越了支配电流和电荷的影响因素。而在此之前，光学还仍然作为一种独立于电学和磁学的主体进行讨论的。这里，我们不再对电磁学的基本知识进行详细的讨论，因为它们在普通物理课程中都有讲述，并且有大量的文献和书籍对其进行了细致的分析。但我们要简要的从波动方程出发，求解旁轴近似下的Maxwell方程组，得到激光传输与变换的基本方程，以方便我们后续的讨论和应用。为了体现Matlab在可视化方面的优势，我们先以一个简单的例子作为本书的开篇，以达到抛砖引玉的效果。在电动力学中，我们会遇到真空电磁场波动方程的旁轴近似解，众所周知，其解为具有高斯分布的电场复振幅：(1.5)式中，为光波传播常数。、是与光束有关的传播参数。分别表示为：(1.6)(0.7)(0.8)光束远场发散角为：(0.9)或者(0.10)我们可以用几行简单matlab程序就可以画出具有高斯分布的电场强度，如图1.1所示，图形美观，方便对光强的分布有一个感性的视觉认识。程序代码为：1- clear;2- clc;3- w0=0.5;4- r=linspace(0,3*w0,200);5- eta=linspace(0,2*pi,200);6- rho,theta=meshgrid(r,eta);7- x,y=pol2cart(theta,rho);8- Iopt=exp(-2*rho.2/w02);9- surf(x,y,Iopt);10- shading interp;11- xlabel(位置 /mm);12- ylabel(位置 /mm);13- zlabel(相对强度 /a.u.);14- title(高斯强度分布);15- axis(-3*w0,3*w0,-3*w0,3*w0,0,1);16- colorbar;17- colormap(hot);18- box on;19- grid off;图1.1 高斯光强分布另外，我们还可以画出高斯光束在自由传输过程中的强度变化，如图1.2所示，程序代码如下：1- clear;2- clc;3- lambda=1.064e-3;4- w0=0.5;5- ZR=pi*w02/lambda;6- z=linspace(-2*ZR,2*ZR,200);7- y=linspace(-4*w0,4*w0,200);8- py,pz=meshgrid(y,z);9- wz=w0*sqrt(1+(lambda*pz/pi/w02).2);10- Iopt=w02./wz.2.*exp(-2*py.2./wz.2);11- surf(pz,py,Iopt);12- shading interp;13- xlabel(位置 /mm);14- ylabel(位置 /mm);15- zlabel(相对强度 /a.u.);16- title(高斯强度分布的传输);17- colorbar;18- colormap(hot);19- box on;20- grid off;图1.2 高斯光束自由传输强度变化以上我们以简单的例子展示了Matlab在可视化方面的强大功能，但本文不再对Matlab的基本功能和语法常识进行介绍，我们认为本书的读者已经具备了基本的Matlab编程技巧。或者说，我们所做的只是将我们的实际运用跟读者进行交流讨论，促进大家共同进步。当然，我们会在一些比较关键的地方指出编程过程中需要注意的问题。1.2 波动方程当Maxwell统一了电磁理论以后，他马上意识到，波动可能是该方程组的解的形式。事实上，他希望找到一组满足波动形式的方程组，以辅助他完成找到真正的波动方程。既然已经知道了光是以波动方式传播的，基尔霍夫首先注意到了正好给出了精确的光速(之前就已经被测量过)，并且法拉第和克尔已经观测到强磁场和强电场会影响光在晶体中的传播。对初始接触Maxwell方程组的人来说，并不能一眼就看出它的解具有波动形式。但是经过适当的数学操作，我们就可以将它变为波动方程的形式。我们来推到电场E的波动方程，磁场B的波动方程的推到过程是类似的。我们将方程进行卷积，可得：(1.11)该方程可以由矢量微分恒等式简化：(1.12)卷积可由式代换，由此得到：(1.13)再由式代入上式，经过整理就可得到：(1.14)需要指出的是，上式中没有考虑到介质的极化。若考虑到介质的极化和实际一般光学问题中自由电荷为零的条件，上式修正为：(1.15)式中，为极化强度矢量。这样我们得到了一般的电场传播方程，该方程在非线性光学中有很重要的地位。当光在真空中传播时，式中的右边所有项均为零，方程简化为：(1.16)这样我们就得到了电场传播的波动方程形式。当然在有些实际问题中，式中右边的项并不是都为零，至少会有一项不为零，这与介质的性质有关。1.3 衍射考虑一个振动频率为的光场，其复振幅可以表述为，则它也必须满足波动方程：(1.17)由于(是假设)电场振幅的含时部分是显式给出的，则方程可以简化为：(1.18)式中是波矢量的大小。式就是所谓的赫姆霍兹方程。如果我们忽略波动的矢量特性，而只考虑它的振幅(这里不再详细讨论其过程)，那么在标量近似下，就得到了标量赫姆霍兹方程：(1.19)然后，我们考虑一束沿z轴传播的光束，它的电场复振幅写成的形式。我们将它代入标量赫姆霍兹方程式，得到：(1.20)在旁轴近似下，有。即是说，我们假设了电场的复振幅沿z轴传播方向是缓慢变化的，与平面波类似。但是我们允许振幅沿z轴在远大于波长量级的范围上有明显的变化。这样就得到了旁轴波方程：(1.21)求解方程式，得到：(1.22)于是电场的表达式为：(1.23)值得一提的是，基尔霍夫早在1887年就提出了著名的菲涅耳-基尔霍夫衍射公式：(1.24)式和在分母时具有一致性，并在指数上：(1.25)同时，式是式在满足条件下的菲尼尔旁轴近似。另外，如果进一步满足远场条件，就得到夫琅禾费衍射近似：(1.26) 小孔衍射假设光场透过一个圆柱对称的小孔，这时，孔径上的场分布可以写为：(1.27)这样，二维衍射积分可以简化为一维衍射积分。将式代入到菲尼尔衍射积分公式中，得到简化衍射积分式：(1.28)对角度的积分项，我们可以借助下面的公式完成：(1.29)式中，称为零阶Bessel函数。这样，式可以简化为：(1.30)式中的积分项也称为的汉克尔变换。在夫琅禾费衍射近似下，项等于1，积分项变为的汉克尔变换。于是夫琅禾费柱对称圆孔衍射方程为：(1.31)虽然经过了一系列简化，然而是复振幅通常都是不确定的，即使知道了强度分布，相位分布也可能是比较难预测的。当然，也可以通过辅助手段测量强度分布和相位分布。这里，我们以平面波入射为例，讨论圆孔夫琅禾费衍射问题。这时可用常数代替，不妨设为1。利用Bessel函数的递推关系，我们可以得到解析的圆孔夫琅禾费衍射公式：(1.32)图1.3 夫琅禾费圆孔衍射 (孔径，距离1m远)即使是解析式，我们还是不能直观地感受到衍射斑的样式。下面我们利用Matlab给出夫琅禾费圆孔衍射的强度分布。程序代码如下：1- R=0.1;2- lambda=1.064e-3;3- k=2*pi/lambda;4- z=1.0e3;5- r=linspace(0,2*1.22*lambda/2/R*z,201);6- eta=linspace(0,2*pi,201);7- rho,theta=meshgrid(r,eta);8- x,y=pol2cart(theta,rho);9- Bess=besselj(1,rho*R*k/z);10- Ie=4*pi2*R2*Bess.2./(rho*k).2/lambda2;11- surf(x,y,Ie);12- axis(-max(r),max(r),-max(r),max(r),0,max(Ie(:);13- shading interp;box on; grid off;14- figure;plot(x(1,:),Ie(1,:),k,x(101,:),Ie(101,:),k);程序中使用了这样一个参数，因为平面波假设具有最小的衍射角，这个参数就称为衍射极限角，所以在衍射区域里，只取二倍衍射极限范围就可以大致看出小孔的衍射特性。并且，对于解析表达，Matlab提供了Bessel函数工具箱，可以直接调用。当然，为了能够从图上看到跟实际观测相近的衍射环，可以将相对强度分布取四次开根号，如图1.3下图所示。图1.4 数值积分法计算的夫琅禾费圆孔衍射上面，我们使用的解析表达式绘图。接下来，我们采用数值积分方法直接求解式夫琅禾费圆孔衍射分布。这需要将目标平面划分网格，然后用网格格点上的复振幅代替领域内的平均振幅分布求相面上的振幅分布。为了简化计算量，我们不再计算二维分布，只计算一维分布，计算结果如图1.4所示。程序代码如下：1- R=0.1;2- lambda=1.064e-3;3- k=2*pi/lambda;4- z=1.0e3;5- r=linspace(0,2*1.22*lambda/2/R*z,201);6- eta=linspace(0,2*pi,201);7- rho,theta=meshgrid(r,eta);8- x,y=pol2cart(theta,rho);9- r0=linspace(0,R,201);10- eta0=linspace(0,2*pi,201);11- rho0,theta0=meshgrid(r0,eta0);12- x0,y0=pol2cart(theta0,rho0);13- deta=R/200*2*pi/200;14- E2=zeros(201,1);15- for gk=1:20116- for m=1:20017- for n=1:20118- E2(gk)=E2(gk)-j/lambda/z*exp(x(1,gk)2+y(1,gk)2)/z/2+z).19- *j*k)*exp(j*k*(x(1,gk)*x0(m,n)+y(1,gk)*y0(m,n)/z).20- *deta*rho0(m,n);21- end22- end23- end24- Ie=conj(E2).*E2;25- plot(r,Ie,k,-r,Ie,k);比较图1.3和图1.4，可以看出，数值积分法计算出的衍射斑与理论解析结果是一致的。另外，值得注意的是，从图上看到，一级衍射峰值比零级衍射峰值要低很多，其他级次就更低了。那为什么我们在实验中仍然可以清晰地看到许多级次的衍射环呢？那是因为人眼对光子的感光度比较高，而实验所用的激光都具有很高的光子数密度，因此还是可以看到很多的衍射级次。对于旁轴远场条件不够满足的情况，夫琅禾费近似不成立，即使是菲涅耳近似也可能不成立，这时必须使用更为严格的菲涅耳-基尔霍夫标量衍射公式：(1.33)式中，。我们仍然以平面波的圆孔衍射为例，对式数值计算不同相面上的衍射图样，如图1.5所示。程序代码如下：1- %exmp1_3_12- R=0.1;3- lambda=1.064e-3;4- k=2*pi/lambda;5- for z=0.5:5:25;6- r=linspace(0,2*1.22*lambda/2/R*z,201);7- eta=linspace(0,2*pi,201);8- rho,theta=meshgrid(r,eta);9- x,y=pol2cart(theta,rho);10- r0=linspace(0,R,201);11- eta0=linspace(0,2*pi,201);12- rho0,theta0=meshgrid(r0,eta0);13- x0,y0=pol2cart(theta0,rho0);14- deta=R/200*2*pi/200;15- E2=zeros(201,1);16- for gk=1:20117- for m=1:20018- for n=1:20119- Rrho=sqrt(x(1,gk)-x0(m,n)2+(y(1,gk)-y0(m,n)2+z2);20- Rtheta=z/Rrho;21- E2(gk)=E2(gk)-j/lambda/2*exp(Rrho*j*k).22- *(1+Rtheta)/Rrho*deta*rho0(m,n);23- end24- end25- end26- Ie=conj(E2).*E2;27- Ie=Ie/max(Ie);28- plot3(z*ones(size(r),r,Ie,k,z*ones(size(r),-r,Ie,k);29- hold on;30- end图1.5 不同传输距离处的圆孔衍射样式(孔径)可见，在远场条件不满的情况下，衍射样式有较大的差异。远场条件由菲涅耳数决定：(1.34)式中，a为对光束起实际限制作用元件的横向限度(半径)，L为传输距离。当F满足约为1的数量级时，菲涅耳衍射近似成立，当F满足远小于1的数量级时，夫琅禾费衍射近似成立。(a)(b)图1.6 矩孔衍射图样D= (a)z=20,(b)z=500当衍射孔为矩形孔时，同样利用式，可以数值求解不同相面上的衍射样式。如图1.6所示。其中左图为衍射样式的相对强度分布，右图为了突出与实际观测的近似，对强度分布取四次开根号。程序代码如下：1- %exmp1_3_2 矩孔衍射2- clear;3- tic;4- R=0.1;5- lambda=1.064e-3;6- k=2*pi/lambda;7- z=500;8- xmax=2*1.22*lambda/2/R*z;9- x=linspace(-xmax,xmax,61);10- y=x;11- x,y=meshgrid(x,y);12- x0=linspace(-R,R,61);13- y0=x0;14- x0,y0=meshgrid(x0,y0);15- deta=(2*R/60)2;16- E2=zeros(61,61);17- for bk=1:6118- for gk=1:6119- for m=1:6020- for n=1:6021- Rrho=sqrt(x(gk,bk)-x0(m,n)2+(y(gk,bk)-y0(m,n)2+z2);22- Rtheta=z/Rrho;23- E2(gk,bk)=E2(gk,bk)-j/lambda/2*exp(Rrho*j*k)*(1+Rtheta).24- /Rrho*deta;25- end26- end27- end28- end29- Ie=conj(E2).*E2;30- surf(x,y,Ie);31- shading interp;32- axis(-xmax,xmax,-xmax,xmax);33- figure;34- Ie=conj(E2).*E2;35- surf(x,y,Ie.(1/4);36- shading interp;37- axis(-xmax,xmax,-xmax,xmax);38- toc;另外，矩孔衍射也可以求得平面波衍射夫琅禾费近似解析表达式：(1.35)则远场衍射的强度分布为：(1.36)式中，a和b分别为矩孔的边长。它的衍射图样如图1.7所示，单缝衍射与此有相类似的分布，只不过去除一个自由度，这里就不再赘述。程序代码如下：1- %矩孔衍射的解析计算 exmp1_3_32- R=0.1;3- lambda=1.064e-3;4- k=2*pi/lambda;5- z=1.0e3;6- xmax=8*1.22*lambda/2/R*z;7- x=linspace(-xmax,xmax,200);8- y=x;9- x,y=meshgrid(x,y);10- 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;11- surf(x,y,IF.(1/2);12- colormap(hot)13- axis equal14- shading interp;图1.7 夫琅禾费矩孔衍射振幅绝对值图样 双缝衍射对于双缝夫琅禾费衍射，根据相干叠加原理，只需要将稍作调整并叠加：(1.37)式中，a为缝宽，d为缝的中心间距。计算结果如图1.8所示。程序代码如下：1- %双缝衍射解析计算exmp1_3_42- clear;3- R=0.1;4- d=4*R;5- lambda=1.064e-3;6- k=2*pi/lambda;7- z=1.0e3;8- xmax=3*1.22*lambda/2/R*z;9- x=linspace(-xmax,xmax,200);10- y=x;11- x,y=meshgrid(x,y);12- 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;13- 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;14- IF=IF1+IF2;15- surf(y,x,abs(IF);16- colormap(hot)17- shading interp;18- figure;19- plot(abs(IF(100,:).2,x(100,:);图1.8 夫琅禾费双缝衍射(孔径=，缝间距d=)如果不满足夫琅禾费远场条件或者入射光场不是平面波，则不能直接求取目标平面上的衍射场的解析表达式。这时，也只能使用菲涅耳-基尔霍夫标量衍射公式式。例如，假设入射到两狭缝的光场振幅均为超高斯分布：(1.38)当指数项中的n2时，称为超高斯分布，且n越大越接近平顶分布。同时，双缝衍射可以减少一个自由度，积分公式变为：(1.39)式中，。设缝宽a=，=，n=4,缝的间距d=，入射场的传播轴与缝的中心平行。计算过节如图1.9所示，程序代码如下：1- %exmp1_3_5 双缝衍射数值解2- clear;3- tic;4- R=0.1;5- d=0.2;6- ws=0.08;7- lambda=1.064e-3;8- k=2*pi/lambda;9- z=1.0e3;10- xmax=3*1.22*lambda/2/R*z;11- x=linspace(-xmax,xmax,61);12- y=x;13- x,y=meshgrid(x,y);14- x0=linspace(-R-d,R+d,61);15- y0=x0;16- x0,y0=meshgrid(x0,y0);17- E0=exp(-(x0-d).4/ws4)+exp(-(x0+d).4/ws4);18- deta=(2*R/60)2;19- E2=zeros(61,61);20- for bk=1:6121- for gk=1:6122- for m=1:6023- for n=1:6024- Rrho=sqrt(x(gk,bk)-x0(m,n)2+z2);25- Rtheta=z/Rrho; 26- E2(gk,bk)=E2(gk,bk)-j/lambda/2*E0(m,n)*exp(Rrho*j*k).27- *(1+Rtheta)/Rrho*deta;28- end29- end30- end31- end32- Ie=conj(E2).*E2;33- surf(x,y,Ie);34- shading interp;35- axis(-xmax,xmax,-xmax,xmax);36- figure;37- Ie=conj(E2).*E2;38- surf(x,y,Ie.(1/2);39- shading interp;40- axis(-xmax,xmax,-xmax,xmax);41- toc;图1.9 超高斯光束双缝衍射，左图为相对强度，右图为振幅绝对值(缝宽D=缝间距d=,距离z=1m)虽然可以减少一个自由度，但上面的程序仍然使用了二维数值积分，这是为了画出二维强度分布图。如果只是关注一维强度分布的话，则程序可以再简化，以减少计算量，节约程序运行时间。原则上，我们可以模拟任意复振幅分布入射的双缝衍射，甚至其他几何孔径的衍射，只要确定知道了复振幅分布的具体数值。例子很多，这里就不再一一演示了。1.4 （最新版2014年第三版 完整内容在淘宝发售）感兴趣的读者可在淘宝搜索书名找到。(未完待续)