画个太阳系
曾经在 Quora 上看到一个问题：海王星和冥王星会相撞吗？
（上图是题主发的图片）
在冥王星被降级之前，各种太阳系示意图都是有画冥王星的。其中画了轨道的，基本上都跟题主这张一个德性。小时候我也曾经产生过它俩会不会相撞的疑问，不过当时我也只是稍微好奇一下，很快就把这件事情抛在脑后了。
那么到底会不会撞呢？不会。因为这两个轨道其实根本不相交！这种示意图都有一定的误导性，让人以为太阳系九（八）大行星的轨道都是共面的。实际上它们的轨道都有点歪，虽然倾角不大，但也绝对不能算共面。这些图画的都是轨道在黄道面上的投影。如果真要画个几何关系正确的图，那得用透视投影。
本文将介绍一种能近似地算出在 3000BC ~ 3000AD 间任意时刻时各大行星的位置和轨道参数的方法（来自 JPL Solar System Dynamics）。我用这个方法用 HTML5 inline SVG 做了一个画太阳系简单示意图（仍然是平面图）的小玩意。代码已发布在 GitHub, 最终效果看这里（传送门）。
JPL 的原文其实已经很完整了，但缺少推导过程。所以本文的主要目的就是补全推导过程。
在这个方法中，我们近似地认为行星的轨道在短时间内是理想的椭圆。
首先，我们在行星（对地球来说是地月系质心，但可以认为地月系质心就是地球质心。由 \( M_{earth}R_{earth}=M_{moon}R_{moon} \) 可以算出地月系质心跟地球质心的距离还不如地球的平均半径大，这在天文的尺度上完全可以忽略）自己的轨道平面上来研究它的位置。
（图片来源：davidcolarusso.com）
如图，椭圆轨道半长轴为 \(a\), 太阳 \(S\) 是椭圆轨道的一个焦点，行星为 \(X\), 行星到太阳的距离为 \(r\), \(P\) 是近日点。
设椭圆轨道离心率为 \(e\), \(\angle{XSP}\) 为 \(\nu\), 称为真近点角，即行星从近日点转过的角度。
以 \(S\) 为原点，\(\overrightarrow{SP}\) 的方向为 \(x\) 轴正方向，建立一个 \(y\) 轴向「上」（过近日点时行星的速度方向）的平面直角坐标系。现在要求行星 \(X\) 在这个坐标系内的坐标。
容易看出，用 \(a\)、\(e\)、\(\nu\) 三个数很容易就能把 \(X\) 的坐标求出来。但是，理想很丰满，现实很骨感，因为 \(\nu\) 并不好算（没法直接通过时间算出来）。
那么，根据时间，能直接算出什么呢？能直接算出行星的平近点角。假设有一个「平均行星」，以这个行星的周期为周期在轨道平面上绕椭圆中心（不是太阳）做匀速圆周运动。它从行星近日点转过的角度就是行星的平近点角，记作 \(M\). 因为行星的周期几乎是恒定的，所以通过时间容易算出这个角。所以，现在问题就变成了：怎么用 \(M\)、\(a\)、\(e\) 三个参数求出行星的坐标？
机智的开普勒给出了一个方法：你不是椭圆轨道吗？我就把你拉成正圆的。如下图：
（图片来源：Wikipedia）
（因为懒得自己画图所以图中字母跟上一张不大一样……将就一下吧）
\(S\) 为太阳，\(P\) 为行星，\(X\) 是把椭圆轨道通过伸缩变换变成半径为 \(a\) （半长轴）的圆后行星的位置。\(Y\) 为「平均行星」，\(M\) 即平近点角。
设椭圆中心为 \(C\), \(X\) 在椭圆长轴上的投影为 \(D\), 近日点为 \(Z\). 定义 \(\angle{XCZ}\) 为偏近点角，记作 \(E\), 作为一个中间跳板来通过平近点角 \(M\) 求行星 \(P\) 的坐标。
具体思路是这样的：先通过 \(M\) 算出 \(E\), 再通过 \(E\) 算出行星 \(P\) 的坐标。
由于椭圆的参数方程
\( \left\{\begin{matrix} x = a\cos{t} \\ y = b\sin{t} \end{matrix}\right. \)
中，参数 \(t\) 的几何意义即是把椭圆通过伸缩变换变成圆后，点 \((x,\;y)\) 与原点的连线跟 \(x\) 轴正半轴的夹角（\(\mathrm{atan2}{(y,\;x)}\)），所以我们可以把偏近点角 \(E\) 看作参数 \(t\). 这样，只要求出了 \(E\), 坐标就直接出来了：
\( \left\{\begin{matrix} x = a\cos{E}-c \\ y = a\sqrt{1-e^2}\sin{E} \end{matrix}\right. \)
其中 \(c\) 是半焦距，减去它是因为我们原来建的坐标系原点是焦点（太阳），而椭圆的参数方程是以椭圆的中心为原点的。由于后面还要进行一次坐标变换（把所有行星放到同一个坐标系中），所以必须以太阳为原点。
由 \(c=ae\) 可知 \( x = a\cos{E}-c = a\cos{E}-ae = a(\cos{E}-e) \). 这样就得到了只与 \(E\)、\(a\)、\(e\) 有关的行星坐标：
\( P\left(a(\cos{E}-e),\;a\sqrt{1-e^2}\sin{E}\right) \)
知道怎么用 \(E\) 求坐标了，现在的问题就是怎么用 \(M\) 求 \(E\). 这就要解「开普勒方程」：
\( M=E-e\sin{E} \)
那么， 这个方程是怎么来的呢？由开普勒第二定律可知，从过近日点开始计时，「平均行星」扫过的面积 \(S_{YCZ}\) 等于行星实际扫过的面积 \(S_{XSZ}\), 所以：
\( S_{YCZ} = S_{XSZ} = S_{XCZ} - S_{\triangle{XSC}} \)
代入它们的值，得：
\( \frac{1}{2}Ma^2 = \frac{1}{2}Ea^2-\frac{1}{2}c\cdot a\sin{E} \)
代入 \(c=ea\), 得：
\( \frac{1}{2}Ma^2 = \frac{1}{2}Ea^2-\frac{1}{2}ea^2\sin{E} \)
消去 \( \frac{1}{2}a^2 \), 得：
\( M=E-e\sin{E} \)
这个推导只适用于图示的情况（\(E<{90}^{\circ}\)），但实际上结论在各种情况下都成立（推导过程类似）。
很明显，这个方程是个超越方程。那么怎么解呢？用牛顿迭代法。
设 \( f(E) = M - E + e\sin{E} \), 可得 \( f'(E) = e\cos{E} - 1 \).
下证 \(f(E)\) 只有一个零点，可以放心大胆地使用牛顿迭代法。
\( \because 0<e<1 \)
\( \therefore 0 \leqslant ecosE \leqslant 1 \)
\( \therefore -1 \leqslant ecosE - 1 \leqslant 0 \)
\( \therefore \) \(f(E)\) 单调递减
\( \therefore \) \(f(E)\) 只有一个零点
现在，我们就算是求出了行星在它轨道平面内的坐标。现在，把视野放宽，求一下各行星在同一坐标系内的坐标。
（图片来源：davidcolarusso.com, 图中 \(\omega\) 的解释写错了，应该是 "argument of perihelion" 才对）
如图，以太阳为原点（设为 \(O\)）， ( J2000.0 平 ) 春分点方向为 \(x\) 轴正方向，建立空间直角坐标系（右手系）。\(xOy\) 平面在 ( J2000.0 平 ) 黄道所在平面内，\(z\) 轴朝上（北）。近日点为 \(P\).
行星自南向北穿过黄道面的点称过升交点，自北向南穿过黄道面的点称为降交点。升交点的（日心）黄经称为升交点黄经，记作 \(\Omega\). 升交点和近日点之间的角距（沿行星运动方向度量）称为近日点参数（也叫近日点角），记作 \(\omega\). 轨道平面和黄道面之间的夹角（锐角）称为轨道倾角，记作 \(I\).
这三个量的作用就是完整地描述行星轨道和黄道面之间的位置关系，以此把行星在自己轨道上的坐标变换成在我们刚建的坐标系中的坐标。
我们已经求出了行星在自己轨道上的坐标，不过是二维的。现在向北拔出一个 \(z\) 轴来，行星坐标 \( (x,\;y) \) 都改写成 \( (x,\;y,\;0) \), 开始进行坐标变换：
① 将 \(x\) 轴和 \(y\) 轴绕 \(z\) 轴在 \(xOy\) 平面内旋转 \(-\omega\), 使 \(x\) 轴正方向指向近日点，变换矩阵为：
\( \begin{bmatrix} \cos{\omega} & -\sin{\omega} & 0 \\ \sin{\omega} & \cos{\omega} & 0 \\ 0 & 0 & 1 \end{bmatrix} \)
② 将 \(y\) 轴和 \(z\) 轴绕 \(x\) 轴在 \(yOz\) 平面内旋转 \(-I\), 使 \(xOy\) 平面与黄道面重合，变换矩阵为：
\( \begin{bmatrix} 1 & 0 & 0 \\ 0 &\cos{I} & -\sin{I} \\ 0 & \sin{I} & \cos{I} \end{bmatrix} \)
③ 将 \(x\) 轴和 \(y\) 轴绕 \(z\) 轴在 \(xOy\) 平面内旋转 \(-\Omega\), 使 \(x\) 轴正方向指向春分点，变换矩阵为：
\( \begin{bmatrix} \cos{\Omega} & -\sin{\Omega} & 0 \\ \sin{\Omega} & \cos{\Omega} & 0 \\ 0 & 0 & 1 \end{bmatrix} \)
三个矩阵乘到一块儿（③×②×①），得到最终的坐标变换矩阵：
\( \begin{bmatrix} cos{\omega} cos{\Omega}-cos{I} sin{\omega} sin{\Omega} & -cos{\Omega} sin{\omega}-cos{I} cos{\omega} sin{\Omega} & sin{I} sin{\Omega} \\ cos{I} cos{\Omega} sin{\omega}+cos{\omega} sin{\Omega} & cos{I} cos{\omega} cos{\Omega}-sin{\omega} sin{\Omega} & -cos{\Omega} sin{I} \\ sin{I} sin{\omega} & cos{\omega} sin{I} & cos{I} \end{bmatrix} \)
至此，所有行星的坐标就都转换成同一坐标系内的坐标了，就可以画图了。
结尾补充一下怎么用通过时间求平近点角 M ( 偷懒直接引用一段 )：
计算时要注意这两点：
- 近日点黄经 \( \varpi = \omega + \Omega \) 就是这么定义的，并不是一个实际上存在的几何上的角。（要不，不在同一面上的俩角怎么能直接相加呢？）平黄经也是用类似的方法定义的，也不是几何上的角。
- \( T_{eph} \) 是指的 TDB 时间（的儒略日），要算出比较精确的数据是要把 UTC 转换成 TDB 的。但对于画示意图这种「超低精度」的计算，直接代入 UTC 时间即可。
最后是两个传送门：
