计算任意时刻椭圆轨道上行星的位置的插值方法

需求

给定椭圆轨道参数（轨道倾角、升交点经度、轨道周期、轨道半长轴、轨道偏心率、初始相位或初始位置），在游戏中实时演算椭圆轨道中行星的位置。

• 可以忽略轨道进动
• 节约计算性能
• 在大轨道半长轴时精度损失相对小

笔者在网络上找了很久，都没有找到一篇完整的适合游戏中计算任意时刻行星在椭圆轨道上位置的文章。于是决定自己研究一下。

搜集到的方法和问题

直接用牛顿力学公式进行模拟的问题

这是最直觉的方法，但问题在于给出合适的初始参数。

物理常数调整问题

在游戏中一般不会直接使用现实中各个物理量的值，例如恒星质量、行星质量、万有引力常数 G 等，游戏开发者往往需要根据需求对现实世界进行缩放。

初始化参量不直观

另外，利用牛顿力学公式，需要根据给定的轨道计算出行星的初始位置和初始速度。一旦参数调整出现一点点小问题，轨道便会偏离设定（例如巨大的轨道进动），更有甚者甚至无法形成稳定的椭圆轨道。
而这和常用描述轨道的参数例如 轨道倾角、升交点经度、轨道周期、轨道半长轴、轨道偏心率 没有直观的关系——我们并不能从初始位置和初始速度简单地推测出这些参数。

实时求解动力学微分方程

实时解微分方程是什么性能黑洞，更别提精度浪费的问题了。

插值方法

问题简化

我们先不考虑轨道倾角和升交点经度，把行星的轨道放到平面上。

根据高中物理书上的知识：行星的轨道呈椭圆形，其围绕的恒星在椭圆的一个焦点上。
那么，把所围绕的恒星放在坐标原点，用极坐标表示行星的椭圆轨道为：
$$\begin{array}{c}r=\frac{p}{1+e\cos \theta }\end{array}$$
这里 $p$ 是椭圆的通径、$e$ 是椭圆的离心率，$p=b^2/a,0<e<1$.
$$p=a(1-e^2)$$
这样一来，行星在任意时刻的位置问题就变成了 求解 $\theta (t)$ 的问题。

求解 $\theta (t)$ (?)

这里复习一下大学物理知识：角动量的定义和角动量守恒（开普勒第二定律）
$$\begin{array}{c}J:=r^2\dot{\theta }=const.\end{array}$$
这里 $\dot{\theta }=\frac{\mathrm{d}\theta }{\mathrm{d}t}$，转化为 $\theta (t)$ 的方程，并设无量纲化时间 $\tau =\frac{t}{T}$ 使方程核心部分无量纲化，得到
$$\frac{p^2}{JT}\frac{{\theta }^{\prime }(\tau )}{(1+e\cos \theta {)}^2}=1$$
显然这个微分方程中很难解出 $\theta (\tau )$，但没关系，我们的目的是获取 $\theta (\tau )$ 的插值，而且 $\theta (\tau )$ 显然是个单调递增的函数，因此也可以去求解他的反函数 $\tau (\theta )$，有
$$\begin{array}{c}{\tau }^{\prime }(\theta )=\frac{p^2}{JT}\frac{1}{(1+e\cos \theta {)}^2}\end{array}$$
i.e.
$$\begin{array}{c}\tau (\theta )=\frac{p^2}{JT}{\int }_0^{\theta }\frac{1}{(1+e\cos \varphi {)}^2}\mathrm{d}\varphi \end{array}$$
由此，我们可以用积分的方法获得 $\theta (\tau )$ 上的坐标点。现在只剩下一个问题：如何得到式子中的常数项。（虽然有物理直觉的同学可以看出：这个常数项只与偏心率 $e$ 有关，而且大概率含一个 $2\pi$）

常数项求解

因为 $\tau =t/T$ 是无量纲化的时间，所以行星运行一周，即 $\theta =2\pi$ 时，$\tau =1$
这个活就交给 Mathematica 了，答案是 $\frac{2\pi }{{\left(1-e^2\right)}^{3/2}}$，所以
$$\begin{array}{c}\tau (\theta )={\tau }_0+\frac{{\left(1-e^2\right)}^{3/2}}{2\pi }{\int }_0^{\theta }\frac{1}{(1+e\cos \varphi {)}^2}\mathrm{d}\varphi \end{array}$$
这里 ${\tau }_0$ 是行星的初始“相位”。从这个公式就可以直接计算 $\tau (\theta )$ 曲线上坐标点了。

Note：开普勒方程

推完一遍发现这个积分有解析解，Mathematica 解出来是

τ ( θ ) = ( 1