计算椭圆运动轨迹的算法
在中学或大学时代，我们应该在数学中都学过椭圆方程、双曲线方程等等，当然那个时候学习这些知识的目的就是为了考试，为了能够拿个好成绩，上个好大学等。那么除此之外，这些知识对于我们今后的生活或者工作又带来什么便利呢？
巧合的是作为程序员，尤其是偏向算法方面的人员，经常会有这种需求，例如实现一个物体做椭圆运动的效果，或者做一个圆形轨迹的效果，亦或者做一个看似毫无规律的曲线运动，那么这就用到了椭圆方程、圆形方程及贝塞尔曲线方程等等。
在此着重介绍下椭圆方程的应用。
一、 标准椭圆方程公式：
二、 中心点在 ( 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);
