大学毕设论文--文--车灯线光源的优化设计数学建模论文.doc

车灯线光源的优化设计数学建模论文
高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： B
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）： 重 庆 大 学
参赛队员 (打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人 (打印并签名)：

日期： 年 月 日

赛区评阅编号（由赛区组委会评阅前进行编号）：
高教社杯全国大学生数学建模竞赛

编 号 专 用 页

赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

评
阅
人

评
分

备
注

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： B
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）： 重 庆 大 学
参赛队员 (打印并签名) ：1.
2.
3.
指导教师或指导教师组负责人 (打印并签名)：

日期： 年 月 日

赛区评阅编号（由赛区组委会评阅前进行编号）：
高教社杯全国大学生数学建模竞赛

编 号 专 用 页

赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

评
阅
人

评
分

备
注

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

车灯线光源的优化设计
摘要
车灯作为汽车的一个重要部件，不仅影响车子造型的美观，而且是夜间行车时必要的照明工具。如何设计车灯线光源的长度，使它既能达到一定的照明效果又能充分地节约能源，是一个值得探讨的问题。
本文在满足给定设计规范的条件下，对车灯线光源的优化问题进行了分析研究。首先根据几何光学的反射定律，求出了线光源上任意点经过抛物面任意点反射所得光线的方程。其次利用所得直线方程，根据能量守恒推导出点光源的一条入射光的强度与测试屏上相应点的光强度的关系。对于固定的点光源和测试屏上的一点，利用反射方程可以确定有限个反射点，对这些点进行叠加，求出点光源经抛物面一次反射后在屏幕上某点所得的光强度。将线光源离散化成若干个点光源，用点光源的性质计算了线光源在屏幕上给定点B、C的光强度。以设计规范对B、C点光强度的要求为约束条件，线光源的最小功率为目标，利用数值计算方法计算出功率最小的线光源长度4.15毫米。用反射光线方程与测试屏方程，将反射光线在屏幕上的“照亮点”用点光源作标记反射点坐标表示出来（见式子（6））。用一定步长取遍线光源及反射点坐标，逐个计算相应的照亮点，就能得到线光源在测试屏上的反射光的亮区（见图3）。
在考虑反射光传播过程时，本文采用了曲面变换的雅克比行列式，描述了光束截面扩大的效应，简洁的得到了由于面镜对改变光的传播方向所引起的光强度的变化，另外，针对反射光的亮区图案作了定性分析。最后，在文末讨论了该规范设计的合理性，就B、C点位置特征的合理性进行分析，并对于二灯制前照灯设计规范惊醒了改进，以及近远光的考虑。
一问题重述

安装汽车头部的车灯的形状唯一旋转抛物面，车灯的对称轴水平的指向正前方。一种典型车灯的开口半径为69毫米，深度为21.6毫米。经过车灯的焦点，在与对称轴相垂直的水平方向面对成的放置一根长度均匀分布的线光源（中心位于焦点），其每一个发光都相同。要求在某一设计规范标准下确定光源的长度，使其功率最小。
该规范如下：在焦点F正前方25米处的Ａ点竖直放置一测试屏，屏与连线ＦＡ垂直，用以测试车灯的反射光，在屏上过Ａ点引出一条水平的直线，在该直线上Ａ点的同侧取Ｂ点和Ｃ点，使得ＡＣ＝２．６米，ＡＢ＝１．６米。要求Ｃ点的光强度度不小于某一额定值（可取为１个单位），Ｂ点的光强度不小于该额定值的两倍（只须考虑一次反射）。
　　　为了优化车灯光源的设计，使之既满足以上的设计规范又能够达到节能耐用的标准，现在要求如下问题：
（１）在满足给定设计规范的条件下，计算线光源的长度，使线光源的功率最小。
（２）对得到的线光源长度，在有标尺的坐标系中画出测试屏上的反射光的亮区，得到一个直观的灯光照明效果的图像。
（３）讨论给定的设计规范的合理性，看它能否比较好的反映出车灯照明效果的好坏。
二、问题分析
以车灯抛物面定点为原点，抛物面的对称轴为坐标系的ｚ轴，车灯所对正前方为ｚ轴的正方向，建立空间直角坐标系。则根据抛物面开口半径为３６，深度为２１．６，可以求得车灯抛物面的方程如下：
ｚ＝（）
其焦点坐标为Ｇ（０，０，１５），用Ｍａｔｌａｂ绘制出该抛物面的图形如下：

由给定的设计规范，整个“车灯－测试屏”的检测系统简化结构如下：
假设线光源的亮度是均匀分布的，可以认为它是由无数发光强度相同的点光源构成。设线光源长度为Ｌ，总发光功率为Ｐ，在线光源上任取一小段长为dl的线元，其功率为dp，则有：
dp＝dl

将ｄｌ长的线光源看成一个功率为ｄｐ的点光源，求出每个点光源ｄｌ对Ｂ、Ｃ两点光亮度的贡献，然后把所有点光源的光强度贡献积分，就可以得出Ｂ 、Ｃ两点的总光强度。
三、模型假设
１．线光源亮度均匀，它的中心位于旋转抛物面的焦点，水平放置在焦平面上；
２．不考虑灯抛物面对光线的多次反射，只考虑一次反射。　
３．设计规范中要求的Ｃ点的最小光强度为１；

四、模型的建立与求解
１．反射光线与测试屏交点坐标的确定
车灯抛物面焦点的坐标Ｇ（０，０，１５），线光源所在直线ν过G点且平行于ｘ轴，可以求出其直线方程为；

测试屏所在平面β与z轴垂直，求得其方程为：
Z=25015
取B、C两点在A点的X轴正方向一侧，于是得到：
B点坐标：（1300，0，25015） C点坐标：（2600，0，25015）
直线v与抛物面的交点为：（30，0，15）和（-30，0，15），所以线光源的最大长度为60毫米
设点光源上一点T（3t，0，15）发出的光线射到抛物面上某一点P0（x0，y0，z0），则入射光线Tp0的方程为
抛物面上点P0处的法线方程为
有P0
在抛物面上，故有：
60Z=x+y (1)
计T点关于法线的对称点为T（x1，y1，z1），则由线段TT 的中点 在法线上可得的:
（2）
又因直线TT与法线垂直，有（x-t）x+(y-0)y-30(z-15)=0
联立（1）（2）（3）可求出x,y,z的表达式如下：
x=x+[(-900tx-ty)/(900+x+y)]
y=y+[2txy/(900+x+y)]
z=(405000-1800tx+1350x+x+1350y+2xy+y
反射光线PT的方程为在此处键入公式。
pT与测试屏平面β的交点记为Q（x，y，z），其坐标满足下列方程组：

解之得：

把（1），（4）代入（5）得到下式：

2。关于反射点的个数的分析
命题：任何一个确定的点光源T和测试屏上的一点Q，在抛物面上只有有限个反射点使从点光源T发射出的光线经过反射后可以到达D点。
证明：（6）式完整描述了屏上任意一点Q、线光源上一点T与抛物面上对应的反射点的关系。由表达式可知x，y均只与x0，y0，t有关，可以记为；

当给定Q、T时，x，y，t均为常数。于是以上函数关系就转化为下列方程组；

这就是求两条平面曲线的交点。由于这两个曲线方程的形式不同，不可能重合。因此它们的交点最多只有有限个。
三。测试屏上任意一点光强度的确定（6）式确定了从线光源上一点T发出的光线经过抛物面上一点P0反射后在测试屏β任意一点Q*（x，y，z）的坐标值代入，就可以得到所有通过Q*点的反射光线所满足的关于x0，y0，t的关系式。对给定的t求出所有过Q*的反射光线在Q*产生的光强度，再对t的线光源L上积分就可以求出Q*的总反射光强度E0，它是一个关于线光源长度L和线光源总功率P的函数。
点的光强度包含两个部分，除总反射光强度外，还有总直射光强度：
=+
其中直射光强度E1为线光源不经反射直接射到Q*的光强度
4．光强度的物理分析
设车灯光源的总发光功率为Ｐ，光通量为F，则有如下规律：
　　F=ＣＶＰ……………………………………………　（７）
其中：　　Ｃ为比例常数：
　　　　　Ｖ为人眼对于波长光线（白光）的视觉灵敏度。
有发光强度Ｉ的关系式：　
　　　　　　　　　　　　Ｉ＝
其中：为立体角，对于电光源，＝４∏
由于点光源发出的光为球面波，球面上任意一点的光强度都相同，有：Ｉ＝　……………………………………………　（８）
对于照度（光强度）Ｅ有：
　Ｅ＝……………………………………………　（９）

其中：ｒ为发光点到照射面的距离：
　　　为照射面发现方向与照明方向的夹角。
联立（７）（８）（９）解得：
Ｅ＝　……………………………………………　（１０）
由上式可知：光强度（照度）与线光源总功率Ｐ成正比，与投射角的余弦成正比，与距离ｒ的平方成反比。
５．直射光线所产生的光强度分析
由于线光源与测试屏的距离为２５米，而线光源的最大长度为６厘米，可把线光源当成在焦点处的点光源，其产生的光强度：
＝　　…………………………………………　（１１）
其中：＝＋（ｄ为测试屏上Ｑ点到A电的距离） ， F=ＣＶＰ
6.B、两点光强度与线光源的功率长度的关系
把点取为Ｂ点，就可以求出Ｂ点的总反射光强度，它是关于Ｐ、Ｌ的函数，设为：
＝（Ｐ，Ｌ）
又射线光源在B点的总直射光强度为，它也是关于P，L的函数，设为：
＝（Ｐ，Ｌ）
则B点总光强度有：
+=（Ｐ，Ｌ）+（Ｐ，Ｌ）
设计规范条件要求B点的光强度不小于2，即：

=（Ｐ，Ｌ）+（Ｐ，Ｌ）2
同理，对才C点也可以得到：
=（Ｐ，Ｌ）+（Ｐ，Ｌ）1
这样（12）（13）分别确定了一个在以Ｐ，Ｌ为坐标系的平面内地两个区域，两个区域的重叠区即代表了多有满足设计规范的Ｐ，Ｌ的范围，园区该重叠区的Ｐ值最小的店所得到的对应的L值就是所要求的线光源的长度。

7.反射光线与入射光线能量关系的确定
光线从点光源T出发到抛物面上点的传播距离为r，有：
＝＋+
到达处的光强度为：
＝　
式中L为线光源的总长度，F为线光源发出的总的光通量。设一小束光线到达时垂直于T截面面积为d，照射到抛物面上的面积为　d，d在xy平面上的投射面积为d，该光束反射到测试屏上的Q（x，y，z）点处的面积为d，则由几何关系可知：
=
=

第九页空缺

因为由 难以写出x0,y0关于x,y的显式表达式，但是通过固定屏上点Q（x,y,25015）和点光源T，然后搜索抛物面上的反射点P0，结合前面的证明知P0只有有限个点，可以叠加出T点经过抛物面反射到O点的光强。
以下为竖直求解的算法（相应程序见附录1）
考虑T点附近△t的一段线光源可近似为点光源，对Q点，代入（6）中可数值求解出几组（x0,y0）解。2）对解的判别：用数学软件求出所有的实根，首先判断P0是否超出反射面的范围： （因为抛物面开口半径为36毫米）：其次判断反射光线是否交与抛物面开口圆面的内部。
3)由可能的解P0算出T点对Q点光强度的贡献 。因为 的表达式只与 ，P，L 有关，如下式：
=(+)

其中R=为焦点G与Q的距离，为直射光强度：而=，为反射光能损失系数，
为发射光的方向向量，其实只要讨论这一相对量。再对t从-L到L，以为步长，进行累加==q
其中为；q
为Q点处的总光强度。
4）每小段线光源△t的发光功率P△t/L=p△t为恒定的值，对固定A，B，当|t|=0，△t，2△t……n△t，分别求出△，i=0，1，2，……….n，n≤30/△t
对B点：q（B）=∑△（B）
对C点：q（C）=∑△（C）
又因为B要大于额定值的两倍，所以取min（q（B）/2，q（c））而线光源的总发光功率为P=（2n+1）△t þ
所以满足相对于额定值的单位光强度所需的发光功率为：
=
从1开始，改变n的值，算出发光功率Pn。取△t取最小的功率对应的n值，则长度L=（2n+1）△t为最优线光源设计长度。
5）首先以1mm为步长， =10000000则算得 为下表：
n
1
2
3
4
5
Pn
1.0345
0.0025
0.0027
0.0032
0.0038
由表容易看出Pn的最小值对应的n=2，则知最优长度在4毫米附近，再以0.05毫米为步长，同理算得下表：
n
32
33
34
35
36
37
38
39
40
41
42
Pn
2.36
01
1.71
54

1.49
27
1.38
20
1.31
86
1.28
00
1.25
60
1.24
16
1.23
03
1.22
72
1.22
7
易发现在=41时，最小，所以最优长度为：
L=0.05（2n+1）=4.15（毫米）
五.测试屏上反射光亮区的绘制
要求出测试上反射光线的亮区，就是要求所有反射光线与测试屏交点的分布范围。对于每一个点光源在测试屏上都有它们自己对应的亮区，把所有点光源对应的亮区重叠起来形成的总的包络面就是线光源的反射光的亮区。
我们考虑用计算机编程来绘制出亮区。
[算法设计的分析]
由于给定点光源的位置（即值）和抛物面上的任意一个反射点的坐标Xo,Yo（Zo可以有抛物面方程确定，所以不是自由变量），就可以唯一确定一条反射光线，从而可求出该反射光线与测试屏的交点坐标，即这条光线在测试屏上所“打出”的一个亮点的位置可以确定。
再考虑这条反射光线是否只经过了一次反射。一般的做法是用该反射方程与抛物方程联立，看是否有除反射点（Xo,Yo,Zo）以外的交点。但是在这里我们从另一个角度去考虑，即：如果该反射光线所在直线从抛物面开口处的圆中穿过，则它不会经历多于一次的反射。
这样，求出反射光线在z取值为21.6时点的另外两个坐标x,y的值，只要它们满足以下不等式就不会经过多于一次的反射：
+
如此，求解条件和约束条件都确定了，接下来就可以设计程序算法了。由于给定一个反射点（Xo,Yo,Zo）只要给定其Xo,Yo的值就可以了。把抛物面投影到xoy平面，其投影为一圆，方程为：
+
把该圆按照半径r和中心角θ划分成网格或者按照直角坐标x,y划分网格。这样给定一个t值，然后对于每一个网格上的点，确定一条相应的反射光线，求出它与测试屏的交点，即在屏上打一个点；变化t值，重复上述过程。这样最后就得出整个线光源发出的光线反射后在测试屏上对应的照亮区域（用标点出来的）。
[算法实现后的结果]
用以上算法在Matlab 中变成求解（程序见附录2），以直角坐标划分网格，当t在（-2.075，2.075）范围内变动（即线光源长度L为4.15）时，测试屏上亮点的分布如下图（图中标尺的单位为毫米）：

第十三，十四页空缺

发生很大的畸变，而距离焦点越远处的点的亮区的畸变速度越慢；
2．t值为正的点光源在屏上的照亮区域明显地偏向x轴的负方向，即：车灯线光源左边部分的反射光主要照亮汽车右前方，而右边部分的反射光主要照亮汽车左前方，存在光路的交叉。

六．设计规范的合理性分析：
安装在汽车头部的车灯（前照灯）要满足行车前方道路有明亮而均匀的照度，又不致使对方来车驾驶员眩目或有不舒适感，需对其进行精确的调整（对其车灯的光束位置和灯丝的长度设置等）来满足使用需求，模型结合对汽车前照灯配光性能的具体要求，用测试屏对前照灯照射位置和光强度进行检测，来判断此前照灯的好坏。检测车灯时，如果该设计规范是非常合理的，则仅从B、C的光强度，就能反映出整个测试屏的光强度来。
测试屏离车灯25米，屏上过被测试两点B、C形成的线与线光源平行。测量了点B、C的光强度。从测试屏的光强度分布图（图3）来看，车子的前照灯射出的光线在测试屏呈往左右延伸的光带，且往两侧逐渐减弱，让屏上的线BC水平放置，则该设计重点测试了行车前方左右两侧的光强度分布；由于反射镜是个旋转抛物面，且光源呈水平放置的线状，反射光不会过大的扩散，因此在垂直方向没有设置测试点；B、C点在直线上A点的同一侧，也就知道了另一侧相应对称点的光强度，因此B、C在同侧会得到更多的有关光强度分布的消息；B=1.3米，C＝2.6米，其长度与实际汽车在路上奔跑所辖的路宽差不多，其距离也是相当合理。归纳之，这样的测试方法是非常合理的。
从不同的几个测试屏上光强度分布图可发现，B、C分别处在两个不同光强度区域，它们之间有条较明显的明暗截止线。
25 米是汽车前照灯传统测试方法的常用长度。在正常驾车情况之下，与前车保持两秒钟车速的距离是安全的。根据科学研究，一般驾车员的制动反应时间在0.39～0.63秒之间，采取制动的时间在1～1.35秒之间，两者之和为1.39～1.98秒，小于两秒。而通常前后车的车速基本相同，所以在制动时间相同的情况下制动距离也基本相同。两秒钟的间隔，就可以保证前后车辆的有效制动。再考虑交通规则所限制的车速，25米是一个安全合理的车距。而司机必须能看到前方这么远的事物，所以测试面上的点的光强度不能少于某一值。25米是一个比较合理的测试距离。

七．设计规范的进一步分析
1．二灯制前照灯
实际运行的汽车是有两个车灯的，他们的光束互相影响。为了便于分析，不妨令两个车灯的距离为1.3米，且两个灯的类型相同。如下图所示：
C
B

在二灯制前照灯下测试屏上的光强度其实就是每个车灯在测试屏上的光强度叠加。先遮住射出的光线，则只在单灯F的情况下考虑点A、B、C的光强度，同理也得到、、在灯的照射下的光强度，并且C与是对称的。显然，只要测到点B、C的光强度，也就得到整个测试屏的光强度。

2．车灯的远近光分析
每个车灯都会发出远光或近光。前照灯的配光应使其近光具有足够的照明和不眩目，远光具有良好的照明。
二灯光测试：
左、右灯投射的光束互相影响。汽车左、右前照灯之间的间距，根据车型的不同，大约在0.9—1.5m之间。众所周知，汽车前照灯光束，不是平行光，而是以一定角度发散投射出去的。根据实际观察的结果，在同时开亮左右灯的情况下，在距离前照灯3m远处，左灯近光光束右侧向上倾斜的亮区，与右灯近光光束左侧水平线上的暗区产生重叠现象。这种重叠现象在前照灯发光强度比较好的情况下更为明显。这种重叠现象，使得右灯的近光光束的暗区不暗哦，给右灯近光光束的明暗截止线水平部分的准确判别以及明暗截止线拐点的准确判别都带来不可避免的不良影响。而在距离前照灯1m处，则不会出现这种重叠现象。所以在测试近光时，测试屏应在距离前照灯1m处。

八．参考文献
[1] 荆工 史尔 应用光学 国防工业出版社 1973
[2] 吕林根 许子道 解析几何（第三版） 高等教育出版社 1987
[3] 网址 ：1568
[4] 云舟工作室 Matlab数学建模基础教程 人民邮电出版社 2001

附录1
p=10000000;
m1=0;
q=[0,0];
c=0;flag=0;e=1;
[x2,y2]=solve(‘x+(900*t+t*x^2-t*y^2)/(900+x^2)/((405000-1800*t*x+1350*x^2+x^4+350*y^2+2*x^2*y^2+y^4)/(30*(900+x^2+y^2)/60)*(25015-(x^2+y^2)/60)-1300,y+2*t*x*y/(900+x^2+y^2)/(900+x^2+y^2)/((405000-1800*t*x+1350*x^2+x^4+1350*y^2+2*x^2*y^2+y^4)/(30*(900+x^2+y^2))-(x^2+y^2)/60)*(25015-(x^2+y^2)/60)’);
q2=0;
n=length(x2);
for i=1:n
x(i)=x2(i);y(i)=y2(i);
end
for i=1:n
if imag(x(i))==0&imag(y(i))==0&(x(i)^2+y(i)^2)<1296&
f1(x(i),y(i),t,21.6)^2+f2(x(i),y(i),t,21.6)^2<1296
a=[x(i)-t,y(i),(x(i)^2+y(i)^2)/60-15];flag=flag+1;
b=[x(i),y(i),-30];r=norm(a,2);
a=1/norm(a,2)*a;b=1/norm(b,2)*b;
a=1/jcb(x(i),y(i),t);
q2=q2+0.85*abs(g)*p*abs(a*b’)/(4*pi*r^2*abs(b*[0,0.1]’));
else
end
end
m1=p/(4*pi*(1300^2+25000^22))+q2;
end
p=10000000;
c=0;m2=0;e=1;
[x2,y2]=solve(‘x+(9000*t+t*x^2-t*y^2)/(900+x^2+y^2)/((405000-1800*t*x+1350*x^2+x^4+1350*y^2+2*x^2*y^2+y^4)/(30*(900+x^2+y^2))-(x^2+y^2)/60)*(25015-(x^2+y^2)/60)-2600,y+2*t*x*y/(900+x^2+y^2)/((405000-1800*t*x+1350*x^2+1350*y^2+2*x^2*y^2+y^v4)/(30*(900+x^2+y^2))-(x^2+y^2)/60)*(25015-(x^2+y^2)/60)’);
q2=0;
n=length(x2);
for i=1:n
x(i)=x2(i);y(i)=y2(i);
end
for i=1:n
if imag(x(i))==0&imag(y(i))==0&(x(i)^2+y(i)^2)<1296&
f1(x(i),y(i),t,21.6)^2+f2(x(i),y(i),t,21.6)^2<1296
a=[x(i)-t,y(i),(x(i)^2+y(i)^2)/60-15];flag=flag+1;
b=[x(i),y(i),-30];r=norm(b,2)*b;
a=1/norm(a,2)*a;b=1/norm(b,2)*b;
g=1/jcb(x(i),y(i),t);
q2=q2+0.85*abs(g)*p*abs(a*b’)/(4*pi*r^2*abs(b*[0,0.1]’));
else
end
end
m2=p/(4*pi*(2600^2+25000^2))+q2;
end
t
m1
m2

附录2
T=-2.075:0.1:2.075;n=length(T);
r=-36:1:36;m=length(r);
s=-36:1:36;l=length(s);z=25015;h=21.6;
for i=1:n
t=T(i);
for j=1:m
for k=1:l;
if r(j)^2+r(k)^2<36^2
flag=sqrt(f1(r(j),s(k),t,h)^2+f2(r(j),s(k),t,h)^2);
if flag<=36
ox=f1(r(j),s(k),t,z);
oy=f2(r(j),s(k),t,z);
plot(ox,oy);
hold on
end
end
end
end
end