用.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.m

        var ttm = StepDt * StepDt; // t*t/s1.m

        foreach (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/m

        var 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,

    };

    // Earth

    var earth = new Star

    {

        Px = 0, Py = -41,

        Vx = -5, Vy = 0,

        Mass = 100,

        Color = Color.Blue,

    };

    // Moon

    var moon = new Star

    {

        Px = 0, Py = -45,

        Vx = -10, Vy = 0,

        Mass = 10,

    };

    return new StarSystem(new List<Star> { solar, earth, moon });

}

注意所有数据都没使用真实数字模拟（不然地球绕太阳转一圈需要365天才能看完