微透镜的蒙特卡罗仿真研究

1、武汉大学本科毕业论文院系：信息工程学院专业：通信工程班级：07通信本作者：小马哥指导教师：宋子豪完成时间：2011年04月1附录15微透镜的蒙特卡罗仿真研究摘要基于光波的折反射定律和衍射理论，利用蒙特卡罗方法，建立微透镜的蒙特卡罗模型，并针对微透镜的传输特性问题，建立了高斯光束经过微透镜传输的数值模型。该模型的优点在于不需要使用处理球差的波相差或者相位变换概念，是直接在界面使用折射和反射定律，还可以根据两种媒质界面平面波的透射-反射菲涅耳公式考虑反射损失。该模型能同时考虑球差效应、衍射效应对微透镜成像的影响。利用MATLAB软件编写相关程序对该模型进行仿真研究，仿真结果显示：高斯光束通过微透镜

2、,束腰落点位置并不在几何焦点处，而是存在焦移；几何焦平面和束腰平面的光子分布呈现典型的干涉图样,而且几何焦平面上的聚焦光斑半径明显大于束腰平面上;关键词微透镜蒙特卡罗仿真数值模型MATLAB仿真目录摘要弓I言一、概述TOC o 1-5 h z1.1现阶段微透镜研究现状.1 HYPERLINK l bookmark10 1.2毕业设计课题内容1 HYPERLINK l bookmark12 二、Matlab软件和蒙特卡罗方法及基本理论方法1 HYPERLINK l bookmark14 Matlab软件2 HYPERLINK l bookmark16 2.2蒙特卡罗方法3 HYPERLINK l

3、 bookmark18 2.3基本理论方法5折射定律5反射定律52.3.3球差理论6 HYPERLINK l bookmark22 三、蒙特卡罗模型的建立和设计6 HYPERLINK l bookmark24 3.1入射光子的产生7 HYPERLINK l bookmark26 3.2衍射对光子落点坐标的影响7 HYPERLINK l bookmark28 3.3统计有效光子8四、数据结果分析84.1束腰位置的确定及焦移94.2束腰面与几何焦平面104.3束腰面与几何焦平面光子效率五、结束语12 HYPERLINK l bookmark34 参考文献.13英文题目、摘要、关键字14主程序15调

4、用程序21引言微透镜12是指几何尺寸在微米量级的透镜。通常微透镜的面型是球面和柱面，应用最多的微透镜是球面微透镜34。光器件的突破带动了光纤通信的飞速发展。光学微透镜广泛应用于光纤通信系统的各种器件，是一种非常重要的元件,所以对光学微透镜的研究将有助于更深入地了解各种微透镜的性能特点。由于微透镜具有许多优点，诸如衍射效率高导致的充分利用光能；用计算机可设计能产生任意波面的微透镜；薄片状重量轻；可复制价格低廉；可以微型化与阵列化，并且可以同微电子器件一起集成化，提供微光电一体的集成器件，在光通讯、光互连与光交换、光存储光学信息处理和微光学传感器等方面有着广泛的应用。蒙特卡罗方法的应用范围非常广泛

5、，包括：金融工程学，宏观经济学，生物医学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等。本论文所用到的正是蒙特卡罗的在粒子输运中的计算，在解决实际问题的时候应用蒙特卡罗方法主要有两部分工作：用蒙特卡罗方法模拟某一过程时，需要产生各种概率分布的随机变量。用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。该论文基于光波的折反射定律和衍射理论，利用蒙特卡罗方法，建立微透镜的蒙特卡罗模型，并针对微透镜的传输特性问题，建立了高斯光束经过微透镜传输的数值模型。该模型的优点在于不需要使用处理球差的波相差或者相位变换概念，是直接在界面使用折射和反射定律，并同时考虑球差效应、衍射效

6、应对微透镜成像的影响，得到的结果更加符合实际情况。因此，利用蒙特卡罗的方法对高斯光束通过微透镜成像进行仿真研究的结果对实际应用有着重大意义。 一、概述1.1现阶段微透镜研究现状微透镜和微透镜阵列被广泛的应用于各种光学耦合中，如光通信系统中的耦合器，发光器件，耦合装置67等。通常微透镜的面型是球面和柱面，应用最多的微透镜是球面微透镜由于微透镜不同于传统的透镜，它的结构尺寸达到微米量级，衍射效应就非常重要，同时球面结构带来的像差也不可忽略。现在对微透镜的研究有很多，有关球面和双曲面光纤微透镜的理论分析方法，很多都是将透镜视为理想薄透镜，即不考虑象差等对微透镜成像系统的影响。还有的研究方法则是基于几

7、何光学或基于波动光学，即只考率球差效应或衍射效应对光波的影响89。在微透镜定焦方面，利用数字图像处理技术的清晰度评价函数来进行定焦，但是这种方法却无法减小其定焦误差10。而本论文使用的蒙特卡罗方法，所获得的问题的结果实质上更接近物理实验结果，而不是经典的数值计算。所以，结论更加切合实际情况。1.2毕业设计课题内容利用光学知识，学习MATLAB软件。利用蒙特卡罗方法，建立了光纤基模光束通过微透镜成像的数值模型。利用此模型，对微透镜的传输特性进行了仿真研究。分析微透镜参数对传输特性的影响。二、Matlab软件和蒙特卡罗方法及基本理论方法2.1Matlab软件1120世纪70年代，美国新墨西哥大学计

8、算机科学系主任CleveMoler为了减轻学生编程的负担，用FORTRAN编写了最早的MATLAB。1984年由Little、Moler、SteveBangert合作成立了的MathWorks公司正式把MATLAB推向市场。到20世纪90年代，MATLAB已成为国际控制界的标准计算软件。Matlab在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。2.1.1Matlab软件的优势友好的工作平台和编程环境MATLA

9、B由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。简单易用的程序语言Matlab一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。这也是MATL

10、AB能够深入到科学研究及工程计算各个领域的重要原因。强大的科学计算机数据处理能力MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而前经过了各种优化和容错处理。MATLAB的这些函数集包括从最简单最基本的函数到诸如矩阵，特征向量、快速傅立叶变换的复杂函数。出色的图形处理功能MATLAB自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可用于科学计算和工程绘

11、图。新版本的MATLAB对整个图形处理功能作了很大的改进和完善，使它不仅在一般数据可视化软件都具有的功能（例如二维曲线和三维曲面的绘制和处理等）方面更加完善，而且对于一些其他软件所没有的功能（例如图形的光照处理、色度处理以及四维数据的表现等），MATLAB同样表现了出色的处理能力。应用广泛的模块集合工具箱MATLAB对许多专门的领域都开发了功能强大的模块集和工具箱。一般来说，它们都是由特定领域的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。目前，MATLAB已经把工具箱延伸到了科学研究和工程应用的诸多领域实用的程序接口和发布平台允许用户编写可以和MATLAB

12、进行交互的C或C+语言程序。另外，MATLAB网页服务程序还容许在Web应用中使用自己的MATLAB数学和图形程序。MATLAB的一个重要特色就是具有一套程序扩展系统和一组称之为工具箱的特殊应用子程序。应用软件开发(包括用户界面)在开发环境中，使用户更方便地控制多个文件和图形窗口；在编程方面支持了函数嵌套，有条件中断等；在图形化方面，有了更强大的图形标注和处理功能，包括对性对起连接注释等；在输入输出方面，可以直接向Excel和HDF5进行连接。2.2蒙特卡罗方法蒙特卡罗方法是以概率统计理论为主要理论，以随机变量抽样为手段的随机模拟方法。它的基本思想就是当所求解问题是某种随机事件出现的概率，或者

13、是某个随机变量的期望值时，通过某种“实验”的方法，把待求解的问题建立一个概率统计模型，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量(或随机向量)，就成为实现蒙特卡罗方法模拟实验的基本手段，它使得模型的随机变量的概率分布等于待求的解，通过对模型抽样获得所求随机变量的统计特征，从而得到所求问题的近似解。2.2.1分布函数和密度函数如果每次的试验结果可以用一个数&来表示，这个变量是随着试验的结果不同而变化的，称为随机变量。给定随机变量，其可能取值不超过实数x的概率P忆冬x)是x的函数，称为随机变量的累积分布函数，简称为分布函数，记作F(x)，即F(x)=P忆

14、x)-gx(1)分布函数指明了随机变量对所有可能取值的概率，时随机变量的全面描述。如果随机变量所取得值可以一一列举出来，而且&以各种确定的概率取这些不同的值，我们就称&为离散型随机变量。若随机变量&可取某一个区间a,b或(-g，g)中的一切值，而且其分布函数F(x)可以表示为xF(x)=jf(x)dx一g则称E为连续随机变量，称f(x)为的概率密度函数，简称为密度函数。根据密度分布函数利用随机抽样，得到满足该密度函数的解。2.2.2随机数最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布(或称矩形分布)。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单

15、子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。2.2.3蒙特卡罗方法特点蒙特卡罗

16、方法通过构造符合一定规则的随机数来解决数学上的各种问题。其特点是方法和程序结构简单、解题时受问题条件限制的影响较小。直接追踪粒子，物理思路清晰，易于理解。采用随机抽样的方法，较真切的模拟粒子输运的过程，反映了统计涨落的规律。不受系统多维、多因素等复杂性的限制，是解决复杂系统粒子输运问题的好方法。2.3基本理论方法2.3.1折射定律折射定律(lawofrefraction)或斯涅尔定律(SnellsLaw)由荷兰数学家斯涅尔发现，是在光的折射现象中，确定折射光线方向的定律。当光由第一媒质(折射率n1)射入第二媒质(折射率n2)时，在平滑界面上，部分光由第一媒质进入第二媒质后即发生折射。如图2-3

17、-1）实验指出：1）折射光线位于入射光线和界面法线所决定的2）平面内；（2）折射线和入射线分别在法线的两侧；入射角（3）入射角i的正弦和折射角i的正弦的比值，对折射率一定的两种媒质来说是一个常数.折射角图2-3-1射角；从浅显的说，就是光从光速大的介质进入光速小的介质中时，折射角小于光速小的介质进入光速大的介质中时，折射角大于入射角。此定律是几何光学的基本实验定律。它适用于均匀的各向同性的媒质。用来控制光路和用来成象的各种光学仪器，其光路结构原理主要是根据光的折射和反射定2.3.2反射定律光在光滑界面上反射时确定反射光线与入射光线传播方间关系的定律。几何光学的基本定律之一。如图2-3-2，入射

18、光线IO与界面在入射点O的法线ON所构成的平面称入射面，入射光线IO与反射光线律。图2-3-2OR的传播方向可分别用它们与法线ON的夹角0i和0r表示。通常把0i和0r分别称为入射角和反射角。通俗地说，就是反射光线与入射光线、法线在同一平面上；反射光线和入射光线分别位于法线的两侧；反射角等于入射角。可归纳为：“三线共面，两线分居，两角相等”。反射定律为：反射光线与入射光线同在入射面内。反射角等于入射角，即0i=0r。2.3.3球差理论1.球差亦称球面像差。光轴上物点发出的光束经光学系统后，不同孔径的光线与光轴夹不同角度的光线交光轴于不同位置（即像点位置不同），因此，在像面上形成一个圆形弥散斑，

19、这种现象称为球差。（如图2-3-3（a）一般是以实际光线在像方与光轴的交点相对于近轴光线与光轴交点（即高斯像点）的轴向距离来度量它。它可以表示成5l=lL（3）当边缘光线的交点位于近轴像点的左边时，球差为负；当边缘光线的交点位于近轴像点的右边时，球差为正。2.波像差基于波动光学理论，在近轴区内的一个物点发出的球面波经过光学系统后仍然是一球面波(惠更斯原理)，由于衍射现象的存在，一个物点的理想像是一个复杂的艾利斑。对于实际的光学系统，由于像差的存在，经光学系统形成的波面已不是球面，这种实际波面与理想波面的偏差称为波像差。如图2-3-3(b)所示：三、蒙特卡罗模型的建立和设计几何焦平面彳P束I臆

20、光线经过微透镜聚焦，其传输示意图如图3所示。AA平面为入射光的入射平面，透镜为半球透镜，其折射率为1.5。光轴为Z轴，F所在的面为几何焦平面；M所在的面为束腰平面。建立空间坐标系时，以O点为坐标原点，以光轴为Z轴，P所在直线为y轴，垂直纸面向外为x轴。根据高斯光束在微透镜中传输的物理模型，利用蒙特卡洛方法建立其数值模型，建模思路如下：、入射光线有指定的概率分布，通过这个分布，通过随机采样的方法产生入射光线，从而确定了入射光线的初始坐标；、通过光线追迹公式确定出射光线的几何传播方向和坐标，由于微透镜的衍射效应对出射光线的几何传播方向产生了影响，在模拟仿真时，这种影响可以看成是微扰则出射光线的轨迹

21、就是几何轨迹与衍射造成的微扰的叠加；、考虑到出射光线在观察面上的相干效应，因此在相干时间内落到观察平面的光子才被统计。3.1入射光子的产生以高斯光束通过微透镜聚焦为例，介绍使用蒙特卡罗方法，建立数值模型的思路和步骤：入射光束为高斯光束，高斯光束的束腰落在平凸微透镜的平面端面上，高斯光束入射光束的横截面上的场分布呈高斯分布，投射到微透镜平面(AA面)上的初始坐标为(P,0)，它的径向和角向分布分别为：f(P)=Cexp(-w24)5)f(0)=其中w是束腰宽度，C为常数。入射光通过微透镜时可以被认作是通过二个不同的系统，一个系统是球面像差效应系统，一个系统是衍射效应系统，聚焦结果可以被认作是二个

22、系统共同作用的结果。衍射对光子落点坐标的影响利用光线追迹公式可以求出在球差干扰下的光子落点坐标，下面考虑衍射的影响。当入射光束经过透镜时引起菲涅尔衍射，设入射场为U（x1,人），根据透镜位相变换公式，入射光束经过透镜后的光场分布为U(x,y)=11U(x,y)exp111(x2+y2)2f11P(x,y)11(6)式中：P(x1,y1)透镜的光瞳函数，其值在光瞳内为1，在光瞳外为0。根据菲涅尔衍射公式，距离为z的观察面上任意一点的光场分布为exp(jkz)j九zU(x,y)exp11jk(x一x1)2+(y一y1)22zdxdy11(7)根据光线追迹方程，确定离开出射光瞳的出射光线的方向和坐标

23、（X2,y2）。由于出射光瞳的窗口非常小，由于窗口衍射效应的存在，窗口对出射光线的方向产生影响，使得最终的出射光线方向发生改变，这个变化量就是微扰，这个微扰的的产生概率分布满足惠更斯-菲涅尔衍射公式。出射光线的轨迹就是几何轨迹与微扰项叠加。统计有效光子不同光子经过透镜到达观察面上时，不同光子之间存在光程差，在观察面上的某一观察点处，根据干涉理论，首先判断这些光子干涉相长还是干涉相消。在我们的模型中我们将观察面分成很多的面积微元，然后找到落在各个面积微元内的光子，再分别计算各个微元中光子之间的光程差，若该光程差为2兀的整数倍，则该点为干涉相长，若不为2兀的整数倍，则根据干涉公式计算出该点的有效的

24、光子数目。定义观察面上包含光子数目84%的区域为光斑半径，为得到精确结果，每一个模型都进行3次试验，对3次试验结果进行平均，得到最后的数据结果。、数据结果分析根据上述思路，我们编写了一套基于蒙特卡罗法MATLAB程序，用于模拟高斯光束通过微透镜后聚焦以及光线出射透镜后的传输过程。为了简单起见，仿真高斯光束通过平凸透镜进行聚焦。微透镜由于是半球透镜，因此用曲率半径（用R表示）来表示透镜厚度，入射的高斯光束束腰半径用w表示。4.1束腰位置的确定及焦移图4-1-1入射光束的焦移通过程序main1,main2（见附录）取多个曲率半径的多重计算，并对多组数据求平均值以提高数据的精确性，通过对数据的进一步

25、分析，表明在同时考虑球差、衍射、干涉的影响时束腰落点位置并不在几何焦点处，而是存在焦移的（如图4-1-1），而且随着曲率半径的不断增大，其焦移量也随之增大。而其本质原因是由于在保持束腰半径不变的情况下，随着曲率图4-1-2.束腰位置与焦距位置对比示意图从图4-1-2明显看出，焦距是曲率半径的2倍，而束腰位置与焦距位置并不在同一平面，而是存在焦移，并且随着曲率半径的增大，两者的距离也就越大，即焦移量随曲率半径的增大而变大。4.2束腰面与几何焦平面b)R=5几何焦平面(a)R=5聚焦光束束腰面图4-2-1.R=5微米束腰平面和几何焦平面光斑分布图在图4-2-1中，无论是几何焦平面还是束腰平面，光子

26、分布呈现典型的干涉图样。几何焦平面上的聚焦光斑半径明显大于束腰平面上(图(a)(b)，并且几何焦平面上光子分布较束腰面更加弥散不集中。这是由于在球差效应的影响下，当入射光线入射到透镜上的高度不同时,出射光线与光轴的交点位置也不同，导致其束腰位置未落在几何焦平面上，而是发生了偏移。(c)R=10聚焦光束束腰面(d)R=10几何焦平面(e)R=20聚焦光束束腰面(f)R=20几何焦平面图4-2-2.透镜曲率半径R分别为5微米,10微米米，20微米时，几何焦平面和束腰平面上的光斑分布图。从图4(a-c-e)，可以看出当入射光束的宽度一定时，透镜曲率半径越小，其光斑弥散程度增大。这是由于透镜曲率半径越

27、小，其衍射效应越大，则光子分布越分散；而透镜曲率半径越大时其衍射效应越低，光斑弥散程度越小，则光子分部越聚集。束腰面与几何焦平面光子效率五结束语本论文基于高斯光束在微透镜界面的折射与反射和波的衍射原理，综合考虑微透镜的几何形状、透镜的像差和衍射效应，利用蒙特卡罗方法，建立了高斯光束经过微透镜传输的数值模型，利用此模型，对高斯光束通过球型微透镜聚焦进行了仿真研究，仿真结果显示：由于球差的影响，使得聚焦光束束腰位置不在几何焦点处，存在焦移；相对于几何焦平面的光子分布曲线的最大值，束腰平面的光子分布曲线的最大值较大。微透镜不同于传统的透镜，由于它的结构尺寸非常小，达到微米量级，衍射效应就非常重要，同

28、时，球面结构带来的像差也不能忽略。而现有研究方法均没有综合考虑衍射效应和像差的对微透镜成像的影响。该论文模型基于几何光学和波动光学理论，利用蒙特卡罗方法仿真球面微透镜的光传输特性。该模型的优点在于不需要使用处理球差的波相差或者相位变换概念，是直接在界面使用折射和反射定律，并同时考虑球差效应、衍射效应对微透镜成像的影响，得到的结果更加符合实际情况。因此论文中的研究结果对微透镜的设计和应用具有现实指导意义。参考文献1贾正根；微透镜阵列J；微电子技术；1998年03期2SihaiChen,XinjianYi,HongchenWang,LuqinLiu,YingruiWang.MonolithicIn

29、tegrationofDiffractiveMicrolensArraysandInfraredFocalPlaneArraysJ,2002C.W.BarnardandJ.W.Y.Lit.Single-modefibermicrolenswithcontrollablespotsize.Appl.Optics1991(30):1958-1962C.W.BarnardandJ.W.YLit.Modetransformingpropertiesoftaperedsingle-modeefibermicrolensesJ.Appl.Optics1993(32):2090-2094尹增谦；管景峰；张晓

30、宏；曹春梅；蒙特卡罗方法及应用J；物理与工程；2002年03期C.A.Edwards,H.M.Presby,andC.Dragone.IdealmicrolensesforlasertofibercouplingJ.J.LightwaveTechnol1993(11):252-257H.M.Yang,S.Y.Huang,C.W.Leeetal.High-couplingtaperedhyperbolicfibermicrolensandtaperasymmetryeffectJ.J.LightwaveTechnol2004(22):1395-1401A.Yoshida.Shpericalabe

31、rrationinbeamopticalsystemsJAppl.Opitcs1982(21):1812-1816A.MiksandJ.Novak.PropagationofGaussianbeaminopticalsystemwithaberrationJ.Optik2003(114):437-440A.A.Alkelly.SpotsizeandradialintensitydistributionoffocusedGaussianbeamsinsphericalandnon-sphericalaberrationlensesJ.Opt.Commun2007(277):397-405焦勇；周

32、喻虹；基于MATLAB的快速图形化数据处理软件设计J.电子科技.2005(07)TheMicrolensMonteCarloSimulationStudyMonteCarlomethodisbasedonthediscountlawofreflectionanddiffractiontheoryoflightwaves,theMonteCarlomodelofthemicrolens,thesimulationofGaussianbeamthemicrolenstransmissionprocess.Andthetransmissioncharacteristicsofthemicro-len

33、s,numericalmodeloftheGaussianbeamtransferthroughthemicrolens.Theadvantagesofthismodelisthatdonotneedtousetreatmentdifferencebetweenthesphericalaberrationofthewaveorphasetransformationconceptisthedirectuseofrefractionandthelawofreflectionattheinterfacecanalsobebasedonthetwomediainterfaceplanewavetran

34、smission-reflectionFresnelformulatakesintoaccountthereflectionloss.Themodeltakingintoaccounttheeffectofsphericalaberration,diffractioneffectsonthemicro-lensimaging.PreparationoftherelevantproceduresofthemodelsimulationstudyusingMATLABsoftware,andsimulationresultsshowthat:KeywordsTheMicrolensMonteCar

35、loSimulationNumericalModelMATLABSimulation附录相关程序主程序mainl:包括参数设置、根据衍射场光强分布产生光子、衍射场计算clcclearticn1=1.5;%芯折射率即微透镜折射率n2=1.46;%包层折射率此主函数文件main1所有调用的文件为f33n3=1;%真空折射率D1=4.15;%单模光纤芯半径%D2=62.5;%包层半径%n4=2.148;n5=2.138;deff=0.9658;%n4=1.454;n5=1.445;deff=0.867;%bb=2.5;deltabb=1;%波导中心偏离，0:0.2:1向下偏离R=20;%曲率半径DD

36、=sin(asin(n3/n1)*R*2;%通光直径Lz=15;%z方向观察范围dz=0.05;%步长Q=44000;QQ=44000;%取样数N=15000;D=R*2;%D是透镜的高度(半球透镜)%D0=D1*(0.65+1.619/v4(3/2)+2.879/v46);f=R/(n1-1);%f为焦距hh1=R-sqrt(R42-(D/2).A2);%hh1透镜厚度此处为半球zt=-sqrt(R.A2-(D/2).A2);%zt是球心在z轴上的坐标，在笛卡尔坐标系下，坐标原点是光纤端面与光轴的交点。c=3e+8;%光速mu=pi*4e-7;%真空磁导率k0=2*pi/1.55;kn1=n

37、1*k0;%k0为波矢，kn1为在该介质中的波矢=2*pi*f=2*pi*wv=2*pi*4.15*sqrt(n1.A2-n2.A2)/1.55;%归一化频率w=4.0;%入射光束的束腰半径N1=0;N2=0;N3=0;%艾里斑分布rou3=linspace(0.0001,50,100000);cc=k0.*(DD/2).*rou3./f;I0=rou3*(k0.*(DD/2).A2./2./f).A2.*(2.*besselj(1,cc)./cc).A2;%艾里斑光强大小衍射符合一级贝塞尔函数圆孔衍射rou3*(k0.*(DD/2).A2./2./f).A2为中心点处的光强%I0=(k0.*

38、(DD/2).A2./2./f).A2.*(2.*besselj(1,cc)./cc).A2;C3=max(I0);pr3=unifrnd(0,50,QQ,1);%产生随机数(蒙特卡罗解题步骤)yy3=unifrnd(0,C3,QQ,1);N3=0;%以下for循环实现亮暗相间的光斑fori=1:QQN3=N3+1;r03=pr3(i);yy03=yy3(i);cc1=k0.*(DD/2).*r03./f;fr3=r03*(k0*(DD/2)42/2/f)42*(2*besselj(1,cc1)/cc1)42;%光强大小ifyy03fr3%舍选法r0airy(N3)=r03;elseN3=N3

39、-1;endendfori=1:N3theta00=2*pi*rand(1);%映射到x、y轴，将极坐标变成直角坐标xx0(i)=r0airy(i)*cos(theta00);%光子坐标yy0(i)=r0airy(i)*sin(theta00);end%以上程序根据艾里斑分布，产生在不同位置处的光子do=50;%观察半径sz=dz;%z方向观察步长sx=0.01;regz=Lz/sz;rouob=linspace(0,do,do/sx+1);fori=1:(regz+1)forl=1:(do/sx+1);%观察区域z=(1-i)*sz;r=(l-1)*sx;Y2(l,i)=f33(z,r,f,

40、DD,w);%调用子函数f33菲涅尔衍射积分公式Y2是积分号里面的Y5是积分号外面的Y5(l,i)=-j.*kO./(f+z).*exp(j.*kO.*(f+z)+j.*kO.*r.A2./(f+z)./2);endendY=Y2.*Y5;I=abs(Y).A2;rouob1=repmat(rouob,1,regz+1);I1=I.*rouob1;t=tocsaveR9_w5主程序main2:光源产生，包括通过衍射场分布，用舍选法得到衍射场分布、光线追踪、光斑处理，包括光斑大小、光子落点位置、焦移等clearclcticloadR9_w5Q1Q2=size(I1);pr5=unifrnd(0,

41、do,Q,1);yy5=unifrnd(0,max(max(I1),Q,1);N4=0;fori=1:(regz+1)N4=0;forl=1:QN4=N4+1;r05=pr5(l);yy05=yy5(l);ra=floor(r05/sx);IR=r05*(I(ra+1,i)*(ra+2)*sx-r05)+I(ra+2,i)*(r05-(ra+1)*sx)/(sx);ifyy05=(DD/2);%挑选光子j=j-1;continueendx00(j)=r00(j)*cos(theta0);y00(j)=r00(j)*sin(theta0);%光线追迹部分z0=0;x0=x00(j);y0=y00

42、(j);%calculatedtheangleofthewavevectorandx,yaxis波矢garma0=0;kz=k0*n1*cos(garma0);%garma是z轴夹角，alpha是与x轴夹角，beta是与y轴夹角。theta1=2*pi*rand(1);kalpa=pi/2-garma0;%angleofthewavevectoranditsprojectioninthexoyplane与z平面垂直。kxoy=kn1*cos(kalpa);%theprojectionofthewavevectorinthexoyplane？cos(pi/2)kxoy_x=kxoy*cos(th

43、eta1);%在xoy平面的横纵坐标值kxoy_y=kxoy*sin(theta1);kxx=sqrt(kz.A2+kxoy_y42);kyy=sqrt(kz.A2+kxoy_x42);calpha0(j)=(kxoy_x.A2+kn1.A2-kxx.A2)/2/kxoy_x/kn1;%angelbetweenthewavevectorandthexaxiscbeta0(j)=(kxoy_y.A2+kn1.A2-kyy.A2)/2/kxoy_y/kn1;%angelbetweenthewavevectorandtheyaxiscgarma0p(j)=cos(garma0);calpha=cal

44、pha0(j);cbeta=cbeta0(j);cgarmap=cgarma0p(j);x2p,y2p,z2p=F6(x0,y0,z0,calpha,cbeta,cgarmap,D,R);%求入射光线与透镜球面交点l1p=-x2p;l2p=-y2p;l3p=-z2p+zt;pusaip=F4(calpha,cbeta,cgarmap,l1p,l2p,l3p);%透镜的凸面上的入射角Pusaip=asin(sin(pusaip)*n1/n3);%透镜的凸面上的折射角ifpusaip=asin(n3/n1)%判断是否全反theconditionofthetotalreflectionj=j-1;c

45、ontinueendcalpha1p,cbeta1p,cgarma1p=F7(calpha,cbeta,cgarmap,l1p,l2p,l3p,Pusaip,hh1,x2p,y2p,z2p);%求折射光线的空间向量，里面有判定条件，表示的是取呈会聚趋势的折射光线x4p1(j)=x2p+calpha1p*(hh1+f-z2p)/cgarma1p;%像焦平面上的坐标y4p1(j)=y2p+cbeta1p*(hh1+f-z2p)/cgarma1p;calpha1(j)=calpha1p;cbeta1(j)=cbeta1p;cgarma1(j)=cgarma1p;x2(j)=x2p;y2(j)=y2p

46、;z2(j)=z2p;%微透镜凸面上各个光子的坐标endl1=sqrt(x00-x2).A2+(y00-y2).A2+(z2).A2);x4=x4p1;%像焦平面上的坐标y4=y4p1;%确定仅球差作用时候束腰NN=size(x2);L=Lz/dz;xx4=zeros(round(L),NN(2);yy4=zeros(round(L),NN(2);rr4=zeros(round(L),NN(2);b=zeros(round(L),NN(2);fori=1:L+1%任意观察面上的焦点xx4(i,:)=x2+calpha1.*(hh1+f-z2-(i-1)*dz)./cgarma1;%xx4,yy

47、4取第一行是只考虑球差时几何焦平面上光子的坐标，取第i行是只考虑球差时衍射焦平面上的光子坐标yy4(i,:)=y2+cbeta1.*(hh1+f-z2-(i-1)*dz)./cgarma1;rr4(i,:)=sqrt(xx4(i,:).A2+yy4(i,:)42);b(i,:)=sort(rr4(i,:);qiucha(i)=b(i,(round(NN(2)*0.84);endviewaxis=6;fori=1:L+1ifqiucha(i)=min(qiucha)shuyao=f-dz*(i-1);f-shuyao;shuyaobanjing=qiucha(i);x4=x4p1;y4=y4p1

48、;r4=sqrt(x4.A2+y4.A2);r41=sort(r4);i1i2=size(r41);i3=round(i2*0.84);qiucha1=r41(i3);endendr4=sqrt(x4.A2+y4.A2);r41=sort(r4);i1i2=size(r41);qiucha1=r41(round(i2*0.84);rantheta=2*pi*rand(1,j);x4p=zeros(1,j);y4p=zeros(1,j);%考虑球差、衍射同时作用fori=1:jr0airya(i)=r0airy(round(unifrnd(0.5,N3+0.5);endx44=x4+r0airy

49、a.*cos(rantheta);y44=y4+r0airya.*sin(rantheta);r44=sqrt(x44.A2+y44.A2);b=sort(r44);i3i4=size(b);both1=b(round(i2*0.84);both2=b(round(i2*0.8647);calpha11=(x44-x2)./sqrt(x44-x2).A2+(y44-y2).A2+(hh1+f-z2)42);cbeta11=(y44-y2)./sqrt(x44-x2).A2+(y44-y2).A2+(hh1+f-z2)42);cgarma11=(hh1+f-z2)./sqrt(x44-x2).A

50、2+(y44-y2).A2+(hh1+f-z2).A2);%菲涅尔衍射公式算的NN1=size(R0);fori=1:NN1(2)forl=1:jLEN1=find(R0(:,i)0);LEN11=size(LEN1);LEN2=LEN11(1);RO1(i,l)=R0(floor(unifrnd(1,LEN2+1),i);endendrant=2*pi*rand(NN1(2),j);X=xx4+RO1.*cos(rant);Y=yy4+R01.*sin(rant);%X,Y是球差加衍射时光子坐标分别取第一行和第M行是几何焦平面和衍射焦平面上的光子坐标RR=sqrt(X.A2+Y.A2);fo

51、ri=1:NN1(2)%光子出射透镜后所走路程l2(i,:)=sqrt(X(i,:)-x2).A2+(Y(i,:)-y2).A2+(hh1+f-z2-(i-1)*dz).A2);end%考虑干涉fori=1:NN1(2)X2=X(i,:);Y2=Y(i,:);l1_2=l1;l2_2=l2(i,:);j=1;j_1=1;jj=0;jjj=0;whilelength(X2)1len=(X2-X2(1).A2+(Y2-Y2(1).A2;a=find(len=0.6);iflength(a)=1aa1(j)=a+jj+jjj;j=j+1;jj=jj+1;X2(a)=;Y2(a)=;l1_2(a)=;l2_2(a)=;elsel3=l1_2(a)-l1_2(1)+l2_2(a)-l2_2(1);l4=a