计算椭圆运动轨迹的算法

在中学或大学时代，我们应该在数学中都学过椭圆方程、双曲线方程等等，当然那个时候学习这些知识的目的就是为了考试，为了能够拿个好成绩，上个好大学等。那么除此之外，这些知识对于我们今后的生活或者工作又带来什么便利呢？

巧合的是作为程序员，尤其是偏向算法方面的人员，经常会有这种需求，例如实现一个物体做椭圆运动的效果，或者做一个圆形轨迹的效果，亦或者做一个看似毫无规律的曲线运动，那么这就用到了椭圆方程、圆形方程及贝塞尔曲线方程等等。

在此着重介绍下椭圆方程的应用。

一、 标准椭圆方程公式：

二、 中心点在 ( 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_z
47         * 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_z
65         * 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);