用.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天才能看完😂），運行效果如下：
從軌跡可以看出，由於太陽引力的作用，地球會轉着太陽轉，但也同樣由於地球和月球引力的作用，太陽也在以微小的角度在“公轉”。
擴展
如果將太陽質量翻倍(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分鍾不崩潰，還能在我的Windows Phone(😂)上流暢運行。
之前有人將mpir移植到了.NET，但不支持.NET Core（https://github.com/wezeku/Mpir.NET），有人將mpfr移植到了.NET（https://github.com/emphasis87/mpfr.NET/pull/5），.NET Core可以運行，但有坑爹的bug，我提了PR，但作者似乎沒心思Merge😂。
大小數計算在天文、地震、天氣、海洋等科研領域有不可取代的作用，我挺希望.NET能提供一個高性能、高精度的小數計算庫，如BigFloat。有人會問.NET 4.0不是提供了BigInteger嗎？難道不夠？是真不夠！整數計算和小數計算不完全一樣，性能這一關就過不去。但在.NET Core中這個問題官方似乎沒有太大動力去做，我在github上找到了幾個相關issue，都是open狀態：
本文中涉及的所有完整、可運行代碼都已經上傳到我的github博客，各位可以自行下載：https://github.com/sdcb/blog-data/tree/master/2019/20191214-simulate-planet-movement-using-dotnet
喜歡的朋友 請關注我的微信公眾號：【DotNet騷操作】
