用.NET模拟天体运动

这将是一篇罕见而偏极客的文章。

我上大学时就见过一些模拟太阳系等天体运动的软件和网站，觉得非常酷炫，比如这个（http://www.astronoo.com/en/articles/positions-of-the-planets.html）：

其酷炫之处不仅在于天体运动轨迹能画出美妙的弧线，更重要的是其运动规律完全由万有引力定律产生，不需要对其运动轨迹进行编程。所有天体受其它天体的合力，然后按照其加速度运行。只需一个起始坐标和起始速度，就能坐下来欣赏画面。

我从大学毕业后就一直对这个抱有深厚兴趣，工作第一年时我就用 C++做过一版；后来我负责公司前端工作，又用 js/ canvas又做了一个重制版；由于近期爆发的 .NET“革命”，我近期又用 C#/ .NET再次重制了一版。

需要的数学知识

由于是程序看数学知识，此处我将使用代码来表示公式。

万有引力，该力 F与两个天体的质量 m1, m2成正比，如距离 r的平方成反比，用代码表示为： F=m1*m2*G/r^2；
牛顿第二定律，加速度 a等于合力 F除以质量 m，用代码表示为： a=F/m；
速度 v与加速度 a以及时间 t的关系，用代码表示为： v=a*t；
距离 s与速度 v以及时间 t的关系，用代码表示为： s=v*t。

这里面的所有知识都已经在高中物理提过了，但有两点需要注意：

所有的力、加速度、速度以及距离都需要分为 x轴和 y轴两个分量；
所有的时间 t实际上是小段时间 dt，程序将循环模拟小段时间累加起来，来模拟天体运动。

核心代码

天体类：

class Star
{public LinkedList<Vector2> PositionTrack = new LinkedList<SharpDX.Vector2>();public double Px, Py, Vx, Vy;public double Mass;public float Size => (float)Math.Log(Mass) * 2;public Color Color = Color.Black;public void Move(double step){Px += Vx * step;Py += Vy * step;PositionTrack.AddFirst(new Vector2((float)Px, (float)Py));if (PositionTrack.Count > 1000){PositionTrack.RemoveLast();}}
}

注意，我没指定大小 Size，直接给值为其质量的对数乘 2，另外注意我使用了一个 PositionTrack的 链表来存储其运动轨迹。

万有引力、加速度、速度计算

void Step()
{foreach (var s1 in Stars){// star velocity// F = G * m1 * m2 / r^2// F has a direction:double Fdx = 0;double Fdy = 0;const double Gm1 = 100.0f;     // G*s1.mvar ttm = StepDt * StepDt; // t*t/s1.mforeach (var s2 in Stars){if (s1 == s2) continue;var rx = s2.Px - s1.Px;var ry = s2.Py - s1.Py;var rr = rx * rx + ry * ry;var r = Math.Sqrt(rr);var f = Gm1 * s2.Mass / rr;var fdx = f * rx / r;var fdy = f * ry / r;Fdx += fdx;Fdy += fdy;}// Ft = ma    -> a = Ft/m// v  = at    -> v = Ftt/mvar dvx = Fdx * ttm;var dvy = Fdy * ttm;s1.Vx += dvx;s1.Vy += dvy;}foreach (var star in Stars){star.Move(StepDt);}
}

注意其中有个 foreach循环，它将一一计算每个天体对某天体的力，然后通过累加的方式求出合力，最后依照合力计算加速度和速度。如果使用 gmp等高精度计算库，该循环将存在性能热点，因此可以将这个 foreach改成 Parallel.For加 lock的方式修改合力 Fdx和 Fdy，可以提高性能（ C++的代码就是这样写的）。

绘图代码

public void Draw(DeviceContext ctx)
{ctx.Clear(Color.DarkGray);using var solidBrash = new SolidColorBrush(ctx, Color.White);float allHeight = ctx.Size.Height;float allWidth = ctx.Size.Width;float scale = allHeight / 100.0f;foreach (var star in Stars){using var radialBrush = new RadialGradientBrush(ctx, new RadialGradientBrushProperties{Center = Vector2.Zero,RadiusX = 1.0f,RadiusY = 1.0f,}, new SharpDX.Direct2D1.GradientStopCollection(ctx, new[]{new GradientStop{ Color = Color.White, Position = 0f},new GradientStop{ Color = star.Color, Position = 1.0f},}));ctx.Transform =Matrix3x2.Scaling(star.Size) *Matrix3x2.Translation(((float)star.Px + 50) * scale + (allWidth - allHeight) / 2, ((float)star.Py + 50) * scale);ctx.FillEllipse(new Ellipse(Vector2.Zero, 1, 1), radialBrush);ctx.Transform =Matrix3x2.Translation(allHeight / 2 + (allWidth - allHeight) / 2, allHeight / 2);foreach (var line in star.PositionTrack.Zip(star.PositionTrack.Skip(1))){ctx.DrawLine(line.First * scale, line.Second * scale, solidBrash, 1.0f);}}ctx.Transform = Matrix3x2.Identity;
}

注意我在绘图代码逻辑中做了一些矩阵变换，我把所有逻辑做成了窗体分辨率无关的，假定屏幕长和宽的较小值为 100，然后左上角坐标为 -50,-50，右下角坐标为 50,50。

星系模拟

太阳、地球和月亮

这是最容易想到了，地球绕太阳转，月亮绕地球转，创建代码如下：

public static StarSystem CreateSolarEarthMoon()
{var solar = new Star{Px = 0, Py = 0,Vx = 0.6, Vy = 0,Mass = 1000,Color = Color.Red,};// Earthvar earth = new Star{Px = 0, Py = -41,Vx = -5, Vy = 0,Mass = 100,Color = Color.Blue,};// Moonvar moon = new Star{Px = 0, Py = -45,Vx = -10, Vy = 0,Mass = 10,};return new StarSystem(new List<Star> { solar, earth, moon });
}

注意所有数据都没使用真实数字模拟（不然地球绕太阳转一圈需要 365天才能看完????），运行效果如下：

从轨迹可以看出，由于太阳引力的作用，地球会转着太阳转，但也同样由于地球和月球引力的作用，太阳也在以微小的角度在“公转”。

扩展

如果将太阳质量翻倍( 1000-> 2000)，会是何种效果呢？

可见这样一来，由于引力太大，导致地球速度变快，月亮就被地球“甩”出去了，然后地球轨道也变成了实实在在的椭圆。

双子星

宇宙中存在这样一种星系，它的两颗恒星互相围绕对方转，也可以模拟出来：

注意两个天体在接近时速度会变快，远离时速度会变慢，这是由于万有引力与距离平方成反比决定的。

扩展N星系统

static IEnumerable<Star> CreateStars(int N)
{for (var i = 0; i < N; ++i){double angle = 1.0f * i / N * Math.PI * 2;double R = 45;double M = 10000 * 2 / (N * Math.Sqrt(N) * Math.Log(N));double v = 5;double px = R * Math.Sin(angle);double py = R * -Math.Cos(angle);double vx = v * Math.Cos(angle);double vy = v * Math.Sin(angle);yield return new Star{Px = px,Py = py,Vx = vx,Vy = vy,Mass = M,};}
}

通过精密的数学计算，可以让任意多的天体组织为系统，如将 3当作 N传入函数，即可组织为“三星系统”，运行效果如下：

注意，超过 2星的系统都不稳定（因此“三星系统”也不稳定），转过两圈之后所有天体由于 double类型的误差已经累积到不可逆转的程度，“三星系统”会慢慢崩溃解体。

看看四星系统，命运也差不多（又比“三星”稍稳定，需要等待好几圈才崩溃）：

展望与总结

由于误差是 double类型的精度限制而累积的，在 C++中我可以使用 gmp、 mpir、 mpfr等高精度计算库来模拟计算，性能也非常高。我之前使用 C++和 mpir/ boost配合，可以让四星系统稳定运行长达 15分钟不崩溃，还能在我的 WindowsPhone(????)上流畅运行。

之前有人将 mpir移植到了 .NET，但不支持 .NETCore（https://github.com/wezeku/Mpir.NET），有人将 mpfr移植到了 .NET（https://github.com/emphasis87/mpfr.NET/pull/5）， .NETCore可以运行，但有坑爹的 bug，我提了 PR，但作者似乎没心思Merge????。

大小数计算在天文、地震、天气、海洋等科研领域有不可取代的作用，我挺希望 .NET能提供一个高性能、高精度的小数计算库，如 BigFloat。有人会问 .NET4.0不是提供了 BigInteger吗？难道不够？是真不够！整数计算和小数计算不完全一样，性能这一关就过不去。但在 .NETCore中这个问题官方似乎没有太大动力去做，我在 github上找到了几个相关 issue，都是 open状态：