三體三體[代碼開源]

本文轉載自查看原文 2015-02-25 12:51 10354 物理/ 開源

假期就這么結束了,不甘心啊!放假前還打算利用這段時間看看高數,學學微積分,想想如何實現精確的三體模擬.結果整個假期里,書都沒有翻一下.罷了,發布下我寫的使用勢能守衡對三體進行近似模擬的代碼,希望能夠拋磚引玉.歡迎大家進來討論,修改,優化.

空間中三個星體,在萬有引力作用下的運動被稱為三體運動.它的輸入數據是:三個物體的質量,位置和速度.在萬有引力的作用下,三個物體的位置和速度會時時發生變化.需要計算的是它們在某一時刻的位置與速度.

下面三個GIF動畫圖像為三體的錄屏:

代碼是C++的,首先要定義一個星體的結構體,包含星體的質量,位置和速度.還有一個"半徑縮放"它與星體的質量有關,只是表示星體的體積,用於圖形的渲染.

struct Body
{
    Vec3 vPosition;     // 位置
    Vec3 vVelocity;     // 速度
    float fWeight;      // 重量
    float fRadiusScale; // 半徑縮放

    Body()
    {
        fWeight = 0.0f;
        fRadiusScale = 0.0f;
    }

    // 計算動能
    float        CalculateKineticEnergy() const
    {
        float sqv = vVelocity.MagnitudeSquared();
        return 0.5f*fWeight*sqv;
    }

    // 設置動能,即通過動能計算速度
    void        SetKineticEnergy(float energy)
    {
        float v = sqrt(2*energy/fWeight);
        float w = vVelocity.Magnitude();
        if (w < FLT_EPSILON)
        {
            vVelocity.x = v;
        }
        else
        {
            vVelocity *= v/w;
        }
    }

    // 設置質量
    void        SetWeight(float weight)
    {
        fWeight = weight;
        fRadiusScale = powf(weight, 1.0f/3.0f);
    }
};

然后是定義一個三體世界的基類,這是一個抽象接口類.它包含三個星體對象和一個引力系數.如果你有自己的算法,可以繼承該類實現它.

  1 class IThreeBody
  2 {
  3 public:
  4 
  5     IThreeBody()
  6     {
  7         m_fGravitationCoefficient = 100.0f;
  8     }
  9 
 10     // 重置
 11     virtual void        Reset()
 12     {
 13         RandomBody(m_bodyA);
 14         RandomBody(m_bodyB);
 15         RandomBody(m_bodyC);
 16     }
 17 
 18     // 按時更新
 19     virtual void        Update(float deltaTime) = NULL;
 20 
 21     // 引力系數
 22     virtual void        SetGravitationCoefficient(float c)
 23     {
 24         m_fGravitationCoefficient = c;
 25     }
 26 
 27     // 星體A的重量
 28     virtual void        SetBodyAWeight(float weight)
 29     {
 30         m_bodyA.SetWeight(weight);
 31     }
 32 
 33     // 星體B的重量
 34     virtual void        SetBodyBWeight(float weight)
 35     {
 36         m_bodyB.SetWeight(weight);
 37     }
 38 
 39     // 星體C的重量
 40     virtual void        SetBodyCWeight(float weight)
 41     {
 42         m_bodyC.SetWeight(weight);
 43     }
 44 
 45     float               GetGravitationCoefficient() const
 46     {
 47         return m_fGravitationCoefficient;
 48     }
 49 
 50     const Body&         GetBodyA() const
 51     {
 52         return m_bodyA;
 53     }
 54 
 55     const Body&         GetBodyB() const
 56     {
 57         return m_bodyB;
 58     }
 59 
 60     const Body&         GetBodyC() const
 61     {
 62         return m_bodyC;
 63     }
 64 
 65 protected:
 66 
 67     // 重置星體
 68     void        RandomBody(Body& body)
 69     {
 70         body.vPosition.x = rand2(-YD_THREE_BODY_EXTEND, YD_THREE_BODY_EXTEND);
 71         body.vPosition.y = rand2(-YD_THREE_BODY_EXTEND, YD_THREE_BODY_EXTEND);
 72         body.vPosition.z = rand2(-YD_THREE_BODY_EXTEND, YD_THREE_BODY_EXTEND);
 73 
 74         body.vVelocity.x = rand2(-1, 1);
 75         body.vVelocity.y = rand2(-1, 1);
 76         body.vVelocity.z = rand2(-1, 1);
 77         body.vVelocity.Normalize();
 78         body.vVelocity *= rand2(YD_THREE_BODY_VELOCITY*0.1f, YD_THREE_BODY_VELOCITY*0.5f);
 79 
 80         body.fWeight = rand2(1.0f, 8.0f);
 81         body.fRadiusScale = powf(body.fWeight, 1.0f/3.0f);
 82     }
 83 
 84     // 計算某一星體的加速度
 85     void        CalculateBodyAcceleration(const Body& body, const Body& s1, const Body& s2, Vec3& acc)
 86     {
 87         Vec3 v1 = s1.vPosition - body.vPosition;
 88         float sqd1 = v1.MagnitudeSquared();
 89         if (sqd1 < YD_ESP)
 90         {
 91             sqd1 = YD_ESP;
 92         }
 93         float d1 = sqrt(sqd1);
 94         float a1 = m_fGravitationCoefficient*s1.fWeight/sqd1;
 95 
 96         Vec3 v2 = s2.vPosition - body.vPosition;
 97         float sqd2 = v2.MagnitudeSquared();
 98         if (sqd2 < YD_ESP)
 99         {
100             sqd2 = YD_ESP;
101         }
102         float d2 = sqrt(sqd2);
103         float a2 = m_fGravitationCoefficient*s2.fWeight/sqd2;
104 
105         float a1x = a1*v1.x/d1;
106         float a1y = a1*v1.y/d1;
107         float a1z = a1*v1.z/d1;
108 
109         float a2x = a2*v2.x/d2;
110         float a2y = a2*v2.y/d2;
111         float a2z = a2*v2.z/d2;
112 
113         acc.x = a1x + a2x;
114         acc.y = a1y + a2y;
115         acc.z = a1z + a2z;
116     }
117 
118 protected:
119     float m_fGravitationCoefficient; // 引力系數
120 
121     Body m_bodyA;
122     Body m_bodyB;
123     Body m_bodyC;
124 };

接下來就是使用動能勢能守衡對三體進行近似模擬的代碼,即我實現的三體模擬,定義一個類MyThreeBody它是IThreeBody的子類:

 1 class MyThreeBody : public IThreeBody
 2 {
 3 public:
 4     MyThreeBody()
 5     {
 6         m_fPotentialEnergyAB = 0.0f;
 7         m_fPotentialEnergyBC = 0.0f;
 8         m_fPotentialEnergyAC = 0.0f;
 9         m_fAmountOfEnergy = 0.0f;   
10     }
11 
12     // 重置
13     void        Reset();
14 
15     // 按時更新
16     void        Update(float deltaTime);
17 
18     // 引力系數
19     void        SetGravitationCoefficient(float c);
20 
21     // 星體A的重量
22     void        SetBodyAWeight(float weight);
23 
24     // 星體B的重量
25     void        SetBodyBWeight(float weight);
26 
27     // 星體C的重量
28     void        SetBodyCWeight(float weight);
29 
30 protected:
31 
32     // 更新星體的位置,速度
33     void        UpdateBody(Body& body, float t, const Vec3& acc)
34     {
35         body.vPosition.x = body.vPosition.x + body.vVelocity.x*t + 0.5f*acc.x*t*t;
36         body.vPosition.y = body.vPosition.y + body.vVelocity.y*t + 0.5f*acc.y*t*t;
37         body.vPosition.z = body.vPosition.z + body.vVelocity.z*t + 0.5f*acc.z*t*t;
38 
39         body.vVelocity.x += acc.x*t;
40         body.vVelocity.y += acc.y*t;
41         body.vVelocity.z += acc.z*t;
42     }
43 
44     // 計算勢能
45     void        CalculatePotentialEnergy()
46     {
47         Vec3 vAB = m_bodyA.vPosition - m_bodyB.vPosition;
48         float disAB = vAB.Magnitude();
49         if (disAB < YD_ESP)
50         {
51             disAB = YD_ESP;
52         }
53         m_fPotentialEnergyAB = -m_fGravitationCoefficient*m_bodyA.fWeight*m_bodyB.fWeight/disAB;
54 
55         Vec3 vBC = m_bodyC.vPosition - m_bodyB.vPosition;
56         float disBC = vBC.Magnitude();
57         if (disBC < YD_ESP)
58         {
59             disBC = YD_ESP;
60         }
61         m_fPotentialEnergyBC = -m_fGravitationCoefficient*m_bodyC.fWeight*m_bodyB.fWeight/disBC;
62 
63         Vec3 vAC = m_bodyA.vPosition - m_bodyC.vPosition;
64         float disAC = vAC.Magnitude();
65         if (disAC < YD_ESP)
66         {
67             disAC = YD_ESP;
68         }
69         m_fPotentialEnergyAC = -m_fGravitationCoefficient*m_bodyA.fWeight*m_bodyC.fWeight/disAC;
70     }
71 
72     // 計算總能量
73     void        CalculateAmountOfEnergy()
74     {
75         // 動能
76         float fKineticEnergyA = m_bodyA.CalculateKineticEnergy();
77         float fKineticEnergyB = m_bodyB.CalculateKineticEnergy();
78         float fKineticEnergyC = m_bodyC.CalculateKineticEnergy();
79 
80         // 勢能
81         CalculatePotentialEnergy();
82 
83         // 總能量
84         m_fAmountOfEnergy = (m_fPotentialEnergyAB + m_fPotentialEnergyAC + m_fPotentialEnergyBC) + 
85             fKineticEnergyA + fKineticEnergyB + fKineticEnergyC;
86     }
87 
88 private:
89     float m_fPotentialEnergyAB; // 勢能AB
90     float m_fPotentialEnergyBC; // 勢能BC
91     float m_fPotentialEnergyAC; // 勢能AC
92     float m_fAmountOfEnergy;    // 三體總能量
93 };

MyThreeBody需要實現一個關鍵的函數Update(float deltaTime),它用於對三體的更新,輸入一個間隔時間,計算出當前三個物體的位置與速度:

 1 void        MyThreeBody::Update(float deltaTime)
 2 {
 3     Vec3 accA, accB, accC;
 4     CalculateBodyAcceleration(m_bodyA, m_bodyB, m_bodyC, accA);
 5     CalculateBodyAcceleration(m_bodyB, m_bodyC, m_bodyA, accB);
 6     CalculateBodyAcceleration(m_bodyC, m_bodyA, m_bodyB, accC);
 7 
 8     UpdateBody(m_bodyA, deltaTime, accA);
 9     UpdateBody(m_bodyB, deltaTime, accB);
10     UpdateBody(m_bodyC, deltaTime, accC);
11 
12     // 動能
13     float fKineticEnergyA = m_bodyA.CalculateKineticEnergy();
14     float fKineticEnergyB = m_bodyB.CalculateKineticEnergy();
15     float fKineticEnergyC = m_bodyC.CalculateKineticEnergy();
16 
17     // 勢能
18     CalculatePotentialEnergy();
19 
20     // 能量守恆
21     float fKineticEnergy = m_fAmountOfEnergy - (m_fPotentialEnergyAB + m_fPotentialEnergyAC + m_fPotentialEnergyBC);
22     if (fKineticEnergy < 0.0f)
23     {
24         fKineticEnergy = 0.0f;
25     }
26 
27     float fKineticEnergy2 = fKineticEnergyA + fKineticEnergyB + fKineticEnergyC;
28     float r = fKineticEnergy/fKineticEnergy2;
29     if (r > YD_ESP)
30     {
31         fKineticEnergyA *= r;
32         fKineticEnergyB *= r;
33         fKineticEnergyC *= r;
34     }
35     else
36     {
37         fKineticEnergyA = 1.0f;
38         fKineticEnergyB = 1.0f;
39         fKineticEnergyC = 1.0f;
40     }
41 
42     m_bodyA.SetKineticEnergy(fKineticEnergyA);
43     m_bodyB.SetKineticEnergy(fKineticEnergyB);
44     m_bodyC.SetKineticEnergy(fKineticEnergyC);
45 }