用.NET模擬天體運動

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

這裡面的所有知識都已經在高中物理提過了，但有兩點需要注意：

核心代碼天體類：

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天才能看完😂），運行效果如下： 

從軌跡可以看出，由於太陽引力的作用，地球會轉著太陽轉，但也同樣由於地球和月球引力的作用，太陽也在以微小的角度在「公轉」。

擴展

如果將太陽質量翻倍( 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狀態：

本文中涉及的所有完整、可運行代碼都已經上傳到我的 github博客，各位可以自行下載：https://github.com/sdcb/blog-data/tree/master/2019/20191214-simulate-planet-movement-using-dotnet

喜歡的朋友 請關注我的微信公眾號：【DotNet騷操作】