計算橢圓運動軌跡的演算法

標籤：cross   class   大學   code   ima   div   unit   根據   貝茲路徑   

在中學或大學時代，我們應該在數學中都學過橢圓方程、雙曲線方程等等，當然那個時候學習這些知識的目的就是為了考試，為了能夠拿個好成績，上個好大學等。那麼除此之外，這些知識對於我們今後的生活或者工作又帶來什麼便利呢？

 

巧合的是作為程式員，尤其是偏向演算法方面的人員，經常會有這種需求，例如實現一個物體做橢圓運動的效果，或者做一個圓形軌跡的效果，亦或者做一個看似毫無規律的曲線運動，那麼這就用到了橢圓方程、圓形方程及貝茲路徑方程等等。

 

在此著重介紹下橢圓方程的應用。

 

一、 標準橢圓方程公式：

二、 中心點在 ( h , k )，主軸平行於x軸時公式：

三、 常見的名詞概念及橢圓形狀：

       

　　　　　　　　　　圖1　　　　　　　　　　　　　　　　　　　　　　　　圖2　　　　　　

• 長軸：通過串連橢圓上的兩個點所能獲得的最長線段，對應圖2中的major axis。
• 半長軸：長軸的一半，對應橢圓公式中的a。
• 短軸：垂直平分長軸的直線所得弦， 對應圖2中的minor axis。
  
• 半短軸：短軸的一半，對應橢圓公式中的b。
• 焦點：每個橢圓均有兩個焦點，分別上中的F1 , F2。
• 近拱點：指定一個焦點F1 , 距離焦點F1最近點成為近拱點，也就是圖中的A。
• 遠拱點：相對於F1焦點而言，距離F1最遠點成為遠拱點，也就是圖中的B。
• 離心率：用來描述軌道的形狀，橢圓的離心率大小區間(0, 1) ,當離心率為0時表示圓。

四、 離心率的計算公式：

• 根據近拱點和遠拱點距離計算：

其中dp表示近拱點的距離，da表示遠拱點的距離。

• 根據半焦距和半長軸計算：

 

其中c表示半焦距，a表示半長軸。

• 根據半長軸和半短軸計算：

其中，a表示半長軸，b表示半短軸。

 

五、 已知橢圓其中一個焦點、近拱點、離心率以及表示橢圓在三維空間的傾向法向量，可以得到此橢圓的公式。

1. 根據焦點和近拱點，得到近拱點到焦點的距離dp。

2. 根據離心率、近拱點距離可以得到遠拱點到焦點的距離da。

3. 根據近拱點、焦點以及遠拱點距離，可以得到遠拱點的座標值。

4. 根據近拱點、遠拱點即可得到橢圓的中心點座標center。

5. 根據公式：半短軸 = 近拱點距離 * 遠拱點距離，即可得到半長軸b的長度。

6. 根據代表橢圓傾向的法向量n和長軸方向的單位向量，兩者叉乘即可得到短軸方向上的單位向量，並且根據半短軸的長度，即可得到半短軸向量。

7. 已知半短軸向量、半長軸向量，根據橢圓的參數方程，即可得到橢圓上任一點的座標值。

 

具體代碼如下：

 1 //mDirection代表點在橢圓上的運動方向是逆時針還是順時針，angle用於計算橢圓參數方程的角度 2 double angle = (mDirection == OrbitDirection.CLOCKWISE ? -1 : 1) * mAngle * time * MathUtil.PRE_PI_DIV_180; 3  4 //計算近拱點到焦點的距離 5 double periapsisRadius = mPeriapsis.distanceTo(mFocalPoint); 6 //根據離心率、近拱點到焦點的距離、遠拱點到焦點三者的公式即可得到遠拱點距離(近拱點到遠拱點的距離是橢圓的長軸) 7 double apoapsisRadius = periapsisRadius * (1 + mEccentricity) / (1 - mEccentricity); 8  9 // 計算近拱點到焦點的單位向量，注意此處乘以1e8 和除以1e8的目的是 去掉第8個小數點後最不重要的數字，以降低計算誤差。10 double uAx = (Math.round(mFocalPoint.x * 1e8) - Math.round(mPeriapsis.x * 1e8)) / 1e8;11 double uAy = (Math.round(mFocalPoint.y * 1e8) - Math.round(mPeriapsis.y * 1e8)) / 1e8;12 double uAz = (Math.round(mFocalPoint.z * 1e8) - Math.round(mPeriapsis.z * 1e8)) / 1e8;13 double mod = Math.sqrt(uAx * uAx + uAy * uAy + uAz * uAz);//計算近拱點到焦點的距離14 if (mod != 0 && mod != 1) {//單位化15     mod = 1 / mod;16     uAx *= mod;17     uAy *= mod;18     uAz *= mod;19 }20 21 double apoapsisDir_x = Math.round(uAx * apoapsisRadius * 1e8) / 1e8;22 double apoapsisDir_y = Math.round(uAy * apoapsisRadius * 1e8) / 1e8;23 double apoapsisDir_z = Math.round(uAz * apoapsisRadius * 1e8) / 1e8;24 //計算遠拱點座標25 double apoapsisPos_x = Math.round((apoapsisDir_x + mFocalPoint.x) * 1e8) / 1e8;26 double apoapsisPos_y = Math.round((apoapsisDir_y + mFocalPoint.y) * 1e8) / 1e8;27 double apoapsisPos_z = Math.round((apoapsisDir_z + mFocalPoint.z) * 1e8) / 1e8;28 29 //近拱點和遠拱點的中心即橢圓的中心30 double center_x = Math.round(((mPeriapsis.x + apoapsisPos_x) / 2) * 1e8) / 1e8;31 double center_y = Math.round(((mPeriapsis.y + apoapsisPos_y) / 2) * 1e8) / 1e8;32 double center_z = Math.round(((mPeriapsis.z + apoapsisPos_z) / 2) * 1e8) / 1e8;33 34 //計算半短軸的長度35 double b = Math.sqrt(periapsisRadius * apoapsisRadius);36 37 //從中心點到近拱點的向量38 double semimajorAxis_x = Math.round((mPeriapsis.x - center_x) * 1e8) / 1e8;39 double semimajorAxis_y = Math.round((mPeriapsis.y - center_y) * 1e8) / 1e8;40 double semimajorAxis_z = Math.round((mPeriapsis.z - center_z) * 1e8) / 1e8;41 42 //單位化半長軸向量43 double unitSemiMajorAxis_x = semimajorAxis_x;44 double unitSemiMajorAxis_y = semimajorAxis_y;45 double unitSemiMajorAxis_z = semimajorAxis_z;46 mod = Math.sqrt(semimajorAxis_x * semimajorAxis_x + semimajorAxis_y * semimajorAxis_y + semimajorAxis_z47         * semimajorAxis_z);48 if (mod != 0 && mod != 1) {49     mod = 1 / mod;50     unitSemiMajorAxis_x *= mod;51     unitSemiMajorAxis_y *= mod;52     unitSemiMajorAxis_z *= mod;53 }54 55 //將中心點沿法向量方向平移56 Vector3 unitNormal = mNormal.clone();57 unitNormal.normalize();58 double uNx = Math.round(unitNormal.x * 1e8) / 1e8;59 double uNy = Math.round(unitNormal.y * 1e8) / 1e8;60 double uNz = Math.round(unitNormal.z * 1e8) / 1e8;61 double normalCenter_x = center_x + uNx;62 double normalCenter_y = center_y + uNy;63 double normalCenter_z = center_z + uNz;64 mod = Math.sqrt(normalCenter_x * normalCenter_x + normalCenter_y * normalCenter_y + normalCenter_z65         * normalCenter_z);66 if (mod != 0 && mod != 1) {67     mod = 1 / mod;68     normalCenter_x *= mod;69     normalCenter_y *= mod;70     normalCenter_z *= mod;71 }72 73 mScratch1.setAll(unitSemiMajorAxis_x, unitSemiMajorAxis_y, unitSemiMajorAxis_z);74 mScratch2.setAll(normalCenter_x, normalCenter_y, normalCenter_z);75 //叉乘計算半短軸的單位向量76 Vector3 semiminorAxis = mScratch3.crossAndSet(mScratch1, mScratch2);77 //得到半短軸向量78 semiminorAxis.multiply(b);79 80 //3D空間橢圓的參數方程81 double x = center_x + (Math.cos(angle) * semimajorAxis_x) + (Math.sin(angle) * semiminorAxis.x);82 double y = center_y + (Math.cos(angle) * semimajorAxis_y) + (Math.sin(angle) * semiminorAxis.y);83 double z = center_z + (Math.cos(angle) * semimajorAxis_z) + (Math.sin(angle) * semiminorAxis.z);

 

計算橢圓運動軌跡的演算法