ajk2002(数学建模)Word版

1、2002年全国大学生数学建模竞赛( 封面 )选择题号:普通组大专组ABCD( 根据所选竞赛题目在方框内打)学校名称 _××××××××_学生姓名 ××× 、 ××× 、 ××× 指导老师 ××× 全国大学生数学建模竞赛浙江赛区组委会二二年九月1 / 28摘要： 本文主要研究了车灯线光源长度在满足光照强度的设计要求下和功率节能的最优解策略。分别用正向光线追迹、逆向光线追迹、方程组模型求解。得到的结果

2、基本一致，但计算复杂度逐级下降、求解精度逐级上升，最后得出线性光源长度为4.060（mm）。由得出的解绘出测试屏上光强分布图十分附和车灯实际照射情况。最后，对设计规范从照射车距（安全性）、视认度（驾驶员）、车灯功率（节能原则）来评价其合理性。关键字： 光线追迹 离散化 等价光源车灯线光源的优化设计模型××× （计算机系2000） ××× （计算机系2000） ××× （计算机系2000）指导老师 ×××摘 要 本文主要研究了车灯线光源长度在满足光照强度的设计要求和功率节能的最

3、优解策略。分别用正向光线追迹、逆向光线追迹、方程组模型求解。得到的结果基本一致，但计算复杂度逐级下降、求解精度逐级上升，最后得出线性光源长度为4.060（mm）。由得出的解绘出测试屏上光强分布图十分附和车灯实际照射情况。最后，对设计规范从照射车距（安全性）、视认度（驾驶员）、车灯功率（节能原则）来评价其合理性。关键字 光线追迹 离散化 等价光源一、问题重述安装在汽车头部的车灯的形状为一旋转抛物面，车灯的对称轴水平地指向正前方, 其开口半径36毫米，深度21.6毫米。经过车灯的焦点，在与对称轴相垂直的水平方向，对称地放置一定长度的均匀分布的线光源。要求在某一设计规范标准下确定线光源的长度。该设计

4、规范在简化后可描述如下。在焦点F正前方25米处的A点放置一测试屏，屏与FA垂直，用以测试车灯的反射光。在屏上过A点引出一条与地面相平行的直线，在该直线A点的同侧取B点和C点，使AC=2AB=2.6米。要求C点的光强度不小于某一额定值（可取为1个单位），B点的光强度不小于该额定值的两倍（只须考虑一次反射）。二、模型假设假设 1：线光源是理想的线光源，分布均匀且无粗细。假设 2：只考虑经过一次反射。假设 3：所有的反射都符合反射定律。即，入射光线与反射光线的方向与反射面的法线矢量的夹角相等，并且所有的光线都在同一个平面内。并且属无损反射。假设4：电能全部转化成灯丝的光能。假设：不考虑配光镜对光线传

5、输方向上的影响三、问题分析由抛物面的已知数据：开口半径36毫米，深度21.6毫米。可以马上得出抛物面方程：z2y260x （由抛物线：y260x ，z0 ，和焦点坐标是（15,0），沿中轴线z轴旋转而成 ）用METLAB软件可以绘出抛物面的大致立体图形。如图：模型要求分析根据题意要求，在测试屏上要求C点的光强度不小于某一额定值（可取为1个单位），B点的光强度不小于该额定值的两倍。根据上图，分析抛物面对光反射的特性和由假设2可知，B点和C点的光强度取决于两种光的叠加经过反射镜反射后射入B（或C）的光线和线光源直接射入B（或C）的光线。问题 1，从线光源发出，经过反射镜后，照射在B,C上光强度。问

6、题 2，由线光源直接照射到测试屏上B,C点的光强度。汽车灯内部分析对车灯内部结构分析之后我们可以得到如下特性：1 电压U是恒定的。2 灯丝横截面面积为恒定极小S。（由假设1可得）3 在稳定发光时由于温度不变，则电阻率恒定。电阻R=l/S也恒定。因为车灯发光电路是一个非纯电阻电路。不考虑热效应，所以电源供给的能量完全转化为光幅射能。（由假设4可得）。由P=U2/R，我们可以推出P=U2/(l/S) 转换后：Pl=U2S/由上1，2，3，可得U2S/为常量C由此，我们可以得出线光源与功率成反比的关系Pl=C 根据式子要想功率尽可能的小，线光源要尽可能的大，且要满足在测试屏上B和C 点的光强比要求。

7、四、模型的建立和求解1，建立空间坐标系由题目中可知，因为这个状如抛物面的车灯的开口半径为36毫米,深度为21.6。则设车灯抛物面的对称轴为X 轴，车灯的抛物面的顶点作为我们的坐标原点。建立三维空间坐标系。根据已知条件，设在X-Y平面里的抛物线方程为y22px,且过点（21.6，36）（单位：毫米）则，可以求得p=30，那么，抛物线为y260x，z0。该抛物线沿中轴线x轴旋转，可以得出车灯的抛物面：z2y260x。空间图形如图：模型一：进行正向仿真设计光线追迹4主要思路：1计算反射光强度我们把线光源和抛物面都离散化成为离散点，针对离散后的线光源上的每一个点，与抛物面上的各个离散点，可以求出入射光

8、线。由于我们已知了抛物面的曲面方程，也就可以求得该抛物面上点的法向量。根据反射定律，由下面的公式（1），根据入射光线和法向量可以求得反射光线。公式（1）1推导和证明见附录1然后，由反射光线直线和25米处的测试屏相交，可求得交点坐标。用双重循环，遍历线光源上和抛物面上的所有点，记录测试屏上通过B,C点（在误差应许的范围内，也可以是，附近的点）处光线的数量N。计算直射光强度为使直射光线与反射光线含同样的光强量。（也就是跟反射的光线统一单位）。只有这样才能考虑叠加效果。只要把直射的光强也按反射光线同样的离散度转化为光线，这样就能统一单位了。而反射光线是根据抛物面的离散程度来转化的，所以直射光线也需要

9、经过抛物面来转化。转化方法如下：把测试屏上的B（或C）点与抛物面上的各个离散点直接相连。得到直线方程，与线光源线段去交，(方法同上)若有交点,则认为线光源上有一份光线对B(或C)点的光强有贡献。若无,检测抛物面上的下一点 ,重复上述步骤.直到遍历完抛物线上所有点。（记录个数M（直射）求出最终结果所得到的经过B点(或C)的光线数的和M+N，就可以表示线光源直射和反射到B(或C)点的光强。同样,求出C点光强,叠加后,比较B与C点的光强之比。重复上述步骤，线光源线段（15,y,0）范围线光源按步长1>0.001由大变小，逐步搜索，最终得出最优解为3.9280。计算机实现2：1 将抛物反射面离散

10、化，划分为m个面积相同的网格，取得每个网格中心点的空间坐标数据。2 线光源从一定步长d增长，没增加一个步长，将此时长度的线光源离散化为n个点。取的它的空间坐标数据。3 对线光源上的离散点和抛物面上的离散点进行两两连接，这样就取得此时长度下线光源发出的射向抛物面的所有入射光线。4 利用MATLAB上的函数surfnorm求得抛物面（曲面）上每个离散点的法线向量。5 根据入射光线和法线向量，也就是相当于已知入射光线向量和法线向量，通过公式（1），我们可以求得每条入射光线通过抛物面反射面反射后的反射光线。6 根据已知条件，测试屏上的两点B（25000，1300，0）和C（25000，2600，0），

11、所以我们在空间上取x25000这个平面做为测试屏。对5步鄹所求得的所有反射光线求取其在x25000平面上的交点，记录下所有点的坐标数据。7 以B点和C点为圆心，在平面x25000上做一微小半径r的圆，记录在此小圆内的通过6步鄹求得的交点的个数，分别做为此时长度线光源通过反射在B点和C点产生的光强度（即反射光强度）的衡量值Bf和Cf。8 将抛物面上的所有离散点和B点以及C点连线，分别计算连接B点的直线与线光源线段的交点个数和连接C点的直线与线光源线段的交点个数，分别做为此时长度线光源通过直射在B点和C点产生的光强度（直射光强度）的衡量值Bz和Cz。9 因为线光源长度是以一定的步长增加的，所以可以

12、取得满足条件(BfBz)>=2(Cf+Cz)(即B点光强不小于C点光强的两倍)的一系列长度值。取得长度最大的一个长度值做为线光源的最优长度值，此时所需的功率是满足条件(BfBz)>=2(Cf+Cz)下的最小功率。10 解过程中的数据图标如下：B点光强627627627627627C点光强292297297297297强度比2.147262.111112.111112.111112.11111线光源长度3.913.9123.9143.9163.9186276276326326322972972973023022.111112.111112.127952.092722.092723.9

13、23.9223.9243.9263.92811 根据以上数据和算法原理，我们得出线光源长度的最优解为3.9280mm。此时B点的光强度和C点的光强度的比值为2.09272。（源程序见附录2）模型二：逆向光线追迹法主要思路：根据光路可逆原理，从B，C点出发进行逆向光线追迹，即由反射光线去寻找入射光线。入射光线与线光源之间的公共点的存在，即表示所考虑的光线对B，C点的光强有贡献。再由B，C点光强数（光线数）和倍数关系，以及功率越小越好的约束条件，线光源逐步长搜索最优长度。 详细步骤：对于某一定长的线光源，因为线光源上的各个的点光源的光照强度一样，那么，在点，处的光照强度的比值就等于分别对于，点光源

14、有贡献的光线数的比值，分析光线射到B点的照射情况：考虑反射效果由B(或C)点与抛物面上的离散点（x0,y0,z0），那么就可以求出入射光线了。 由于抛物面的方程：F=z2y2-60x，那么，在M(x 0,y0, z0)处的法向量为：（FX，FY，FZ）=（-30，y0, z0）由计算原理中的公式（）可求出反射光线,该反射光线方程与线光源所在的直线去交。并且，y取值范围就是线光源的长度.若有交点,则认为光线对 B(或C)处的光强有贡献。遍历抛物面上的所有离散点，记录下这些对（或）有贡献的光线的个数。（记录个数N（反射）考虑直射效果为使直射光线与反射光线含同样的光强量。（也就是跟反射的光线统一单位

15、）。只有这样才能考虑叠加效果。只要把直射的光强也按反射光线同样的离散度转化为光线，这样就能统一单位了。而反射光线是根据抛物面的离散程度来转化的，所以直射光线也需要经过抛物面来转化。转化方法如下：把测试屏上的B（或C）点与抛物面上的各个离散点直接相连。得到直线方程，与线光源线段去交，(方法同上)若有交点,则认为线光源上有一份光线对B(或C)点的光强有贡献。若无,检测抛物面上的下一点 ,重复上述步骤.直到遍历完抛物线上所有点。（记录个数M（直射）求出最终结果所得到的经过B点(或C)的光线数的和M+N，就可以表示线光源直射和反射到B(或C)点的光强。同样,求出C点光强,叠加后,比较B与C点的光强之比

16、。然后，重复上述步骤，线光源线段（15,y,0）范围线光源长由大变小，逐步搜索，比较.最后可得出最优长度为.06 (mm)。模型二的计算机实现2：1 将反射面（抛物面）分解为721721个等面积的网格，得到721721个网格中心点的坐标数据集合X721,721，Y721,721，Z721,721。2 利用MATLAB中的函数surfnorm计算得到反射面（抛物面）上每个网格中心点对应的法线单位向量NX721,721,NY721,721,NZ721,721。3 在空间中取定测试屏上的两个测试点B点和C点（B(25000,1300,0),C(25000,2600,0)）。4 分别从B点和C点向反射

17、面（抛物面）上的每个网格中心点引向量。经过B点的向量数组为(Bx721,721,By721,721,Bz721,721)，经过C点的向量数组为(Cx721,721,Cy721,721,Cz721,721)。即从线光源发出的光线镜抛物反射面反射后经过B点与C点的反射光线向量B，C.这样就可以利用每个向量和对应网格中心点的空间坐标数据唯一的确定每条反射光线。5 利用公式求出每条经过B点和C点的反射光线的入射光线向量（Brx721,721,Bry721,721,Brz721,721）和(Crx721,721,Cry721,721,Crz721,721)这样我们同样可以利用每个入射光线向量和对应网格中

18、心电表数据确定每一条入射光线。6 分别对经过反射面反射到B点的入射光线和经过反射面反射到C点的入射光线求与线光源所在直线的交点个数，分别记为Mfb和Mfc，作为线光源在此长度下对B点和C点的反射光强度衡量值。7 分别求出经过B点和C点的反射光线与线光源所在直线的交点个数，并分别记为线光源在此长度下对B点和C点的直射光强度衡量值Mzb和Mzc。8 线光源长度以1>0.001的步长从0开始增长，用按步长搜索法求出长度最长的同时又满足(Mfb+Mzb)>=2(Mfc+Mzc)的线光源长度值（我们认为长度越长功率越小）。9 求解过程中的数据表格如下：半长度1.91.911.921.931.

19、941.951.961.97b点光线数309310310310311311311313c点光线数134137139140143143144146光强比2.305972.262772.230222.214292.174832.174832.159722.14384半长度1.981.9922.012.022.032.042.05b点光线数313313317317317317317317c点光线数146149154155158158159161光强比2.143842.100672.058442.045162.006332.006331.993711.96894半长度2.062.072.082.092

20、.12.112.122.13b点光线数317317317319321323323323c点光线数161162164166169173173174光强比1.968941.956791.932931.921691.899411.867051.867051.85632半长度2.142.152.162.172.182.192.2b点光线数323323325325325325325c点光线数176176179180181181181光强比1.835231.835231.815641.805561.795581.795581.7955810 根据以上数据和算法原理，我们得出线光源长度的最优解为4.0600

21、mm。此时B点的光强度和C点的光强度的比值为2.00633。（源程序见附录3）模型三:具体思路考虑反射的情况：反向光束法是用来求由光反射引起测试屏上的点B,C处的照度大小的情况（光线追踪与模型二的情况一样）.反向光束法的思想包括沿着真实光线相反的方向追踪一个想象的光束。射向照明系统的反向光束与光学系统的部件相交(抛物面反射镜)。这些交点可以等效地看成发光源，称为等价光源。它假想发射的光线到点B,C处照度的影响与点光源经过反射到点B,C处产生的照度影响等价(照度虽然与光源到点B,C的距离的平方成反比,但是由于车灯到测试屏的距离很大,在车灯里的点光源由于位置的不同所引起的与点B,C的距离不同的因素

22、就忽略不计)。考虑直射的情况。对于某一个线光源上的点光源,它直射引起的点光源个数为1(即线光源上的该点光源),。由于，照度是有同一个线光源产生的。又由于线光源上的每个点光源对，的照度的影响是一样的。则，点，处的照度也就与点光源(包括等价光源)的个数成正比.那么,点，的照度比就是可以反射光线到，的点光源（包括等价光源）的个数比。.具体的做法如下:先表示线光源在点B处时的照度情况1):在抛物面上设有一点M(x 0,y0, z0),那么到点B(25015,1300,0)(单位:毫米)的反射光线的向量为(x 0 -25015,y0 -1300, z0 ).由于抛物面的方程为：F=z2y2-60x那么在

23、M(x 0,y0, z0)处的法向量为：（FX，FY，FZ）=（-30，y0, z0）设入射向量为（Rx，Ry，Rz）那么由公式（） 求得：Rx =-1799 x 0 +1799*25015+60 y02-60*1300 y0+60z2 ()Ry=60 x 0y0 60*25015 y0 + y0-1300-2 y03 +2600 y02-2 y0 z02 ()Rz=60 x 0z0 60*25015 z0 +60 y02-60*1300 y0-2 z03+ z0 ()2):这条入射光线的方程如下： =由于这条入射光线是由线光源上的点光源产生的，所以这条入射光线与直线x=15z=0相交。把x=

24、15，z=0代入上面的式子，可以得出如下方程= ()3）又因为点M(x 0,y0, z0 )是在抛物面上，所以 y02+ z02 =60 x 0 ()以上由式子（）（）就可以得出y0f(x 0), z0 =g(x 0)4）有= ，那么，与线光源z=0,x=15的交点处的y 值可以用下表示 y=(0-z0)+y0 = q(x 0) ()从焦点开始，以一定的步长增加，对于每个y值yi，根据上面的六个方程，六个变量，就可以算出符合条件的所有点M(x 0,y0, z0)， 这个点的个数就是等价光源的点的个数。再加上一个直射的光源点个数1，设总个数为mi. 同理，可以求出反射光线经过点C时候的情况，对应

25、与每个与到点情况时同样的y值yi， 都有总的光源数（包括等价点光源）ni, 那么根据题目的要求，使得线光源在点B，C处的光照强度符合1：2的关系，只要逐步地以一定步长增加y的值，只要符合如下条件（取了y的个数为N）：（ ）/（250152+13002）2（）/（250152+26002）记录这个时候的所有符合条件的yN。在这些值里找yN最大的情况那么，yN也就是线光源的长度。此模型是前面两个模型的改进，理论上会有更加精确的解，由于时间有限，我们未来得及用程序实现。五、模型的评价和改进方法：模型一的评价： 这种方法思路清晰,模型简单且利用计算机编程容易实现,误差取决于抛物面和线光源的离散精度.步

26、长越小,精度越高,误差越小,误差分析简单. 但这种模型计算量随着精度的提高会成倍增加.在程序实现过程中,抛物面离散点数为m,线光源离散点数为n,则用循环计算入射光线、出射光线的时间复杂度就为0(m*n)。而要得到较高的精度，我们一般需要把m设为50000以上，n设为100以上，那么时间复杂度至少为5.0e+6，因此此中分析算法的缺点是计算机处理时速度慢，效率低，需要较长的时间来处理。 此模型不仅适用于反射面为抛物面,还适用于椭圆面,反射镜,复合反射面反射镜,甚至曲面反射镜)。对光源方面而言，此模型不仅适用于光源是线光源的光源设计，也适用于多线光源自由摆放的光源设计。所以此模型适用性较广，但在高

27、精度要求下，计算机运行速度必须相当快，才可能较快的得出理想的数据。所以此模型在推广过程中较依赖于运行计算机的硬件环境。模型二的评价：逆向光线追迹法只对抛物面进行了点的离散，然后直接让利用通过B、C点的反射光线求得的入射光线与线光源直线段去相交。这样做既避免了线光源的离散，提高了计算的精确度，又提高了计算机实现的程序的执行效率，因为循环次数的成倍减少大大缩短了程序的执行时间。另外，逆向光线追迹法只针对测试屏上的B点和C点进行考虑。题目的要求：衡量线光源长度的取值是否合理，满足条件之一就是B、C点上光强度的比值关系，这使我们这样的做法成为合理和可行。这样做，我们就避免了对一些无用点的计算，大大地减

28、少了所需统计和计算的数据，而又不降低反而是提高了精确度。我们尝试让用正向光线追迹的算法程序(1)和逆向光线追迹的算法(2)程序同时执行求解满足题意要求的线光源长度的任务，程序（1）执行了将近1个半小时，而程序（2）只执行了不到2分钟的时间。模型三的评价：这个模型中，使得点光源能够通过反射到点，处时，对于抛物面上的反射点M(x 0,y0, z0)是有条件限制的。这个限制可由式子（）（）得出y0f(x 0), z0 =g(x 0)，而模型二里只是离散地找了一些点，就会漏掉其它的同样符合条件的反射点。所以，模型三是比模型二更加地精确，误差越小。六、在有标尺的坐标系中画出测试屏上反射光的亮区的实现：1

29、、MATLAB软件实现：取10m×10m的测试屏，将测试屏划分为41×41个方格，取得每个方格中心点的坐标数据。对每个方格中心点利用逆向光线追迹的方法求得每个方格中心点的反射光强度（衡量值为光线条数），将光强度（光线条数）存贮为一个二维矩阵，这个二维矩阵就相当于一幅图象的颜色映射表，利用MATLAB里面的pcolor函数绘出该颜色映射表所对应的图象，此图象就是测试屏上反射光强度的分布图，颜色越浅表明光强度越大。（源程序见附录4）2、Borland Delphi实现在确定了线光源的长度以后，我们以如下的方法来寻找反射光的亮区:1、在线光源上取一段微小的长度，并将其看作点光源S

30、，其坐标为（Sx，Sy，Sz）；2、把反射抛物面看作数量极大的反射点组成的集合；3、按一定的顺序遍历所有的反射点，在遍历的过程中执行以下操作：A设每次得到的反射点坐标为（X0,Y0,Z0），由抛物面方程可以得到该点的切平面（也就是相当于反射时的镜面）的法向量为（X0/30，Y0/30，-1）B求出光源S关于切平面的对称点,也就是求出S的虚像。这可以通过先求出直线L的方程，再得到L与切平面焦点A，通过点的对称来得到S点的坐标。C从S点和反射点（X0,Y0,Z0），就可以确定反射光线所在直线M方程；D光屏所在平面方程为Z25015mm,将其代入直线M的方程就可以得到反射光与屏的交点。E大量交点的合

31、成就是所求的亮区。按上述的方法找到线光源上所有点对于反射抛物面上的所有点的反射光线与屏的交点。程序运行结果为下图。其中颜色越亮，表明该点的亮度值越大。其中横纵座标每个坐标单位的长度为726.2mm，光亮区域横坐标范围为（-3268.03268.0），纵坐标范围为(896.0896.0)。七、考虑规范的合理性当汽车在夜间行驶时，车灯必须有良好的人工照明，以使驾驶员有安全的视认条件。所谓安全的视认条件就是:对道路下目标的视距要大于车辆常速行驶的自动距离。 把测试屏放于车前25m处可以认为是常速行驶时的安全距离。 所以，测试屏放置的距离是合理的. 采用汽车前照灯照明时所得到的视认性3与前照灯照射所产

32、生的目标亮痕L0和背景亮度Lb对比度 C= （L0Lb）L0有关。此值越大，则视认越好。因此，增强前照灯的功率，使目标亮度L0增大，对提高视认性是有益的。通过线光源的不断按步长变化（也就是功率不断变化），绘测出一系列视认度C的点。用METLAB经过多项式拟合，绘出图形如下：B点图型： 关于取B点与C点的光强比要求为2比1的合理性的解释：根据模型，我们求出符合要求的线光源长度为4.06mm，此点对应的视认度大小为79.05。在图形中此视认度点的特征为，在顶点周围斜率最大的点。结合物理意义，他的合理性解释就是：首先，车前灯设计，要求考虑对于司机是视认度最佳的。（或非常接近最佳的视认度）其次，再考虑

33、车灯的功率要尽可能的小，以节省能源。 所以，设计时要取顶点周围斜率最大的点。也就是说，在保证视认度在一定高度的情况下，沿着功率减小的方向，取灵敏读最高的点。 因此，把这样的点作为设计规范点，是很合理的。C点图形 在不同功率下，对点B和点C的光线中含有光强量统一起来。（源程序见附录5）参考文献：【1】周军。FFR汽车前照灯的光学设计。照明工程学报，1999，10（2）：4649。【2】苏金明，阮沈勇。MATLAB6.1实用指南。电子科学出版社，2002。【3】周泰明，昭红，刘荣安。汽车前照灯的进展。复旦大学出版社，1989，6：3338。 【4】Han-Wen Tsai, Chien-Ping

34、Kung. Computer Assisted headlight Design and Research. SAE Technical Paper Series 950590. 1995.附录清单：附录1：公式（1）的推导与证明附录2：模型一中求解线光源最优长度值的源程序附录3：模型二中求解线光源最优长度值的源代码附录4：用MATLAB实现在有标尺的坐标系中画出测试屏上反射光的亮区的实现附录5：不同长度线光源在C点、B点的光强分布图的绘制附录1：公式（1）的推导与证明附录2：模型一中求解线光源最优长度值的源程序1、文件名：qiche.mY,Z=meshgrid(-36:0.05:36);X=

35、(Y.2+Z.2)/60;NX,NY,NZ=surfnorm(X,Y,Z);i=1;j=1;k=1for changdu=1:1:5 for dian1=1:1441 for dian2=1:1441 if X(dian1,dian2)<=21.6 rushefx(i,j,k)=X(dian1,dian2)-15; rushefy(i,j,k)=Y(dian1,dian2)-(changdu-3); rushefz(i,j,k)=Z(dian1,dian2); i=i+1 %else %rushefx(i,j,k)=nan; %rushefy(i,j,k)=nan; %rushefz(i

36、,j,k)=nan; %i=i+1 end end j=j+1,i=1 end k=k+1,j=1endPrushefx_i=rushefx./(sqrt(rushefx.2+rushefy.2+rushefz.2)+eps);Prushefy_i=rushefy./(sqrt(rushefx.2+rushefy.2+rushefz.2)+eps);Prushefz_i=rushefz./(sqrt(rushefx.2+rushefy.2+rushefz.2)+eps);2、文件名 dianc.mcdianc=1;for k=1:5 for j=1:721 for i=1:721 R=-2*N

37、X(i,j)2+1 -2*NX(i,j)*NY(i,j) -2*NX(i,j)*NZ(i,j);-2*NY(i,j)*NX(i,j) -2*NY(i,j)2+1 -2*NY(i,j)*NZ(i,j);-2*NZ(i,j)*NX(i,j) -2*NX(i,j)*NY(i,j) -2*NZ(i,j)2+1*Prushefx_i(i,j,k) Prushefy_i(i,j,k) Prushefz_i(i,j,k)' y=(25000-X(i,j)*R(2)/(R(1)+eps)+Y(i,j); z=(25000-X(i,j)*R(3)/(R(1)+eps)+Z(i,j); D=(y-2600

38、)2+z2 if D<=10000 cdianc=cdianc+1 end end endend附录3：模型二中求解线光源最优长度值的源代码文件名：Dc_L_JY.mY,Z=meshgrid(-36:0.5:36);X=(Y.2+Z.2)/60;NX,NY,NZ=surfnorm(X,Y,Z);i=1;j=1;for dian1=1:145 for dian2=1:145 if X(dian1,dian2)<=21.6 ffanshefx(i,j)=2.5e+4-X(dian1,dian2); ffanshefy(i,j)=2.6e+3-Y(dian1,dian2); bffans

39、hefy(i,j)=1.3e+3-Y(dian1,dian2); ffanshefz(i,j)=0-Z(dian1,dian2); i=i+1; %else %rushefx(i,j,k)=nan; %rushefy(i,j,k)=nan; %rushefz(i,j,k)=nan; %i=i+1 end end j=j+1;i=1; end cPfanshefx_i=ffanshefx./(sqrt(ffanshefx.2+ffanshefy.2+ffanshefz.2)+eps); cPfanshefy_i=ffanshefy./(sqrt(ffanshefx.2+ffanshefy.2+f

40、fanshefz.2)+eps); cPfanshefz_i=ffanshefz./(sqrt(ffanshefx.2+ffanshefy.2+ffanshefz.2)+eps); bPfanshefy_i=bffanshefy./(sqrt(ffanshefx.2+bffanshefy.2+ffanshefz.2)+eps); xiandian=1;bxiandian=1;bdian=1;cdian=1;dianb=1;dianc=1;yyy=0;p=1;dianb(p)=1;%b光线数dianc(p)=1;%c光线数while abs(2-dianb(p)/(dianc(p)+eps)&g

41、t;=0.0001&yyy<=10 bdian=0; cdian=0; xiandian=0; bxiandian=0;for j=1:145 for i=1:145 R=-2*NX(i,j)2+1 -2*NX(i,j)*NY(i,j) -2*NX(i,j)*NZ(i,j);-2*NY(i,j)*NX(i,j) -2*NY(i,j)2+1 -2*NY(i,j)*NZ(i,j);-2*NZ(i,j)*NX(i,j) -2*NX(i,j)*NY(i,j) -2*NZ(i,j)2+1*cPfanshefx_i(i,j) cPfanshefy_i(i,j) cPfanshefz_i(i,

42、j)' BR=-2*NX(i,j)2+1 -2*NX(i,j)*NY(i,j) -2*NX(i,j)*NZ(i,j);-2*NY(i,j)*NX(i,j) -2*NY(i,j)2+1 -2*NY(i,j)*NZ(i,j);-2*NZ(i,j)*NX(i,j) -2*NX(i,j)*NY(i,j) -2*NZ(i,j)2+1*cPfanshefx_i(i,j) bPfanshefy_i(i,j) cPfanshefz_i(i,j)' xiany1=(15-X(i,j)*R(2)/(R(1)+eps)+Y(i,j); xiany2=-Z(i,j)*R(2)/(R(3)+eps)+Y

43、(i,j); bxiany1=(15-X(i,j)*BR(2)/(BR(1)+eps)+Y(i,j); bxiany2=(-Z(i,j)*BR(2)/(BR(3)+eps)+Y(i,j); by=1300-24985*bPfanshefy_i(i,j)/(cPfanshefx_i(i,j)+eps); cy=2600-24985*cPfanshefy_i(i,j)/(cPfanshefx_i(i,j)+eps); if abs(xiany1-xiany2)<=1.0e-3&xiany1<=yyy&xiany1>=-yyy&xiany2<=yyy&

44、amp;xiany2>=-yyy xiandian=xiandian+1; end if abs(bxiany1-bxiany2)<=1.0e-3&bxiany1<=yyy&bxiany1>=-yyy&bxiany2<=yyy&bxiany2>=-yyy bxiandian=bxiandian+1; end if by<=yyy&by>=-yyy&Z(i,j)=0 bdian=bdian+1; end if by<=yyy&by>=-yyy&Z(i,j)=0 cdian=

45、cdian+1; endendendp=p+1;dianc(p)=cdian+xiandiandianb(p)=bdian+bxiandianfanc(p)=xiandianfanb(p)=bxiandianzhic(p)=cdianzhib(p)=bdianyyy=yyy+0.25end附录4：用MATLAB实现在有标尺的坐标系中画出测试屏上反射光的亮区的实现文件名：Gfb_S.m Y,Z=meshgrid(-36:1:36);X=(Y.2+Z.2)/60;NX,NY,NZ=surfnorm(X,Y,Z);cheshiy,cheshiz=meshgrid(-5000:250:5000);cheshix=cheshiy.*0+25000;li=1;lj=1;k=1;for wi=1:41 for wj=1:41 for dian1=1:73 for dian2=1:73 if X(dian1,dian2)<=21.6 ffanshefx(li,lj,k)=25000-X(dian1,dian2); ffanshefy(li,lj,k)=cheshiy(wi,wj)-Y(dian1,dian2); ffanshefz(li