光学原理测量硬盘振动.docx

光学原理测量硬盘振动几何光学测量硬盘振动模型 摘 要 本文通过理论分析和与实验数据的对比, 研究了通过几何光学测量硬盘振动问题, 并求出了硬盘盘面的局部表面方程，分析了不同分辨率对局部表面方程的影响，较好的解释了硬盘振动的规律。 首先研究不考虑盘片中心的上下振动的情形. 针对问题做出了几何模型, 利用所给数据,并参照查阅资料，最终求出考虑盘片中心上下振动情况下的结果，并对求解结果采用代数检验符合要求, 提出模型的改进方向。 对于问题一, 查阅了硬盘的工作原理和内部结构，从几何上进行理论分析，讨论了一条和多条光线对模型测量的影响，通过分析光路、成本、技术等因素得出结论采用两束光线的方案为最优。 对于问题二, 首先建立了不考虑盘片中心的上下振动的理想模型，通过几何关系和数学定理,顺利求出了盘片局部表面方程的参数a、b，并得到c=0；其次在前面计算基础上建立了考虑盘片中心的上下振动的理想模型，通过角度关系有效地解出了c的值。 对于问题三，通过分析得知，屏幕分辨率对于测量结果的影响其实也就是对于a、b参数值的影响。利用MATLAB生成几组反射面倾斜角度的随机值，可以求出一系列B1和B2在接受屏上的坐标值，再用fix函数对该值进行某一精度下的四舍五入取整，这样获得的盘片局部表面方程就包含了接受屏分辨率的影响，再通过最小二乘法求得实际值与改变值差的平方最小的值，其对应的分辨率为最佳分辨率。 在结果的分析与检验中, 改变数据与实际数据的对比图, 直观地反映出分辨率对于测量结果的影响，验证了模型的正确性与方法的可靠性，但两者之间仍存在偏差, 对此我们作出了合理的解释。 在模型的改进中, 着重分析了文章前一部分未考虑在内的各种不定因素对硬盘振动的影响, 并提出了两种改进方案：第一种利用盘片水平瞬间求出入射光线的方向向量；第二种是找到接受屏上的读数点坐标与入射光线的对应关系，便可通过代值法求的参数值；第三种是问题三的求解中应考虑分辨率对c的值的影响及n值的有限性。 关键词: 硬盘振动 几何光学 局部表面方程 分辨率 最小二乘 1问题的重述 硬盘是计算机上的重要部件，正向着更小、更快、密度更高的方向发展。如何减小盘片的振动，成为关键问题，于是怎样检测盘片的振动更加重要。由于硬盘转速很高，不适宜接触式测量，容易想到光学测量。 1 有一种新技术测量物体表面获得振动情况，是利用类似镜面反射的几何光学原理工作的：选用同一点光源发出两束光线，照射有反射能力的被测物体平面，产生的两条反射光线再照射到水平放置的接受屏上。接受屏能够检测出光线照射点在接受屏平面上的坐标位置，问题是怎样利用图2中Oxy平面的两个坐标值反算出被测物体当前的平面方程，就可利用平面方程的变化得到振动情况。 图1是用两条光线检测的原理图，点光源O为坐标原点，E1和E2是盘片瞬间位置的反射点，此时盘片局部表面方程为z=ax+by+(c+D)，其中a,b,c为未知量，由于硬盘在振动，平面就偏离了初始位置z=D(D<0)；接受屏方程为z=h，其中h>0为接受屏的高度,B1和B2是其上的照射点，这样利用它们确定a、b和c获得盘片局部表面方程。 问题 1：请解释这种测量技术采用两束光线的原因。 问题 2：请建数学模型分析盘片局部表面方程的测量原理。 问题 3：试分析接受屏分辨率对盘片局部表面方程的影响。 提示：可利用计算机模拟的办法求出一系列B1和B2在接受屏上的坐标值，再对该值进行某一精度下的四舍五入取整，这样获得的盘片局部表面方程就包含了接受屏分辨率的影响。 2问题的分析 21 问题一 根据硬盘的工作原理和内部结构可以知道，硬盘工作时盘片时刻在振动，不同时刻对应有不同的方程，若采用一条光线则无法求解，同时从成本、技术出发，两条光线为最佳方案。 22 问题二 z 点光源 y O x 接受屏 B1 B2 z=hhy , D B1 B2 O 被测物体平面 x , E1 E2 图2 接受屏示意图 二光束法测量硬盘振动时，两条入射光线与其反射光线所形成的平面均垂直于反射面，又经过共同的点光源O，且过O点有且只有一条直线垂直于反射面，故可以求出反射面的法线。根据镜面反射原理，入射角等于反射角，同时已知OE1、法线、B1E1的方向向量，可以求出E1点的坐标，利用点法式可得出平面方程，即a、b、c的值。 图1 原理图 2 23问题三 通过分析得知，屏幕分辨率对于测量结果的影响其实也就是对于a、b参数值的影响。利用MATLAB生成几组入射角的随机值，可以求出一系列B1和B2在接受屏上的坐标值，再用fix函数对该值进行某一精度下的四舍五入取整，这样获得的盘片局部表面方程就包含了接受屏分辨率的影响，再通过最小二乘法求得实际值与改变值差的平方最小的值，其对应的分辨率为最佳分辨率。 3模型的假设与符号说明 31 模型的假设 (1)两条入射光线的方向向量为已知，且固定不变; (2)接受屏上的坐标原点O'与O点所形成的直线与Z轴重合； (3)硬盘盘片振动过程中不会发生变形，局部表面始终为平面； (4)光线发生反射时为全反射，且传播过程中没有干涉、衍射等现象； (5)硬盘平放，盘片振动连续、稳定。 32 符号的说明 (m1,n1,p1) 表示直线OE1的方向向量; (m2,n2,p2) 表示直线OE2的方向向量; (x1,y1,h) 表示B1点的坐标; (x2,y2,h) 表示B2点的坐标; (x11,y11,z11) 表示E1点的坐标; (x21,y21,z21) 表示E2点的坐标; L1 表示面OE1B1的法线; L2 表示面OE2B2的法线; A 表示反射面的法线的方向向量。 4模型的建立与求解 41求解局部平面方程 根据几何定理得知，若平面A与平面B均垂直于平面C，假设AB相交于直线L，则L必垂直于平面C。又知入射光线与反射光线所形成的平面与反射面垂直，故题中的入3 射问题可简化为图3，如图。 图3 光线入射 根据数学定理，面OE1B1和面OE2B2的法线L1、L2分别为： L1=(m1,n1,p1)´(x1,y1,h) L2=(m2,n2,p2)´(x2,y2,h) 由于两个法线相交所形成的平面与反射面平行，故可以求得反射面的法线的方向向量A： A=L1´L2=(m1,n1,p1)´(x1,y1,h)´(m2,n2,p2)´(x2,y2,h) 因为反射面的方程为z=ax+by+(c+D)，因此其法向量为(a,b,-1)，又两平行向量差乘的绝对值等于0，即： A´(a,b,-1)=0 由上便可求出未知参数a、b， a=(p1´x1-m1´h)´(m2´y2-n2´x2)-(m1´y1-n1´x1)´(p2´x2-m2´h)(-(n1´h-p1´y1)´(p2´x2-m2´h)+(p1´x1-m1´h)´(n2´h-p2´y2)(m1´y1-n1´x1)´(n2´h-p2´y2)-(n1´h-p1´y1)´(m2´y2-n2´x2)(-(n1´h-p1´y1)´(p2´x2-m2´h)+(p1´x1-m1´h)´(n2´h-p2´y2)b=通过上面的求解过程可以得知用一束光无法得出局部平面方程，故题中选用二束光，下面求c的值。 41.1不考虑盘片中心的上下振动 硬盘振动过程中盘片中心始终固定在主轴上，即振动平面均通过(0,0,D)点，带入反4 射面方程即可求得c的值，结果为： c=0 41.2考虑盘片中心的上下振动 在图4中， OE1为入射光线，O'E1为法线，E1B1为反射光线，根据反射定理，入射角等于反射角，故ÐOE1O'=ÐO'E1B1,由题OE1、O'E1、E1B1的方向向量分别为：(m1,n1,p1)、(a,b,-1)、(x1-x11,y1-y11,h-z11) ，根据空间两直线的夹角公式得： m1a+n1b-p1m12+n12+p12a2+b2=+a(x1-x11)+b(y1-y11)-(h-z11)a2+1b2+1(x1-x11)+(y1-y11)+(h-z11)222图 4 等角示意图 此外，E1在直线OE1上，同时也在反射平面上，即： x11m1=y11n1=z11p1ax11+by11-z11+c+D=0 由方程可求的c的值，由于c值较长，详情见附录1。 42接受屏分辨率对局部表面方程的影响 42.1问题的分析 盘片振动过程中以偏离原平面的角度为考察对象，对于盘片表面上任一点，该点的z值反映了该点偏离原平面的情况。而z=ax+by+(c+D)，所以接受屏分辨率对局部表面方程的影响即是考察不同分辨率下z值的变化量的情况。即对其表面方程系数a，b的影响。 盘片时刻在振动，设其偏离初始状态z=D的角度为q，当q取不同值时，得到不同的B1，B2，对B1，B2取不同的精度10-4，10-5，10-6，10，10，并在这些精度下四舍五入取整，最后利用MATLAB软件求解，通过a，b前后变化来分析接受屏分辨率对盘片局部表面方程的影响。 42.2问题的求解 -7-85 图5盘面振动解析图 依据正弦定理，可得到接受屏上B1点在O-XYZ坐标上xB1的计算公式，如下： OCsin(p=OBsin(2OD+q-a)=p2OB-q)sin(2a-2q)sin(p2-a+2q)B1F=DFcot(a-2q)据题意知道B1点坐标为 x1=y1=B1F+OD 其中a为入射光线在平面OXZ或平面OYZ上的投影与OZ轴的夹角，当a在不同平面上取不同值时，就可求得相应的在接受屏上的反射光线点的位置。例如可分别取入射光线OE1在平面OXZ或平面OYZ上的投影与OZ轴的夹角a=p4,a=p3，同时令则求得接受屏上B1点在O-XYZ内的坐标x1,y1，同理可得B2点的坐D=-0.05,h=0.01，标。 然后，对B1，B2点的坐标取不同的精度10-4，10-5，10-6，10-7，10-8，并在这些精度下四舍五入取整，得到五组不同分辨率下的B1，B2点的坐标根据问题二可求的相应的五组a、b值，利用最小二乘法从中得到最优的分辨率。 42.3算法的实现 由资料得知，普通盘片的振动半径为44.5mm，振动范围小于50纳米，故得到q的取值范围为ç0,èæ5044.45´10-6ö÷，用øMATLAB在取值范围内随机抽取100组值,见表一。 6 表一 100组q的随机值 0.0891´100.0260´100.0692´100.0065´100.1069´100.0213´100.0384´100.0960´100.0924´100.0967´100.0320´100.0513´100.0426´100.0743´100.0566´100.0471´100.0156´100.0639´100.0683´100.0943´10-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5p0.0943´100.0840´100.0397´100.0017´100.0500´100.0668´100.0781´100.0343´100.1031´100.0326´100.0305´100.1052´100.0615´100.0790´100.0170´100.0228´100.0547´100.0830´100.0915´100.0924´10-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5 0.0767´100.0559´100.0022´100.0218´100.1005´100.0601´100.0985´100.0228´100.0679´100.0500´100.0021´100.0224´100.0385´100.0195´100.0785´100.0456´100.0011´100.0609´100.1037´100.0341´10-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5 0.0894´100.0348´100.1012´100.0818´100.0224´100.0699´100.0073´100.0306´100.0591´100.0284´100.1003´100.0952´100.0198´100.0336´100.0920´100.0524´100.0427´100.1048´100.0766´100.0501´10-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5 0.1006´100.0013´100.0154´100.1076´100.0756´100.0829´100.1112´100.0461´100.0482´100.0528´100.0857´100.0798´100.0744´100.1102´100.0417´100.0936´100.0725´100.0990´100.0224´100.0588´10-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5-5当a=4时，对每一个q值根据公式即可求得B1点坐标(x1,y1)，其中根据公式x1=y1。 同理当a=p3时，可求得B2点的坐标(x2,y2),且x2=y2。 根据B1，B2的坐标，由公式可以求得a、b的准确值，再把B1，B2取不同的精度，分别为：10-4，10-5，10-6，10，10，并在这些精度下四舍五入取整，-7-8分别求得不同精度下ai,bi。最后利用MATLAB软件，求得相应的五组a、b值，将得到的五组a、b的值与未经四舍五入的a、b的值做差的平方运算，å(a-ai)2+(b-bi)2，再取平均数，比较其大小即可得出结果，选出最优的分辨率，见表二。 表二 最小二乘法结果 n = 4 n = 5 n = 6 n = 7 n = 8 7 å(a-a)i2+(b-bi)2100 0.4211´10-43 0.2548´10-44 0.1825´10-46 0.6026´10-48 0.5928´10-50 5模型结果的分析 å(a-a)以n为横坐标，i2+(b-bi)2100为纵坐标作图如下，当n=8时所对应的分辨率对结果的影响最小，并且可以作出推测：随着n的不断增大，即分辨率越高，所求出的结果就越准确。 图 6 图 7 8 6模型评价与改进 61 模型的评价 6.1.1 模型的优点 (1) 问题2构建模型过程中，充分考虑了盘片中心是否上下振动两种情况，使讨论更加全面，结论更加可靠。 (2)问题3构建模型时，查阅了大量资料，得出了准确有效的已知条件，为后面的计算提供了保障。 (3) 采用图示、表格的方法直观地反映了所建模型求解数据与实测数据的对应关系; (4) 模型具有全面性，容易推广。 6.1.2 模型的缺点 (1) 该模型数据处理比较复杂, 运算较为繁琐; (2) 由于数据的不充足, 对某些参数的确定具有主观成分; (3) 没有考虑干扰因素的影响, a、b、c的值仍然存在误差。 62 模型的改进方向 6.2.1改进方向一: 利用盘片水平瞬间求出入射光线的方向向量 假设中，入射光线的方向向量是已知的，并不是十分合理，可通过下述方法求得其方向向量。如图8，控制h=0，取振动初始情况进行研究，此时刻反射面为水平，即Z=D。通过B点、O点的坐标及反射定理即可方便的求解。 图8水平反射图 6.2.2改进方向二: 找到接受屏上的读数点坐标与入射光线的对应关系 对于通过处理大量数据和多种情况讨论来解决问题的题目，代数法无疑是一种很好的选择，本题中若是已知B点及入射光线的对应关系，则可以快速的假设出O点、B点坐标，带入求出的式便可解出c的值，比较简单。 9 6.2.2改进方向三:问题三的求解中应考虑分辨率对c的值的影响及n值的有限性 硬盘在振动过程中，实际上是不规则的倾斜振动，计算过程中应该同样求出一系列不同分辨率下的c的值，并通过最小二乘法求出最优分辨率，这样将使结果更加准确，更具有可信性。 此外，n的取值具有受限性，题中只模拟了5个不同分辨率对于a、b、c结果的影响就得出结论，其实是不够充分的，应该多取几个不同的n值来模拟，方能使结论更具说服力。 参考文献 1 全国大学生数学建模竞赛组委会,全国大学生数学建模竞赛优秀论文汇编,中国物价出版社, 2002.3. 2 王庚 王敏生,现代数学建模方法,科学出版社,2008. 3 同济大学应用数学系, 高等数学, 北京: 高等教育出版社, 2001. 10 附录 c=1/2*p1*(-a2*x12*n12+m12*a2*h2+m12*a2*y12+n12*b2*x12-b2*y12*m12+n12*b2*h2-a2*x12*p12-b2*y12*p12-2*m1*a*p1*x12-2*m1*a*p1*y12-2*m1*a*p1*h2-2*n1*b*p1*x12-2*n1*b*p1*y12-2*n1*b*p1*h2+2*b*y1*h*m12+2*b*y1*h*n12+2*b*y1*h*p12+p12*x12+p12*y12-h2*m12-h2*n12+2*m1*a*n1*b*x12+2*m1*a*n1*b*y12+2*m1*a*n1*b*h2+2*a*x1*h*m12+2*a*x1*h*n12+2*a*x1*h*p12-2*a*x1*b*y1*m12-2*a*x1*b*y1*n12-2*a*x1*b*y1*p12)/(-a*m1*b*y1*p12-p1*h*m12+p12*m1*x1+p12*n1*y1-p1*h*n12+a*x1*p13+b*y1*p13+a*m13*h+b*n13*h+m12*a2*n1*y1+m12*a2*p1*h+m12*a*n1*b*x1-2*m1*a*p1*n1*y1+m1*a*n12*b*y1+2*m1*a*n1*b*p1*h-m12*a*p1*x1-2*n1*b*p1*m1*x1-m1*a*p12*h+n12*b2*m1*x1+n12*b2*p1*h-a*x1*b*n13-n12*b*p1*y1-n1*b*p12*h-a*m13*b*y1-b2*n1*y1*m12-b2*n1*y1*p12-a*x1*b*n1*p12-a2*m1*x1*n12-a2*m1*x1*p12+a*x1*p1*n12+b*y1*p1*m12+a*m1*h*n12+b*n1*h*m12)-1/2*(p1*x1-m1*h)*(m2*y2-n2*x2)-(m1*y1-n1*x1)*(p2*x2-m2*h)/(-(n1*h-p1*y1)*(p2*x2-m2*h)+(p1*x1-m1*h)*(n2*h-p2*y2)*m1*(-a2*x12*n12+m12*a2*h2+m12*a2*y12+n12*b2*x12-b2*y12*m12+n12*b2*h2-a2*x12*p12-b2*y12*p12-2*m1*a*p1*x12-2*m1*a*p1*y12-2*m1*a*p1*h2-2*n1*b*p1*x12-2*n1*b*p1*y12-2*n1*b*p1*h2+2*b*y1*h*m12+2*b*y1*h*n12+2*b*y1*h*p12+p12*x12+p12*y12-h2*m12-h2*n12+2*m1*a*n1*b*x12+2*m1*a*n1*b*y12+2*m1*a*n1*b*h2+2*a*x1*h*m12+2*a*x1*h*n12+2*a*x1*h*p12-2*a*x1*b*y1*m12-2*a*x1*b*y1*n12-2*a*x1*b*y1*p12)/(-a*m1*b*y1*p12-p1*h*m12+p12*m1*x1+p12*n1*y1-p1*h*n12+a*x1*p13+b*y1*p13+a*m13*h+b*n13*h+m12*a2*n1*y1+m12*a2*p1*h+m12*a*n1*b*x1-2*m1*a*p1*n1*y1+m1*a*n12*b*y1+2*m1*a*n1*b*p1*h-m12*a*p1*x1-2*n1*b*p1*m1*x1-m1*a*p12*h+n12*b2*m1*x1+n12*b2*p1*h-a*x1*b*n13-n12*b*p1*y1-n1*b*p12*h-a*m13*b*y1-b2*n1*y1*m12-b2*n1*y1*p12-a*x1*b*n1*p12-a2*m1*x1*n12-a2*m1*x1*p12+a*x1*p1*n12+b*y1*p1*m12+a*m1*h*n12+b*n1*h*m12)-1/2*(m1*y1-n1*x1)*(n2*h-p2*y2)-(n1*h-p1*y1)*(m2*y2-n2*x2)/(-(n1*h-p1*y1)*(p2*x2-m2*h)+(p1*x1-m1*h)*(n2*h-p2*y2)*n1*(-a2*x12*n12+m12*a2*h2+m12*a2*y12+n12*b2*x12-b2*y12*m12+n12*b2*h2-a2*x12*p12-b2*y12*p12-2*m1*a*p1*x12-2*m1*a*p1*y12-2*m1*a*p1*h2-2*n1*b*p1*x12-2*n1*b*p1*y12-2*n1*b*p1*h2+2*b*y1*h*m12+2*b*y1*h*n12+2*b*y1*h*p12+p12*x12+p12*y12-h2*m12-h2*n12+2*m1*a*n1*b*x12+2*m1*a*n1*b*y12+2*m1*a*n1*b*h2+2*a*x1*h*m12+2*a*x1*h*n12+2*a*x1*h*p12-2*a*x1*b*y1*m12-2*a*x1*b*y1*n12-2*a*x1*b*y1*p12)/(-a*m1*b*y1*p12-p1*h*m12+p12*m1*x1+p12*n1*y1-p1*h*n12+a*x1*p13+b*y1*p13+a*m13*h+b*n13*h+m12*a2*n1*y1+m12*a2*p1*h+m12*a*n1*b*x1-2*m1*a*p1*n1*y1+m1*a*n12*b*y1+2*m1*a*n1*b*p1*h-m12*a*p1*x1-2*n1*b*p1*m1*x1-m1*a*p12*h+n12*b2*m1*x1+n12*b2*p1*h-a*x1*b*n13-n12*b*p1*y1-n1*b*p12*h-a*m13*b*y1-b2*n1*y1*m12-b2*n1*y1*p12-a*x1*b*n1*p12-a2*m1*x1*n12-a2*m1*x1*p12+a*x1*p1*n12+b*y1*p1*m12+a*m1*h*n12+b*n1*h*m12)-D syms m1 n1 p1 m2 n2 p2 %oe1方向向量 已知 syms x1 y1 h x2 y2 %B1 B2坐标 已知 syms x11 y11 z11 x21 y21 z21%E1 E2 坐标 未知 syms a b c L D%所求参数 L为法线 未知 syms Loe1 Loe2 Leb1 Leb2 %各条线 L=a b -1 11 Loe1=m1 n1 p1; Lob1=x1 y1 h; Loe2=m2 n2 p2; Lob2=x2 y2 h; f=cross(cross(Loe1,Lob1),cross(Loe2,Lob2); a=f(1)/-f(3); b=f(2)/-f(3); syms d; l1='(m1*a+n1*b-p1)2/(m12+n12+p12)=(a*(m1*d-x1)+b*(n1*d-y1)-(p1*d-h)2/(m1*d-x1)2+(n1*d-y1)2+(p1*d-h)2)' d=solve(l1,d); x11=m1*d; y11=n1*d; z11=p1*d; c=z11-a*x11-b*y11-D; syms d1 d2 m n o b1f b2f; %h D l1='D/cos(m-o)=n/cos(o)' l2='n/cos(m-2*o)=d1/sin(2*(m-o)' l3='b1f*cot(m-2*o)=h' b1f d1 n=solve(l1,l2,l3,b1f,d1,n); b2f=b1f; d2=d1; y=rand(10); b1f=vpa(eval(subs(b1f,m,o,D,h,pi/4,5000/4455*10(-6)*y,-0.05,0.01); d1=vpa(eval(subs(d1,m,o,D,h,pi/4,5000/4455*10(-6)*y,-0.05,0.01); ans1=b1f+d1; %x1 y1 b2f=vpa(eval(subs(b2f,m,o,D,h,pi/3,5000/4455*10(-6)*y,-0.05,0.01); d2=vpa(eval(subs(d2,m,o,D,h,pi/3,5000/4455*10(-6)*y,-0.05,0.01); ans2=b2f+d2; %x2 y2 ans14=vpa(round(10000*ans1)/10000); ans15=vpa(round(100000*ans1)/100000); ans16=vpa(round(1000000*ans1)/1000000); ans17=vpa(round(10000000*ans1)/10000000); ans18=vpa(round(100000000*ans1)/100000000); ans24=vpa(round(10000*ans2)/10000); ans25=vpa(round(100000*ans2)/100000); ans26=vpa(round(1000000*ans2)/1000000); ans27=vpa(round(10000000*ans2)/10000000); ans28=vpa(round(100000000*ans2)/100000000); 12 %a1=vpa(subs(a,x1,y1,x2,y2,m1,n1,p1,m2,n2,p2,h,D,ans1,ans1,ans2,ans2,1,1,1,-sqrt(3),-sqrt(3),1,0.01-0.05) a10=vpa(subs(a,m1,n1,p1,m2,n2,p2,h,D,1,1,1,-sqrt(3),-sqrt(3),1,0.01,-0.05) a10=collect(numden(a10) a10=collect(vpa(subs(a10,y1,y2,x1,x2) a1=subs(a10,x1,x2,ans1,ans2) a14=subs(a10,x1,x2,ans14,ans24) a15=subs(a10,x1,x2,ans15,ans25) a16=subs(a10,x1,x2,ans16,ans26) a17=subs(a10,x1,x2,ans17,ans27) a18=subs(a10,x1,x2,ans18,ans28) b10=vpa(subs(b,m1,n1,p1,m2,n2,p2,h,D,1,1,1,-sqrt(3),-sqrt(3),1,0.01,-0.05) b10=collect(numden(b10) b10=collect(vpa(subs(b10,y1,y2,x1,x2) b1=subs(b10,x1,x2,ans1,ans2) b14=subs(b10,x1,x2,ans14,ans24) b15=subs(b10,x1,x2,ans15,ans25) b16=subs(b10,x1,x2,ans16,ans26) b17=subs(b10,x1,x2,ans17,ans27) b18=subs(b10,x1,x2,ans18,ans28) a24=sum(sum(a14-a1).2) a25=sum(sum(a15-a1).2) a26=sum(sum(a16-a1).2) a27=sum(sum(a17-a1).2) a28=sum(sum(a18-a1).2) b24=sum(sum(b14-b1).2) b25=sum(sum(b15-b1).2) b26=sum(sum(b16-b1).2) b27=sum(sum(b17-b1).2) b28=sum(sum(b18-b1).2) e=double(a24+b24)/100) double(a25+b25)/100) double(a26+b26)/100) double(a28+b28)/100); nn=4:8;bar(nn,e) 13 double(a24+b27)/100)