python太阳代码_用86行Python代码模拟太阳系

Python代码模拟的太阳系，包括了水星(Mercury), 金星(Venus)，地球(Earth)，月球(Moon)，火星(Mars)

上面的动画是我用86行Python代码模拟的一个比较真实的太阳系，用到的基本原理就是N体问题，这是计算天体物理里面的一个重要问题和方法。计算天体物理是用计算机来模拟宇宙中的天体和天文现象，从而对它们进行科学研究。真正用于科研的模拟程序大多是非常复杂的，有着上万行甚至更多代码。这里我写了一个非常简短的代码模拟太阳系行星的运动轨迹，目的只是为了好玩以及展现一下计算天体物理的奇妙之处。

我们知道，行星沿着主要由太阳引力决定的轨道运动。要计算行星的轨道，我们需要计算来自太阳的引力并将其积分以获得速度和位置。我们在高中的时候学过牛顿万有引力定律和牛顿运动定律：

然后通过下面的公式来演化行星的位置和速度：

这里用的是欧拉积分法的一个改进方法--反向欧拉法。这个方法不同于欧拉法的是，它是稳定的，所以行星的轨道是闭合的。这个方法既简单又比较精确性。

上面的三个向量方程可以写成分量形式：

用Python编程的一个好处就是Python可以很优雅地处理向量，不同于C/C++。上面的九个方程式(三个向量方程式，每个具有三个分量)可以很容易地用三行代码实现：

p.r += p.v * dt

acc = -2.959e-4 * p.r / np.sum(p.r**2)**(3./2) # in units of AU/day^2

p.v += acc * dt

现在我们可以演化行星的轨道了，我们还缺的是初始条件，也就是某个时间所有行星的位置和速度。我没有用随意假设的初始条件去做一个艺术性的动画，而是从NASA的JPL实验室的Horizons程序里调用太阳系的数据

from astroquery.jplhorizons import Horizons

obj = Horizons(id=1, location="@sun", epochs=Time("2017-01-01").jd, id_type='id'