本文主要是来描述一套新参数的应用的，所以公式都是直接给出结论，不会详细给出每一步的数学推导。有兴趣的读者可以自行推导（暴算
序章 开普勒轨道参数简介
众所周知，小质量行星在牛顿引力作用下环绕大质量天体的运动是开普勒轨道运动，而开普勒轨道可以用以下六个参数描述：
1. 倾角 \(i\)
2. 升交点经度 \(\Omega\)
3. 近点辐角 \(\omega\)
4. 半长轴 \(a\)（或近点 \(p\)）
5. 偏心率 \(e\)
6. 平近点角 \(M\)（或真近点角 \(T\)）
这六个参数的自由度等价于状态矢量，也即行星位置 \(\vec r\) 和轨道速度 \(\vec v\)。轨道运动中，行星状态矢量的每个分量都随时间非线性变化，而轨道参数大部分均为常数，只有平近点角 \(M\) 随时间线性变化，十分便于计算给定时刻的行星轨道状态，这也是轨道参数化的意义所在。
但以上六个参数在描述轨道时会遇到很多“奇点”，在某些情况下，轨道状态的微小改变会导致轨道参数的剧烈变化，甚至有些轨道完全无法用轨道参数描述。比如抛物线轨道不存在半长轴 \(a\)；轨道趋近于抛物线时，平近点角 \(M\) 对任何轨道位置的取值都会趋于 \(0\) ，并在抛物线轨道中恒为 \(0\)，无法确定行星位置。
如果简单地将参数 \(a\) 与 \(M\) 使用近点 \(p\) 与真近点角 \(T\) 代替，那么描述近抛物线轨道的奇点会消失，但会对其他轨道产生新的奇点：自由落体轨道中 \(p=0, e=1\)，无法据此计算出半长轴 \(a\)、且真近点角永远为 \(\pi\)，无法确定行星位置；抛物线及双曲线轨道中，在行星极其遥远时，真近点角的微小误差会导致轨道位置的完全不确定。
总之，为了同时消去上述奇点（包括抛物线奇点、自由落体奇点及无限远奇点），需要一套新的参数来描述开普勒轨道，这也是本文的主要讨论内容。
第一章 经典开普勒轨道参数化
首先简介一下，在不考虑奇点的情况下，行星轨道位置和速度如何与经典开普勒参数互相转换。
假定一个标准重力参数（即质量与万有引力常量的积）为 \(\mu\) 的中心天体静止在坐标原点，一个质量可忽略的行星绕中心天体作开普勒运动，位移为 \(\vec r\)，速度为 \(\vec v\)。轨道所在平面的法向可以由行星角动量矢量确定： \[\vec j=\vec r\times\vec v\] 而轨道形状可以由偏心率矢量（即无量纲的拉普拉斯-龙格-楞次矢量）确定： \[\vec e=\frac{\vec v\times\vec j}{\mu}-\hat r\]
以行星偏心率矢量方向为 \(x\) 轴，角动量方向为 \(z\) 轴建立坐标系，就可以在 \(xy\) 平面上以如下极坐标表达式描述行星轨道： \[ r=\frac{j^2/\mu}{1+e\cos{\theta}}\] \[v_r=\frac{\mu e}{j}\sin\theta,\quad v_\theta=\frac{j}{r}\] \[v_x=-\frac{\mu}{j}\sin\theta,\quad v_y=\frac{\mu}{j}(e+\cos\theta)\] 其中 \(r=|\vec r|\) 为行星与中心天体的距离；\((v_x,v_y)=\vec v\) 分别为行星速度的 \(x, y\) 分量；\(\theta\) 为真近点角，也即行星当前所在的极坐标角度；\(e=|\vec e|\) 为偏心率；\(j=|\vec j|\) 为角动量。当 \(e=0\) 时，轨道为圆；当 \( 0 \lt e \lt 1 \) 时，轨道为椭圆；当 \(e=1\) 时，轨道为抛物线；当 \(e\gt 1\) 时，轨道为双曲线。
如果需要计算其他时刻行星的轨道位置，则需要将真近点角 \(\theta\) 依次转换为偏近点角 \(E\) 和平近点角 \(M\)，对 \(e\lt 1\)，有：\[M=E-e\sin E\] \[\sqrt{1-e}\tan\frac{\theta}{2}=\sqrt{1+e}\tan\frac{E}{2}\] \[\theta=\arctan\frac{\sqrt{1-e^2}\sin E}{\cos E-e}\] \[E=\arctan\frac{\sqrt{1-e^2}\sin\theta}{e+\cos\theta}\]
对 \(e\gt 1\)，有： \[M=e\sinh E-E\] \[\sqrt{e-1}\tan\frac{\theta}{2}=\sqrt{e+1}\tanh\frac{E}{2}\] \[\theta=\arctan\frac{\sqrt{e^2-1}\sinh E}{e-\cosh E}\] \[E=\mathrm{arctanh}\frac{\sqrt{e^2-1}\sin\theta}{e+\cos\theta}\]
必须注意，本文所有的形如 \(\arctan\frac{y}{x}\) 的反正切都定义为 \( \arg(x+\mathrm{i}y)\)，即为二元 \(\arctan\) 函数。
平近点角 \(M\) 随时间的变化是线性的：\[\frac{\mathrm{d}M}{\mathrm{d}t}=\frac{\mu^2}{j^3}\left|1-e^2\right|^{3/2}\]
对椭圆轨道（\(e\lt 1\)），平近点角具有 \(2\pi\) 的周期；对双曲线轨道（\(e\gt 1\)），平近点角取值范围无限且没有周期性。将行星的位置 \(\vec r\) 和速度 \(\vec v\) 转换成轨道参数和平近点角后，要计算任意其他时刻行星位置，只需根据以上线性关系计算对应时刻的平近点角值，再依次转换回偏近点角、真近点角，将真近点角代入行星位置和速度的表达式，最后转换参考系即可。
行星轨道中还有一些其他较为常用的参数及公式，虽然在本章中没有使用到，但为了完整性还是列在下面。
半长轴 \(a\)、近点 \(p\)、远点 \(p’\)：
\[a\left(1-e^2\right)=p(1+e)=p'(1-e)=\frac{j^2}{\mu}\]
活力公式，一个很著名的轨道速度与距离的关系：
\[\frac{v^2}{\mu}=\frac{2}{r}-\frac{1}{a}\]
第二章 零角动量奇点与零偏心率奇点
第一章的算法对正常轨道来说已经够用了，不过没有解决序章中提到的奇点问题，从本章开始将逐步讨论可能遇到的各种奇点。
首先就是在第一步计算角动量矢量 \(\vec j\) 时可能遇到的奇点： \(\vec j=0\)。这一奇点出现在 \(\vec v\) 平行于 \(\vec r\) （上抛、下抛）时，或者 \(\vec v=0\)（自由落体）时。
其实这个奇点的处理方法是很朴素的。当 \(\vec j=0\) 时，只需寻找一个任意的垂直于 \(\vec r\) 的方向 \(\hat z\)，并赋予行星一个极小的虚拟角动量 \(\vec j=j\hat z\)，这一角动量应足够小，使得对轨道的影响不超过浮点误差的量级 \(\epsilon\)（对IEEE double类型，这一量级是 \(\epsilon\approx 10^{-16}\)）。为了满足这一要求，需保证 \(j/\sqrt{\mu r}\le \epsilon\) 。计算中可选取满足以上要求的最大的 \(j\) 值，即 \(j=\epsilon\sqrt{\mu r}\)，尽量避免取值过小可能引起的浮点溢出。
实际计算中，由于浮点精度有限，当 \(\vec v\) 与 \(\vec r\) 的平行度达到一定程度时，角动量矢量 \(\vec j\) 虽然不严格为 \(0\)，但其方向会因为浮点误差而不再严格垂直于 \(\vec v\) 与 \(\vec r\)，影响后续 \(\vec e\) 的计算，导致行星完全偏离根据 \(\vec j\) 与 \(\vec e\) 所建立坐标系的 \(xy\) 平面，最终导致根据轨道参数计算的行星方向严重偏离实际方向。因此，当计算出 \(\vec r\times\vec v\) 后，应对其进行正交化以降低数值误差，即实际进行如下计算： \[\vec j_0=\vec r\times\vec v \\ \vec j=\vec j_0-\frac{\vec j_0\cdot\vec r}{\vec r\cdot\vec r}\vec r\] 这样计算的 \(\vec j\) 保证了与 \(\vec r\) 正交，避免了上述问题。
第二个可能的奇点是 \(\vec e=0\)，此时无法确定轨道参考系的 \(x\) 轴。另外与在计算 \(\vec j\) 时所遇问题类似，在 \(\vec e\) 接近但不严格为 \(0\) 时，也有可能因浮点误差而出现其方向不再垂直于 \(\vec j\) 的情况，若将其直接设为 \(x\) 轴，就会导致坐标系不正交，进而产生一系列错误。为了避免这种情况，在确定 \(\vec j\) 的方向之后，应该首先确定一个垂直于 \(\hat j\) 的单位矢量作为参考方向，如升交点 \(\hat \Omega\)。确定升交点后，可以从 \(\vec e\) 计算近点辐角 \(\omega\) ，再由 \(\vec j\) 、 \(\hat \Omega\) 与 \(\omega\) 确定 \(x\) 轴：
\[\omega=\arctan\frac{\vec e\cdot\left(\hat j\times\hat \Omega\right)}{\vec e\cdot\hat \Omega}\\
\hat x=\hat \Omega\cos\omega+\left(\hat j\times\hat \Omega\right)\sin\omega\]
即使因 \(\vec e\) 接近 \(0\) 而 \(\omega\) 偏差较大，以此方法确定的 \(x\) 轴也能保证垂直于 \(\hat j\) 。而在 \(\vec e=0\) 时， \(\omega\) 可取任意值，不影响后续计算的正确性。
第三章 抛物线奇点与无限远奇点
啊……终于到了这一章了……之前写第二章的时候我检查了一下以前写的轨道参数化和反参数化代码，发现了整整五个bug……当然现在已经修复了，本文描述的也都是无bug的算法。不过在我之前自认为完美的代码里找出五个严重bug还是我没想到的……果然写文章有助于全面思考。
好了继续讨论。第一章说到，为了计算行星在其他时刻的位置，需要将真近点角 \(\theta\) 依次转换为偏近点角 \(E\) 和平近点角 \(M\)，并给出了在 \(e\lt 1\) 和 \(e\gt 1\) 两种情况下的算法。然而对 \(e=1\)，也就是抛物线轨道，这些算法完全失效。不仅如此，当 \(e\rightarrow 1\) 时，无论趋近方向如何，偏近点角的计算式中都会出现趋于 \(0\) 的 \(|e^2-1|\) 或者 \(|e-1|\)。这在浮点计算中会产生严重的有效精度损失。随着轨道无限接近抛物线，必须使用无限精度的浮点数存储偏心率才能给出具有有效精度的偏近点角，这是实际计算中无法忍受的。
为了解决上述问题，需要定义一个新的轨道形状参数 \(q\) 代替偏心率： \[q=(e-1)(e+1)=e^2-1\] \[e=\sqrt{1+q}\]
这样，当轨道趋于抛物线时，\(q\) 的取值趋于 \(0\)，故而可以有效利用完整的浮点精度。在实际计算 \(q\) 时，如果 \(e\) 与 \(1\) 相距甚远，则从 \(e\) 计算 \(q\) 没有精度损失，是可行的。但对 \(e\approx 1\)，则不能如此计算，必须寻找计算 \(q\) 的其他方法。
定义： \[\rho=\frac{r}{j^2/\mu}=\frac{1}{1+e\cos\theta}\] \[x=\frac{\vec r\cdot\hat x}{j^2/\mu}=\frac{\cos\theta}{1+e\cos\theta}\] \[y=\frac{\vec r\cdot\hat y}{j^2/\mu}=\frac{\sin\theta}{1+e\cos\theta}\]
可以证明如下三个恒等式成立： \[y=\frac{\vec v\cdot\vec r}{ej}\] \[(e-qx)^2=1+qy^2\] \[q=\frac{1+y^2-2\rho}{x^2}\]
于是 \(q\) 也可以从以上第三个恒等关系中计算得出。需要注意，当 \(\theta\approx\pm\pi/2\) 或 \(e\rightarrow\infty\) 时，\(x\approx 0\)，这一恒等式的精度也会受到影响。实际计算时需要比较该恒等式的相对精度 \(|qx^2|/\left(1+y^2\right)\) 与 \(q\) 的定义式的相对精度 \(|q/e^2|\)，仅当该式能给出更加精确的结果，即 \((ex)^2\gt 1+y^2\) 时，使用该恒等式计算 \(q\) 值。否则仍然使用定义式 \(q=e^2-1\) 计算。
在使用上述恒等式计算 \(q\) 值之前，需要首先计算 \(\rho\)、\(x\) 与 \(y\) 的值。\(\rho\) 与 \(x\) 只能从定义式计算，但 \(y\) 除定义式外还可以从前述第一个恒等式计算。当 \(e\rightarrow 1\) 且 \(\theta\rightarrow\pm\pi\) 时，\(y/\rho\rightarrow 0\)，从定义式计算的 \(y\) 值有效精度损失严重。另一方面，当 \(e\rightarrow 0\) 时，行星轨道为近圆形，\(\vec v\) 与 \(\vec r\) 几乎垂直，恒等式的有效精度也会损失严重。故实际计算中，当 \(e\le 1/2\) 时，可使用定义式计算 \(y\)，而当 \(e\gt 1/2\) 时使用恒等式。这样可以保证 \(y\) 总是拥有足够的精度。
计算得到 \(q\) 之后，已经可以精确地描述趋近于抛物线的轨道了，接下来需要解决平近点角 \(M\) 对抛物线轨道无定义的问题。首先利用以下数学关系：\[\begin{aligned}\sin z=\frac{\exp(\mathrm i z)-\exp(-\mathrm i z)}{2\mathrm i}\quad&\quad\sinh z=\frac{\exp(z)-\exp(-z)}2\\\arctan\frac y x=-\mathrm i\ln\frac{x+\mathrm i y}{\sqrt{x^2+y^2}}\quad&\quad\mathrm{arctanh}\frac y x=\ln\frac{x+y}{\sqrt{x^2-y^2}}\end{aligned}\]
将 \(e\lt 1\) 与 \(e\gt 1\) 两种情形下的平近点角 \(M\) 分别展开为真近点角 \(\theta\) 的函数。当 \(e\lt 1\) 时：
\[\begin{aligned}
M&=E-e\sin E\\
&=\arctan\frac{\sqrt{1-e^2}\sin\theta}{e+\cos\theta}-e\sin\arctan\frac{\sqrt{1-e^2}\sin\theta}{e+\cos\theta}\\
&=-\frac{e\sqrt{1-e^2}\sin\theta}{1+e\cos\theta}-\mathrm i\ln\frac{e+\cos\theta+\mathrm i\sqrt{1-e^2}\sin\theta}{1+e\cos\theta}
\end{aligned}\]
当 \(e\gt 1\) 时：
\[\begin{aligned}
M&=e\sinh E-E\\
&=e\sinh\mathrm{arctanh}\frac{\sqrt{e^2-1}\sin\theta}{e+\cos\theta}-\mathrm{arctanh}\frac{\sqrt{e^2-1}\sin\theta}{e+\cos\theta}\\
&=\frac{e\sqrt{e^2-1}\sin\theta}{1+e\cos\theta}-\ln\frac{e+\cos\theta+\sqrt{e^2-1}\sin\theta}{1+e\cos\theta}
\end{aligned}\]
可见两种情况下 \(M\) 的表达式具有极高的相似性。经过仔细观察，显然可以引入约化平近点角：\[m=\frac{M}{\left|e^2-1\right|^{3/2}}\]
使得两种情形下，\(m\) 具有统一的形式：\[\begin{aligned}m&=\frac{e\sin\theta}{(e^2-1)(1+e\cos\theta)}-\frac{1}{\left(e^2-1\right)^{3/2}}\ln\frac{e+\cos\theta+\sqrt{e^2-1}\sin\theta}{1+e\cos\theta}\\
&=\frac{ey}q-\frac{1}{q^{3/2}}\ln\frac{e(1+e\cos\theta)-(e^2-1)\cos\theta+\sqrt{e^2-1}\sin\theta}{1+e\cos\theta}\\
&=\frac{ey}q-\frac{1}{q^{3/2}}\ln\left(e-qx+y\sqrt{q}\right)\end{aligned}\]
此即 \(m\) 的定义式。现在可以使用约化平近点角 \(m\) 代替平近点角 \(M\) 作为描述行星在轨道上的位置的参数。\(m\) 对任何轨道形状的定义式都是统一的，但实际计算中，对椭圆轨道（\(q\lt 0\)），上式计算中会出现虚数，可以使用 \(\arctan\) 替代 \(\ln\) 来避免虚数计算；对双曲线轨道（\(q\gt 0\)），上式中的 \(e-qx\) 恒为正，可以利用恒等关系使用 \(\sqrt{1+qy^2}\) 代替，避免在 \(q\rightarrow\infty\) 的极端情况下损失精度；此时可以将整个 \(\ln\) 式写成 \(\mathrm{arcsinh}\) 的形式，避免在 \(y\sqrt{q}\rightarrow -\infty\) 及 \(y\sqrt{q}\rightarrow 0\) 的两种情况下分别出现的 \(\infty-\infty\) 与 \(\ln 1\) 造成的精度损失；而轨道接近抛物线时，\(q\rightarrow 0\)，上式将成为不定式 \(\infty-\infty\)，无法直接用来计算 \(m\)。
考虑 \(q\rightarrow 0, e\rightarrow 1\)，即轨道趋于抛物线的情形。此时 \(e-qx\rightarrow 1\gt 0\)，也可利用恒等关系使用 \(\sqrt{1+qy^2}\) 代替。将 \(m\) 的表达式写成仅与 \(y\) 和 \(q\) 有关的形式，并在 \(q=0\) 处进行泰勒展开，可得：
\[\begin{aligned}m&=\frac{y\sqrt{1+q}}q-\frac{1}{q^{3/2}}\ln\left(\sqrt{1+qy^2}+y\sqrt{q}\right)\\
&=y(\frac{y^2+3}{6}-\frac{3y^4+5}{40}q+\frac{5y^6+7}{112}q^2+O(q^3)),\quad(q\rightarrow 0)\end{aligned}\]
可见，约化平近点角 \(m\) 在 \(q=0\) 处是平滑的，去除了平近点角 \(M\) 在抛物线轨道上的奇点。实际计算中，当 \(|q|\) 小于一定值时，就应该将轨道视为近抛物线，使用展开式近似计算，避免损失精度。经过一些误差分析，这一判断条件可以设置为 \(|(2-x)q|\lt 6\epsilon^{1/4}\) 或 \(|(y^2+3)q|\lt 12\epsilon^{1/4}\)。这两者是等价的，因为对抛物线轨道有 \(y^2=1-2x\)。在根据轨道状态计算 \(m\) 时应该使用前者，因为在 \(e\rightarrow 1\) 的椭圆轨道中，远点和近点附近均有 \(y\approx 0\)，但只有近点附近满足 \(e-qx\gt 0\)，可以近似为抛物线处理。使用 \(x\) 作为判据可以很好地区分这两种情形。
综上，约化平近点角的计算可实现为（已被更新方法替代）：
\[m=\left\{\begin{array}{ll}
y(\frac{y^2+3}{6}-\frac{3y^4+5}{40}q+\frac{5y^6+7}{112}q^2), &\mathrm{if}\;|(2-x)q|\lt 6\epsilon^{1/4} \\
\frac{ey}q-\frac{1}{q\sqrt{-q}}\arctan\frac{y\sqrt{-q}}{e-qx}, &\mathrm{else\;if}\;q\lt 0 \\
\frac{ey}q-\frac{1}{q^{3/2}}\mathrm{arcsinh}\left(y\sqrt{q}\right), &\mathrm{else}
\end{array}
\right.
\]
——– 2023.7.17 更新 ——–
以上计算 \(m\) 的方法虽然避免了在轨道接近抛物线时产生显著精度丢失，但由于公式仅仅使用了无穷级数中有限的3项来近似抛物线情形，这导致并不能将结果计算到完整浮点精度。且分支选择需要一个与浮点精度 \(\epsilon\) 相关的判据。这里引入一个新的单变量初等函数 \(\mathrm{kep}\) ，以解决上述问题。其定义为
\[\begin{aligned}
\mathrm{kep}~x&=\frac{\arccos x}{(1-x^2)^{3/2}}-\frac{1}{1-x^2}\\
&=-\frac{\mathrm{arccosh}~x}{(x^2-1)^{3/2}}+\frac{1}{x^2-1}\\
&=\frac{1}{(x+1)^3}\left(\frac{1}{3}+x+3\sum_{i=0}^\infty\frac{i!(1-x)^{i+1}}{(2i+5)!!}\right)
\end{aligned}
\]
函数 \(\mathrm{kep}\) 在定义域 \(x\gt -1\) 内连续光滑，单调递减且任意阶可导，且以上三种定义在其定义域内完全等价。但是，定义1与定义2分别在 \(x\gt 1\) 与 \(x\lt 1\) 时会遇到虚数，且在 \(x\rightarrow 1\) 时成为 \(\infty-\infty\) ，导致计算结果显著失精；而定义3需要进行无穷级数求和，在 \(|x-1|\) 较大时收敛极慢，甚至发散。故此，为了计算 \(\mathrm{kep}~x\) 到完整浮点精度，应分别在 \(x\le\frac{1}{2}\), \(x\ge\frac{3}{2}\), \(|x-1|\le\frac{1}{2}\) 时使用定义1、2、3。在使用定义3计算时，无穷级数应求和至级数项足够小，小于浮点误差为止。若需计算 \(\frac{\mathrm d}{\mathrm dx}\mathrm{kep}~x\)，直接按照 \(x\) 所在区间，计算对应定义的导函数即可。
有了 \(\mathrm{kep}\) 函数， \(m\) 可以写成相当简洁的形式：
\[m=\frac{y}{1+e}+y^3\mathrm{kep}(e-qx)\]
要证明这一点，只需要将上式按照 \(\mathrm{kep}\) 的定义2展开，再利用恒等式 \((e-qx)^2=1+qy^2\) 即可：
\[\begin{aligned}
&\frac{y}{1+e}+y^3\mathrm{kep}(e-qx)\\
=&y\frac{e-1}{e^2-1}+y^3\left(-\frac{\mathrm{arccosh}(e-qx)}{((e-qx)^2-1)^{3/2}}+\frac{1}{(e-qx)^2-1}\right)\\
=&y\frac{e-1}{q}+y^3\left(-\frac{\ln(e-qx+\sqrt{(e-qx)^2-1})}{(qy^2)^{3/2}}+\frac{1}{qy^2}\right)\\
=&y\left(\frac{e-1}{q}+\frac{1}{q}\right)-\frac{y^3}{(y^2)^{3/2}}\left(\frac{\ln(e-qx+\sqrt{qy^2})}{q^{3/2}}\right)\\
=&\frac{ey}{q}-\frac{1}{q^{3/2}}\mathrm{sign}~y\cdot\ln(e-qx+\mathrm{sign}~y\cdot y\sqrt{q})\\
=&\frac{ey}{q}-\frac{1}{q^{3/2}}\ln(e-qx+y\sqrt{q})
\end{aligned}\]
注意到这与 \(m\) 的定义式完全相同。上述过程中 \(\mathrm{sign}~y\) 为 \(y\) 的符号，最后一个等号无论 \(\mathrm{sign}~y\) 为何均成立，这是因为 \((e-qx-y\sqrt q)(e-qx+y\sqrt q)=(e-qx)^2-qy^2=1\)。
虽然新的形式非常简洁，但其使用范围也有限制：如果 \(e-qx\rightarrow -1^+\)，即行星位于椭圆轨道的远点附近，此时 \(y\approx 0\) 而 \(\mathrm{kep}(e-qx)\rightarrow\infty\)，出现 \(0\times\infty\) 型不定式；另外，对 \(e\) 与 \(q\) 非常大的双曲线轨道，\(e-qx\gt 0\) 可能成为 \(\infty-\infty\) 型不定式，应利用恒等关系改为计算 \(\sqrt{1+qy^2}\)。
综上，新的计算 \(m\) 的方法可写为：\[
c=\left\{\begin{array}{ll}
e-qx, &q\lt 0\\\sqrt{1+qy^2}, &q\ge 0\end{array}\right.\quad,\quad
m=\left\{\begin{array}{ll}
\frac{ey}q-\frac{1}{q\sqrt{-q}}\arctan\frac{y\sqrt{-q}}{c}, &c\lt 0\\
\frac{y}{1+e}+y^3\mathrm{kep}~c, &c\ge 0\end{array}\right.\quad.
\]
——– 更新完毕 ——–
使用本章所描述的轨道形状参数 \(q\) 和约化平近点角 \(m\) 替代经典轨道参数偏心率 \(e\) 和平近点角 \(M\)，就能完全消除开普勒轨道的抛物线奇点。顺便一提，由于在以上整个过程中，真近点角 \(\theta\) 只在理论推导中出现，从未被实际计算使用，所以也就回避了行星在无限远时，真近点角的微小误差导致的行星位置的不确定性，也即同时解决了无限远奇点。
现在开普勒轨道中的所有奇点均已被解决，但还剩最后一个问题。考虑到 \( q=e^2-1 \)，当 \(e\rightarrow 0\) ，也即轨道趋近于圆时，\( \mathrm de/\mathrm dq=1/(2e)\rightarrow \infty\)，这意味着从 \(q\) 计算的偏心率 \(e\) 的浮点误差会被无限扩大。实际中如果使用IEEE double存储 \(q\) 值，其相对精度为 \(\epsilon\approx 10^{-16}\)。那么当 \(e\lt\sqrt{\epsilon}\approx 10^{-8}\) 时，\( q=e^2-1 \) 将丢失 \(e\) 的全部精度，导致轨道计算的精度不能达到 \(\epsilon\)，只能达到 \(\sqrt{\epsilon}\)。为了避免对圆轨道损失精度，可以定义修正的形状参数 \(q_0\)： \[q_0=\frac{q}{e+1}\]
显然 \(q_0\) 满足 \(q_0=e-1\)，但为了避免在 \(e\approx 1\) 时损失精度，不能用 \(e-1\) 来计算 \(q_0\)。实际计算出 \(q\) 与 \(m\) 后，可以保存 \(q_0\) 而非 \(q\) 作为轨道形状参数；在需要使用 \(q\) 时，利用 \(q=q_0(q_0+2)\) 计算出 \(q\) 即可。由于 \(e=q_0+1\)，\(\mathrm de/\mathrm dq_0=1\)，故从 \(q_0\) 计算偏心率就避免了在 \(e\rightarrow 0\) 时，从 \(q\) 计算偏心率遇到的浮点精度损失。
约化平近点角 \(m\) 随时间的变化也是线性的：\[\frac{\mathrm{d}m}{\mathrm{d}t}=\frac{\mu^2}{j^3}\]
对椭圆轨道（\(q\lt 0\)），约化平近点角具有 \(2\pi/(-q)^{3/2}\) 的周期；对抛物线或双曲线轨道（\(q\ge 0\)），约化平近点角取值范围无限且没有周期性。将行星的位置 \(\vec r\) 和速度 \(\vec v\) 转换成轨道参数和约化平近点角后，要计算任意其他时刻的行星位置，只需根据以上线性关系计算对应时刻的约化平近点角值，再从约化平近点角反解出 \(y, x\)，继而由此计算出 \(\vec r\) 与 \(\vec v\) 即可。
间章 小结
在继续讨论之前，首先将上面讲过的新一套轨道参数和它们的计算方法重新完整地描述一遍。
假定一个标准重力参数（即质量与万有引力常量的积）为 \(\mu\) 的中心天体静止在坐标原点，一个质量可忽略的行星绕中心天体作开普勒运动，位移为 \(\vec r\)，速度为 \(\vec v\)。
则可计算角动量矢量 \(\vec j\)：\[\vec j_0=\vec r\times\vec v \\ \vec j=\vec j_0-\frac{\vec j_0\cdot\vec r}{\vec r\cdot\vec r}\vec r\]
如果上式所得 \(\vec j=0\)，则需任选一个垂直于 \(\vec r\) 的方向 \(\hat z\)，重新令 \(\vec j=\epsilon\hat z\sqrt{\mu r}\) 。其中 \(\epsilon\) 为浮点相对误差，对 IEEE double，\(\epsilon\approx 10^{-16}\)。保证角动量 \(j=|\vec j|\gt 0\)。在确定 \(\vec j\) 的方向后，需要选取一个垂直于 \(\vec j\) 的单位矢量作为参考方向，此处可选择升交点 \(\hat \Omega\)。
继续计算偏心率矢量 \(\vec e\) ：\[\vec e=\frac{\vec v\times\vec j}{\mu}-\hat r\]
以 \(\vec j\) 方向为 \(z\) 轴，\(\vec e\) 方向为 \(x\) 轴建立坐标系 \(xyz\)。为了避免 \(\vec e\) 接近 \(0\) 导致 \(x\) 轴方向误差过大，可按下式确定 \(x\) 轴方向：
\[\omega=\arctan\frac{\vec e\cdot\left(\hat j\times\hat \Omega\right)}{\vec e\cdot\hat \Omega}\\
\hat x=\hat \Omega\cos\omega+\left(\hat j\times\hat \Omega\right)\sin\omega\]
其中 \(\omega\) 为近点辐角。确定 \(xz\) 方向后，\(y\) 轴由右手法则确定。该坐标系可以使用轨道倾角 \(i\)、升交点经度 \(\Omega\) 与近点辐角 \(\omega\) 三个参数来描述，其 \(z\) 轴确定了轨道的法向，\(x\) 轴确定了轨道近点的方向。
继续计算 \[\rho=\frac{r}{j^2/\mu},\quad x=\frac{\vec r\cdot\hat x}{j^2/\mu}\] \[y=\left\{\begin{array}{ll}\frac{\vec r\cdot\hat y}{j^2/\mu}, & e\le 1/2\\
\frac{\vec v\cdot\vec r}{ej}, & e \gt 1/2
\end{array}
\right.
\]
计算轨道形状参数 \(q\)：
\[q=\left\{\begin{array}{ll}
e^2-1, & (ex)^2\le 1+y^2\\
\frac{1+y^2-2\rho}{x^2}, & (ex)^2\gt 1+y^2
\end{array}
\right.
\]
计算修正形状参数 \(q_0\)：
\[q_0=\frac{q}{e+1}\]
计算约化平近点角 \(m\)：
\[
c=\left\{\begin{array}{ll}
e-qx, &q\lt 0\\\sqrt{1+qy^2}, &q\ge 0\end{array}\right.\quad,\quad
m=\left\{\begin{array}{ll}
\frac{ey}q-\frac{1}{q\sqrt{-q}}\arctan\frac{y\sqrt{-q}}{c}, &c\lt 0\\
\frac{y}{1+e}+y^3\mathrm{kep}~c, &c\ge 0\end{array}\right.\quad.
\]
其中函数 \(\mathrm{kep}\) 的定义与计算方法已在本文中给出；\(\arctan\frac{y}{x}\) 为二元 \(\arctan\) 函数，定义为 \( \arg(x+\mathrm{i}y)\)。
在行星运动中，约化平近点角 \(m\) 满足 \[\frac{\mathrm{d}m}{\mathrm{d}t}=\frac{\mu^2}{j^3}\]
对椭圆轨道（\(q\lt 0\)），约化平近点角具有 \(2\pi/(-q)^{3/2}\) 的周期；对抛物线或双曲线轨道（\(q\ge 0\)），约化平近点角取值范围无限且没有周期性。
最终，开普勒轨道可用如下六个参数描述：
1. 倾角 \(i\)
2. 升交点经度 \(\Omega\)
3. 近点辐角 \(\omega\)
4. 角动量 \(j\)
5. 修正形状参数 \(q_0\)
6. 约化平近点角 \(m\)
这六个参数可以完备且无奇点地以完整浮点精度描述一切开普勒轨道。解决了经典开普勒轨道参数无法描述抛物线轨道；无法描述上抛、下抛、自由落体轨道；描述接近抛物线的椭圆与双曲线轨道时丢失精度；描述行星在双曲线轨道中趋于无穷远时的状态时丢失精度等问题。
啊……这一篇东西又已经很长了……还是得分成两个部分。那第一部分就到这里为止吧。从这六个轨道参数反解行星位置的具体实现，以及对轨道绘制的一些计算，放到第二部分。
第二部分链接（咕咕咕）
