用.NET模擬天體運(yùn)動(dòng)

這將是一篇罕見而偏極客的文章。

我上大學(xué)時(shí)就見過一些模擬太陽系等天體運(yùn)動(dòng)的軟件和網(wǎng)站，覺得非常酷炫，比如這個(gè)（http://www.astronoo.com/en/articles/positions-of-the-planets.html）：?

其酷炫之處不僅在于天體運(yùn)動(dòng)軌跡能畫出美妙的弧線，更重要的是其運(yùn)動(dòng)規(guī)律完全由萬有引力定律產(chǎn)生，不需要對其運(yùn)動(dòng)軌跡進(jìn)行編程。所有天體受其它天體的合力，然后按照其加速度運(yùn)行。只需一個(gè)起始坐標(biāo)和起始速度，就能坐下來欣賞畫面。

我從大學(xué)畢業(yè)后就一直對這個(gè)抱有深厚興趣，工作第一年時(shí)我就用?C++做過一版；后來我負(fù)責(zé)公司前端工作，又用?js/?canvas又做了一個(gè)重制版；由于近期爆發(fā)的?.NET“革命”，我近期又用?C#/?.NET再次重制了一版。

需要的數(shù)學(xué)知識(shí)

由于是程序看數(shù)學(xué)知識(shí)，此處我將使用代碼來表示公式。

• 萬有引力，該力?F與兩個(gè)天體的質(zhì)量?m1,?m2成正比，如距離?r的平方成反比，用代碼表示為：?F=m1*m2*G/r^2；

• 牛頓第二定律，加速度?a等于合力?F除以質(zhì)量?m，用代碼表示為：?a=F/m；

• 速度?v與加速度?a以及時(shí)間?t的關(guān)系，用代碼表示為：?v=a*t；

• 距離?s與速度?v以及時(shí)間?t的關(guān)系，用代碼表示為：?s=v*t。

這里面的所有知識(shí)都已經(jīng)在高中物理提過了，但有兩點(diǎn)需要注意：

• 所有的力、加速度、速度以及距離都需要分為?x軸和?y軸兩個(gè)分量；

• 所有的時(shí)間?t實(shí)際上是小段時(shí)間?dt，程序?qū)⒀h(huán)模擬小段時(shí)間累加起來，來模擬天體運(yùn)動(dòng)。

核心代碼

天體類：

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，直接給值為其質(zhì)量的對數(shù)乘?2，另外注意我使用了一個(gè)?PositionTrack的?鏈表來存儲(chǔ)其運(yùn)動(dòng)軌跡。

萬有引力、加速度、速度計(jì)算

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);} }

注意其中有個(gè)?foreach循環(huán)，它將一一計(jì)算每個(gè)天體對某天體的力，然后通過累加的方式求出合力，最后依照合力計(jì)算加速度和速度。如果使用?gmp等高精度計(jì)算庫，該循環(huán)將存在性能熱點(diǎn)，因此可以將這個(gè)?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; }

注意我在繪圖代碼邏輯中做了一些矩陣變換，我把所有邏輯做成了窗體分辨率無關(guān)的，假定屏幕長和寬的較小值為?100，然后左上角坐標(biāo)為?-50,-50，右下角坐標(biāo)為?50,50。

星系模擬

太陽、地球和月亮

這是最容易想到了，地球繞太陽轉(zhuǎn)，月亮繞地球轉(zhuǎn)，創(chuàng)建代碼如下：

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 }); }

注意所有數(shù)據(jù)都沒使用真實(shí)數(shù)字模擬（不然地球繞太陽轉(zhuǎn)一圈需要?365天才能看完????），運(yùn)行效果如下：?

從軌跡可以看出，由于太陽引力的作用，地球會(huì)轉(zhuǎn)著太陽轉(zhuǎn)，但也同樣由于地球和月球引力的作用，太陽也在以微小的角度在“公轉(zhuǎn)”。

擴(kuò)展

如果將太陽質(zhì)量翻倍(?1000->?2000)，會(huì)是何種效果呢？

可見這樣一來，由于引力太大，導(dǎo)致地球速度變快，月亮就被地球“甩”出去了，然后地球軌道也變成了實(shí)實(shí)在在的橢圓。

雙子星

宇宙中存在這樣一種星系，它的兩顆恒星互相圍繞對方轉(zhuǎn)，也可以模擬出來：?

注意兩個(gè)天體在接近時(shí)速度會(huì)變快，遠(yuǎn)離時(shí)速度會(huì)變慢，這是由于萬有引力與距離平方成反比決定的。

擴(kuò)展N星系統(tǒng)

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,};} }

通過精密的數(shù)學(xué)計(jì)算，可以讓任意多的天體組織為系統(tǒng)，如將?3當(dāng)作?N傳入函數(shù)，即可組織為“三星系統(tǒng)”，運(yùn)行效果如下：?

注意，超過?2星的系統(tǒng)都不穩(wěn)定（因此“三星系統(tǒng)”也不穩(wěn)定），轉(zhuǎn)過兩圈之后所有天體由于?double類型的誤差已經(jīng)累積到不可逆轉(zhuǎn)的程度，“三星系統(tǒng)”會(huì)慢慢崩潰解體。

看看四星系統(tǒng)，命運(yùn)也差不多（又比“三星”稍穩(wěn)定，需要等待好幾圈才崩潰）：?

展望與總結(jié)

由于誤差是?double類型的精度限制而累積的，在?C++中我可以使用?gmp、?mpir、?mpfr等高精度計(jì)算庫來模擬計(jì)算，性能也非常高。我之前使用?C++和?mpir/?boost配合，可以讓四星系統(tǒng)穩(wěn)定運(yùn)行長達(dá)?15分鐘不崩潰，還能在我的?WindowsPhone(????)上流暢運(yùn)行。

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

大小數(shù)計(jì)算在天文、地震、天氣、海洋等科研領(lǐng)域有不可取代的作用，我挺希望?.NET能提供一個(gè)高性能、高精度的小數(shù)計(jì)算庫，如?BigFloat。有人會(huì)問?.NET4.0不是提供了?BigInteger嗎？難道不夠？是真不夠！整數(shù)計(jì)算和小數(shù)計(jì)算不完全一樣，性能這一關(guān)就過不去。但在?.NETCore中這個(gè)問題官方似乎沒有太大動(dòng)力去做，我在?github上找到了幾個(gè)相關(guān)?issue，都是?open狀態(tài)：

• https://github.com/dotnet/corefx/issues/17267

• https://github.com/dotnet/corefxlab/issues/2635

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

喜歡的朋友 請關(guān)注我的微信公眾號(hào)：【DotNet騷操作】

總結(jié)

以上是生活随笔為你收集整理的用.NET模拟天体运动的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。