光學(xué)原理測量硬盤振動(dòng)參考模板

1、幾何光學(xué)測量硬盤振動(dòng)模型摘 要本文通過理論分析和與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比, 研究了通過幾何光學(xué)測量硬盤振動(dòng)問題, 并求出了硬盤盤面的局部表面方程，分析了不同分辨率對(duì)局部表面方程的影響，較好的解釋了硬盤振動(dòng)的規(guī)律。首先研究不考慮盤片中心的上下振動(dòng)的情形. 針對(duì)問題做出了幾何模型, 利用所給數(shù)據(jù),并參照查閱資料，最終求出考慮盤片中心上下振動(dòng)情況下的結(jié)果，并對(duì)求解結(jié)果采用代數(shù)檢驗(yàn)符合要求, 提出模型的改進(jìn)方向。對(duì)于問題一, 查閱了硬盤的工作原理和內(nèi)部結(jié)構(gòu)，從幾何上進(jìn)行理論分析，討論了一條和多條光線對(duì)模型測量的影響，通過分析光路、成本、技術(shù)等因素得出結(jié)論采用兩束光線的方案為最優(yōu)。對(duì)于問題二, 首先建立了不考慮

2、盤片中心的上下振動(dòng)的理想模型，通過幾何關(guān)系和數(shù)學(xué)定理,順利求出了盤片局部表面方程的參數(shù)a、b，并得到c=0；其次在前面計(jì)算基礎(chǔ)上建立了考慮盤片中心的上下振動(dòng)的理想模型，通過角度關(guān)系有效地解出了c的值（見附錄）。對(duì)于問題三，通過分析得知，屏幕分辨率對(duì)于測量結(jié)果的影響其實(shí)也就是對(duì)于a、b參數(shù)值的影響。利用MATLAB生成幾組反射面傾斜角度的隨機(jī)值，可以求出一系列和在接受屏上的坐標(biāo)值，再用fix函數(shù)對(duì)該值進(jìn)行某一精度（，等等）下的四舍五入取整，這樣獲得的盤片局部表面方程就包含了接受屏分辨率的影響，再通過最小二乘法求得實(shí)際值與改變值差的平方最小的值，其對(duì)應(yīng)的分辨率為最佳分辨率。在結(jié)果的分析與檢驗(yàn)中,

3、改變數(shù)據(jù)與實(shí)際數(shù)據(jù)的對(duì)比圖, 直觀地反映出分辨率對(duì)于測量結(jié)果的影響，驗(yàn)證了模型的正確性與方法的可靠性，但兩者之間仍存在偏差, 對(duì)此我們作出了合理的解釋。在模型的改進(jìn)中, 著重分析了文章前一部分未考慮在內(nèi)的各種不定因素對(duì)硬盤振動(dòng)的影響, 并提出了兩種改進(jìn)方案：第一種利用盤片水平瞬間求出入射光線的方向向量；第二種是找到接受屏上的讀數(shù)點(diǎn)坐標(biāo)與入射光線的對(duì)應(yīng)關(guān)系，便可通過代值法求的參數(shù)值；第三種是問題三的求解中應(yīng)考慮分辨率對(duì)c的值的影響及n值的有限性。關(guān)鍵詞: 硬盤振動(dòng) 幾何光學(xué) 局部表面方程 分辨率 最小二乘 1問題的重述硬盤是計(jì)算機(jī)上的重要部件，正向著更小、更快、密度更高的方向發(fā)展。如何減小盤片的

4、振動(dòng)，成為關(guān)鍵問題，于是怎樣檢測盤片的振動(dòng)更加重要。由于硬盤轉(zhuǎn)速很高，不適宜接觸式測量，容易想到光學(xué)測量。2 / 18有一種新技術(shù)測量物體表面獲得振動(dòng)情況，是利用類似鏡面反射的幾何光學(xué)原理工作的：選用同一點(diǎn)光源發(fā)出兩束光線，照射有反射能力的被測物體平面，產(chǎn)生的兩條反射光線再照射到水平放置的接受屏上。接受屏能夠檢測出光線照射點(diǎn)在接受屏平面上的坐標(biāo)位置，問題是怎樣利用圖2中Oxy平面的兩個(gè)坐標(biāo)值反算出被測物體當(dāng)前的平面方程，就可利用平面方程的變化得到振動(dòng)情況。圖1是用兩條光線檢測的原理圖，點(diǎn)光源O為坐標(biāo)原點(diǎn)，E1和E2是盤片瞬間位置的反射點(diǎn)，此時(shí)盤片局部表面方程為z=ax+by+(c+D)，其中a

5、,b,c為未知量，由于硬盤在振動(dòng)，平面就偏離了初始位置z=D(D<0)；接受屏方程為z=h，其中h>0為接受屏的高度,B1和B2是其上的照射點(diǎn)，這樣利用它們確定、和獲得盤片局部表面方程。（本題作到確定出、和即可，不必分析振動(dòng)情況）問題 1：請(qǐng)解釋這種測量技術(shù)采用兩束光線的原因。問題 2：請(qǐng)建數(shù)學(xué)模型分析盤片局部表面方程的測量原理。問題 3：試分析接受屏分辨率（測量精度）對(duì)盤片局部表面方程的影響。提示：可利用計(jì)算機(jī)模擬的辦法求出一系列B1和B2在接受屏上的坐標(biāo)值，再對(duì)該值進(jìn)行某一精度（10-4，10-5，10-6等等）下的四舍五入取整，這樣獲得的盤片局部表面方程就包含了接受屏分辨率的

6、影響。2問題的分析21 問題一根據(jù)硬盤的工作原理和內(nèi)部結(jié)構(gòu)可以知道，硬盤工作時(shí)盤片時(shí)刻在振動(dòng)，不同時(shí)刻對(duì)應(yīng)有不同的方程，若采用一條光線則無法求解，同時(shí)從成本、技術(shù)出發(fā)，兩條光線為最佳方案。被測物體平面接受屏圖1 原理圖點(diǎn)光源E1E2B2B1zxy22 問題二OB1圖2 接受屏示意圖B2x ,y ,二光束法測量硬盤振動(dòng)時(shí)，兩條入射光線與其反射光線所形成的平面均垂直于反射面，又經(jīng)過共同的點(diǎn)光源，且過點(diǎn)有且只有一條直線垂直于反射面，故可以求出反射面的法線。根據(jù)鏡面反射原理，入射角等于反射角，同時(shí)已知、法線、的方向向量，可以求出點(diǎn)的坐標(biāo)，利用點(diǎn)法式可得出平面方程，即a、b、c的值。23問題三通過分析得

7、知，屏幕分辨率對(duì)于測量結(jié)果的影響其實(shí)也就是對(duì)于a、b參數(shù)值的影響。利用MATLAB生成幾組入射角的隨機(jī)值，可以求出一系列和在接受屏上的坐標(biāo)值，再用fix函數(shù)對(duì)該值進(jìn)行某一精度（，等等）下的四舍五入取整，這樣獲得的盤片局部表面方程就包含了接受屏分辨率的影響，再通過最小二乘法求得實(shí)際值與改變值差的平方最小的值，其對(duì)應(yīng)的分辨率為最佳分辨率。3模型的假設(shè)與符號(hào)說明31 模型的假設(shè)(1)兩條入射光線的方向向量為已知，且固定不變; (2)接受屏上的坐標(biāo)原點(diǎn)與點(diǎn)所形成的直線與軸重合；(3)硬盤盤片振動(dòng)過程中不會(huì)發(fā)生變形，局部表面始終為平面；(4)光線發(fā)生反射時(shí)為全反射，且傳播過程中沒有干涉、衍射等現(xiàn)象；(5

8、)硬盤平放，盤片振動(dòng)連續(xù)、穩(wěn)定。32 符號(hào)的說明 表示直線的方向向量; 表示直線的方向向量; 表示點(diǎn)的坐標(biāo); 表示點(diǎn)的坐標(biāo); 表示點(diǎn)的坐標(biāo); 表示點(diǎn)的坐標(biāo); 表示面的法線; 表示面的法線; 表示反射面的法線的方向向量。4模型的建立與求解41求解局部平面方程根據(jù)幾何定理得知，若平面A與平面B均垂直于平面C，假設(shè)AB相交于直線L，則L必垂直于平面C。又知入射光線與反射光線所形成的平面與反射面垂直，故題中的入射問題可簡化為圖3，如圖。圖3 光線入射根據(jù)數(shù)學(xué)定理，面和面的法線分別為：= （1）= （2）由于兩個(gè)法線相交所形成的平面與反射面平行，故可以求得反射面的法線的方向向量：= （3）因?yàn)榉瓷涿娴姆?/p>

9、程為z=ax+by+(c+D)，因此其法向量為，又兩平行向量差乘的絕對(duì)值等于0，即： （4）由上便可求出未知參數(shù)a、b， （5） （6）通過上面的求解過程可以得知用一束光無法得出局部平面方程，故題中選用二束光，下面求c的值。41.1不考慮盤片中心的上下振動(dòng)硬盤振動(dòng)過程中盤片中心始終固定在主軸上，即振動(dòng)平面均通過點(diǎn)，帶入反射面方程即可求得c的值，結(jié)果為： （7）41.2考慮盤片中心的上下振動(dòng)在圖4中， 為入射光線，為法線，為反射光線，根據(jù)反射定理，入射角等于反射角，故,由題、的方向向量分別為：、 ，根據(jù)空間兩直線的夾角公式得： （8）圖 4 等角示意圖此外，在直線上，同時(shí)也在反射平面上，即： （

10、9） （10）由方程（8）（9）（10）可求的c的值，由于c值較長，詳情見附錄1。42接受屏分辨率對(duì)局部表面方程的影響42.1問題的分析盤片振動(dòng)過程中以偏離原平面的角度為考察對(duì)象，對(duì)于盤片表面上任一點(diǎn)，該點(diǎn)的z值反映了該點(diǎn)偏離原平面的情況。而z=ax+by+(c+D)，所以接受屏分辨率對(duì)局部表面方程的影響即是考察不同分辨率下z值的變化量的情況。即對(duì)其表面方程系數(shù)a，b的影響。盤片時(shí)刻在振動(dòng)，設(shè)其偏離初始狀態(tài)z=D的角度為，當(dāng)取不同值時(shí)，得到不同的，對(duì)取不同的精度，并在這些精度下四舍五入取整，最后利用MATLAB軟件求解，通過a，b前后變化來分析接受屏分辨率（測量精度）對(duì)盤片局部表面方程的影響。

11、42.2問題的求解圖5盤面振動(dòng)解析圖依據(jù)正弦定理，可得到接受屏上點(diǎn)在坐標(biāo)上的計(jì)算公式，如下： （11）據(jù)題意知道點(diǎn)坐標(biāo)為 （12）其中為入射光線在平面OXZ或平面OYZ上的投影與OZ軸的夾角，當(dāng)在不同平面上取不同值時(shí)，就可求得相應(yīng)的在接受屏上的反射光線點(diǎn)的位置。例如可分別取入射光線在平面OXZ或平面OYZ上的投影與OZ軸的夾角，同時(shí)令，則求得接受屏上點(diǎn)在O-XYZ內(nèi)的坐標(biāo)，同理可得點(diǎn)的坐標(biāo)。然后，對(duì)點(diǎn)的坐標(biāo)取不同的精度，并在這些精度下四舍五入取整，得到五組不同分辨率下的點(diǎn)的坐標(biāo)根據(jù)問題二可求的相應(yīng)的五組a、b值，利用最小二乘法從中得到最優(yōu)的分辨率。42.3算法的實(shí)現(xiàn)（1）由資料得知，普通盤片的

12、振動(dòng)半徑為44.5mm，振動(dòng)范圍小于50納米，故得到的取值范圍為，用MATLAB在取值范圍內(nèi)隨機(jī)抽取100組值,見表一。表一 100組的隨機(jī)值0.08910.09430.07670.08940.10060.02600.08400.05590.03480.00130.06920.03970.00220.10120.01540.00650.00170.02180.08180.10760.10690.05000.10050.02240.07560.02130.06680.06010.06990.08290.03840.07810.09850.00730.11120.09600.03430.02280

13、.03060.04610.09240.10310.06790.05910.04820.09670.03260.05000.02840.05280.03200.03050.00210.10030.08570.05130.10520.02240.09520.07980.04260.06150.03850.01980.07440.07430.07900.01950.03360.11020.05660.01700.07850.09200.04170.04710.02280.04560.05240.09360.01560.05470.00110.04270.07250.06390.08300.06090

14、.10480.09900.06830.09150.10370.07660.02240.09430.09240.03410.05010.0588（2）當(dāng)時(shí)，對(duì)每一個(gè)值根據(jù)公式（11）（12）即可求得點(diǎn)坐標(biāo)（見附錄），其中根據(jù)公式（12）。同理當(dāng)時(shí)，可求得點(diǎn)的坐標(biāo),且。（3）根據(jù)的坐標(biāo)，由公式（5）（6）可以求得a、b的準(zhǔn)確值，再把取不同的精度，分別為：，并在這些精度下四舍五入取整，分別求得不同精度下。最后利用MATLAB軟件，求得相應(yīng)的五組a、b值，將得到的五組a、b的值與未經(jīng)四舍五入的a、b的值做差的平方運(yùn)算，再取平均數(shù)，比較其大小即可得出結(jié)果，選出最優(yōu)的分辨率，見表二。表二 最小二乘法結(jié)果n

15、 = 4n = 5n = 6n = 7n = 85模型結(jié)果的分析以n為橫坐標(biāo)，為縱坐標(biāo)作圖如下，當(dāng)n=8時(shí)所對(duì)應(yīng)的分辨率對(duì)結(jié)果的影響最小，并且可以作出推測：隨著n的不斷增大，即分辨率越高，所求出的結(jié)果就越準(zhǔn)確。圖 6 圖 76模型評(píng)價(jià)與改進(jìn)61 模型的評(píng)價(jià)6.1.1 模型的優(yōu)點(diǎn)(1) 問題2構(gòu)建模型過程中，充分考慮了盤片中心是否上下振動(dòng)兩種情況，使討論更加全面，結(jié)論更加可靠。(2)問題3構(gòu)建模型時(shí)，查閱了大量資料，得出了準(zhǔn)確有效的已知條件，為后面的計(jì)算提供了保障。(3) 采用圖示、表格的方法直觀地反映了所建模型求解數(shù)據(jù)與實(shí)測數(shù)據(jù)的對(duì)應(yīng)關(guān)系; (4) 模型具有全面性，容易推廣。 6.1.2 模型

16、的缺點(diǎn)(1) 該模型數(shù)據(jù)處理比較復(fù)雜, 運(yùn)算較為繁瑣; (2) 由于數(shù)據(jù)的不充足, 對(duì)某些參數(shù)的確定具有主觀成分;(3) 沒有考慮干擾因素的影響, a、b、c的值仍然存在誤差。62 模型的改進(jìn)方向6.2.1改進(jìn)方向一: 利用盤片水平瞬間求出入射光線的方向向量假設(shè)中，入射光線的方向向量是已知的，并不是十分合理，可通過下述方法求得其方向向量。如圖8，控制h=0，取振動(dòng)初始情況進(jìn)行研究，此時(shí)刻反射面為水平，即Z=D。通過B點(diǎn)、O點(diǎn)的坐標(biāo)及反射定理即可方便的求解。圖8水平反射圖6.2.2改進(jìn)方向二: 找到接受屏上的讀數(shù)點(diǎn)坐標(biāo)與入射光線的對(duì)應(yīng)關(guān)系對(duì)于通過處理大量數(shù)據(jù)和多種情況討論來解決問題的題目，代數(shù)法

17、無疑是一種很好的選擇，本題中若是已知B點(diǎn)及入射光線的對(duì)應(yīng)關(guān)系，則可以快速的假設(shè)出O點(diǎn)、B點(diǎn)坐標(biāo)，帶入求出的（5）（6）式便可解出c的值，比較簡單。6.2.2改進(jìn)方向三:問題三的求解中應(yīng)考慮分辨率對(duì)c的值的影響及n值的有限性 硬盤在振動(dòng)過程中，實(shí)際上是不規(guī)則的傾斜振動(dòng)，計(jì)算過程中應(yīng)該同樣求出一系列不同分辨率下的c的值，并通過最小二乘法求出最優(yōu)分辨率，這樣將使結(jié)果更加準(zhǔn)確，更具有可信性。此外，n的取值具有受限性，題中只模擬了5個(gè)不同分辨率對(duì)于a、b、c結(jié)果的影響就得出結(jié)論，其實(shí)是不夠充分的，應(yīng)該多取幾個(gè)不同的n值來模擬，方能使結(jié)論更具說服力。參考文獻(xiàn)1 全國大學(xué)生數(shù)學(xué)建模競賽組委會(huì),全國大學(xué)生數(shù)學(xué)

18、建模競賽優(yōu)秀論文匯編,中國物價(jià)出版社, 2002.3. 2 王庚 王敏生,現(xiàn)代數(shù)學(xué)建模方法,科學(xué)出版社,2008.3 同濟(jì)大學(xué)應(yīng)用數(shù)學(xué)系, 高等數(shù)學(xué), 北京: 高等教育出版社, 2001. 附錄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*n1

19、2+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*n

20、1*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-

21、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+

22、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*p

23、1*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

24、*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

25、-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*n1

26、2*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)-Dsyms m1 n1 p1 m2 n2 p2 %oe1方向向量 已知syms x1 y1 h x2 y2

27、%B1 B2坐標(biāo) 已知syms x11 y11 z11 x21 y21 z21%E1 E2 坐標(biāo) 未知syms a b c L D%所求參數(shù) L為法線 未知syms Loe1 Loe2 Leb1 Leb2 %各條線L=a b -1Loe1=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

28、*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 Dl1='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

29、(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 y1b2f=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

30、 y2ans14=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);%a1=vpa(subs(a,x1,y1,x2,y2,m1,n1,p1,m2,n2,p2,h,D,ans1,ans1,ans2,a