这次呢，主要是解决上次遗留的问题之一：行星运动
为了找这个行星运动方程，我翻遍了各个百科的角落旮旯，终于在果壳的某个回答里找到了答案。 注意其中的 e 不是自然常熟而是偏心率，也就是焦距 / 半长轴。 大道至简啊，不过这个公式依旧不好用，想知道任意 t 时刻的r必须先求出参数 ξ 。不过我们可以用二分法来求这个参数，然后就可以求出 t 时刻的r和v了
这个函数输入弧度，也就是参数 ξ 就可以求出对应的 t ，注意这两个量都是全域的，但我们稍后需要限制一下。
//以弧度为参数的时间方程 private float Rad2Time(float rad) { return w * Mathf.Sqrt(w / GM) * (rad - f / w * Mathf.Sin(rad)); }这是一个利用二分递归求弧度的函数，弧度限制在 l 和 r 之间，精度为0.00001。
//用二分法求弧度，输入时间，输出弧度 private float Dichotomy(float l, float r, float time) { if (r - l < 0.00001) return l; float ltime = Rad2Time(l * Mathf.Deg2Rad); float rtime = Rad2Time(r * Mathf.Deg2Rad); float mtime = Rad2Time((l + r) / 2 * Mathf.Deg2Rad); if((ltime - mtime) * (time - mtime) >= 0) return Dichotomy(l, (l + r) / 2, time); return Dichotomy((l + r) / 2, r, time); }OK，我们的准备工作只剩最后一步了，就是时间的初始化。对于参数方程中来说，t = 0就是近日点，而对于我们的行星模拟来说，一开始行星所处的可不是近日点，因此需要初始化一下。
void FixedUpdate() { ... if (t == 0) InitTime(); t += Time.deltaTime; if(Vector3.Cross(r, v).z > 0) SetPosition(t0 + t); else SetPosition(t0 - t); } private void InitTime() { //先计算出周期T T = 2 * Mathf.PI * w * Mathf.Sqrt(w / GM); //初始t0由r和alpha决定 float epsilon = Mathf.Acos((w - r.magnitude) / f); if (alpha < 0) epsilon = -epsilon; t0 = Rad2Time(-epsilon); }当整个程序开始运行时计算出t0，也就是通过位置运动方程来算出参数 ξ ，再算出时间。
接着就是重头戏了！求出任意 t 时刻的r和v。 在这个函数中我们用 beta 和 w 来算出r的绝对值，再用准线法求出x0和y0，注意判断方向，最后求出v，然后转化到三维坐标系中。
private void SetPosition(float time) { //将t调整到周期内 while (time > T / 2) time -= T; while (time < -T / 2) time += T; //求出时间t对应的beta float beta = 0; if (time > 0) beta = Dichotomy(0, Mathf.PI, time); else beta = Dichotomy(-Mathf.PI, 0, time); //由beta求出在XOY坐标系的r和v float absR = w * (1 - f / w * Mathf.Cos(beta)); float x0 = w * w / f - w / f * absR; float y0 = Mathf.Sqrt(Mathf.Abs(absR * absR - (x0 - f) * (x0 - f))); if (beta < 0) y0 = -y0; Vector3 rt = new Vector3(x0 - f, y0, 0); Vector3 vt = new Vector3(-y0 * w / h, x0 * h / w, 0); vt = vt / vt.magnitude * Mathf.Sqrt(GM * (2 / rt.magnitude - 1 / w)); if (Vector3.Cross(rt, vt).z * Vector3.Cross(r, v).z < 0) vt = -vt; //在XOY中绕焦点旋转 Quaternion q = Quaternion.AngleAxis(alpha, new Vector3(0, 0, 1)); rt = q * rt; vt = q * vt; //转换到三维坐标系 rt = rt.x * Xi + rt.y * Yj; vt = vt.x * Xi + vt.y * Yj; this.transform.position = rt + center.transform.position; velocity.position = vt + this.transform.position; }其实在之前的调试中我都是吧Analyse()函数和AddPoints()函数放在Awake()的（不清楚这两个函数的朋友可以去看看我的前一篇博文），这样一旦开始运行 w 和 h 等轨道参数就不会变了，不然物体的位置和速度时刻变化，不知道会出什么幺蛾子。 但是现在我希望可以在程序的运行过程中实时改变物体的位置和速度，产生新的轨道，那就要把它们放回FixUpdate()中了。同时我又不希望它一直重新计算轨道产生资源浪费，就加了一些小限制。
void FixedUpdate() { if(r != this.transform.position - center.position || v != velocity.position - this.transform.position || t ==0) { Analyse(); AddPoints(); t = 0; InitTime(); } t += Time.deltaTime; if(Vector3.Cross(r, v).z > 0) SetPosition(t0 + t); else SetPosition(t0 - t); }显然，程序对r和v的改变只会在SetPosition()的最后进行，因此如果r和v变化了一定是受了外界影响，这时就需要重新计算轨道参数，时间也要重新初始化。
