zoukankan      html  css  js  c++  java
  • N体运动的程序模拟

          这依然是与《三体》有关的一篇文章。空间中三个星体在万有引力作用下的运动被称之为三体问题,参见我的上一篇文章:三体运动的程序模拟。而这一节,对三体问题进行了扩展,实现了空间中N个星体在万有引力下的运动模拟。

    程序中使用了两个物理定律:

    (1)万有引力定律

    这是牛顿发现的:任意两个质点有通过连心线方向上的力相互吸引。该引力的大小与它们的质量乘积成正比,与它们距离的平方成反比,与两物体的化学本质或物理状态以及中介物质无关。

    以数学表示为:F=Gfrac{m_1m_2}{r^2}

    其中:

    • F: 两个物体之间的引力
    • G: 万有引力常数
    • {{m}_{1}}: 物体1的质量
    • {{m}_{2}}: 物体2的质量
    • r: 两个物体之间的距离

    依照国际单位制,F的单位为牛顿(N),m1m2的单位为千克(kg),r 的单位为米(m),常数G近似地等于6.67 × 10−11 N m2 kg−2(牛顿米的平方每千克的平方)。当然程序中G不可能设置为这么小的一个数,用户可以自己设置万有引力常数,以实现合理的运动.

    (2)势能与动能守恒定律

    引力位能或引力势能是指物体因为大质量物体的万有引力而具有的位能,其大小与其到大质量的距离有关。

    E_p = -frac{GMm}{r}

    其中G为万有引力常数,M、m分别为两物体质量,r为两者距离。

    动能,简单的说就是指物体因运动而具有的能量。数值上等于(1/2)Mv^2. 动能是能量的一种,它的国际单位制下单位是焦耳(J),简称焦。

    势能与动能守恒即在整个运动过程中,其总能量是不变的,即N体的动能与势能之和为固定值.

         本来我只打算使用万有引力定律,没有考虑势能动能守恒定律.但程序运行后发现,由于没有使用微积分,即使采样间隔时间再小,误差也很大.这让我想起当年学大学物理的时候,大一上学期开高数,然后下学期开了物理.现在对大学物理讲的东西,只记得是用高数中的微积分解物理问题.可那时,我高数学得就是个渣,所以物理也没办法学好.当然就算物理学好了,现在也早忘光了.没有用微积分,这个程序中算法误差不小,模拟一下吧,不要太计较.

          核心代码也发一下,希望有兴趣的人可以对其做出优化:

      1 /****************************************************************
      2 
      3   File name   :  NBodyKernel.h
      4   Author      :  叶峰
      5   Version     :  1.0a
      6   Create Date :  2014/09/19   
      7   Description :  
      8 
      9 *****************************************************************/
     10 
     11 // --------------------------------------------------------------------------------------
     12 
     13 #ifndef _NBodyKernel_H_
     14 #define _NBodyKernel_H_
     15 
     16 // INCLUDES -----------------------------------------------------------------------------
     17 
     18 #include "......WhyYCommon_hWhyMath.h"
     19 #include "......WhyYInterfaceOthersYdD3D9Math.h"
     20 
     21 // --------------------------------------------------------------------------------------
     22 
     23 #define YD_ORBIT_VERTICES_COUNT         1024
     24 
     25 // --------------------------------------------------------------------------------------
     26 
     27 class Star
     28 {
     29 public:
     30     Vector3 vPosition;          // 位置
     31     Vector3 vVelocity;          // 速度
     32     Vector3 vForce;             // 受力
     33     Yreal   fWeight;            // 质量
     34     Yreal   fRadiusScale;       // 缩放
     35     Yreal   fKineticEnergy;     // 动能
     36     Ycolor  color;              // 颜色
     37     Ybool   bOrbitVisible;      // 轨迹可见
     38     Yuint   id;
     39     Vector3 verticesOrbit[YD_ORBIT_VERTICES_COUNT];
     40     Star*   pNext;
     41 
     42     Star()
     43     {
     44         fWeight = 0.0f;
     45         fRadiusScale = 1.0f;  
     46         fKineticEnergy = 0.0f;
     47         color = 0xffffffff;
     48         bOrbitVisible = YD_FALSE;
     49         id = -1;
     50         memset(verticesOrbit, 0, sizeof(verticesOrbit));
     51         pNext = NULL;
     52     }
     53 
     54     // 根据速度计算动能
     55     void CalculateKineticEnergyByVelocity()
     56     {
     57         Yreal sqv = D3DXVec3LengthSq(&vVelocity);
     58         fKineticEnergy = 0.5f*fWeight*sqv;
     59     }
     60 
     61     // 根据动能计算速度
     62     void CalculateVelocityByKineticEnergy()
     63     {
     64         float v = sqrt(2*fKineticEnergy/fWeight);
     65         float w = D3DXVec3Length(&vVelocity);
     66         vVelocity *= v/w;
     67     };
     68 };
     69 
     70 // --------------------------------------------------------------------------------------
     71 
     72 class CNBodyKernel
     73 {
     74 public:
     75     CNBodyKernel();
     76 
     77     ~CNBodyKernel();
     78 
     79     // 更新N体
     80     void        UpdateNBody(Yreal deltaTime);
     81 
     82     // 清空N体
     83     void        Clear();
     84 
     85     // 设置引力系数
     86     void        SetGravitationCoefficient(Yreal c);
     87 
     88     // 获得引力系数
     89     Yreal       GetGravitationCoefficient() const
     90     {
     91         return m_fGravitationCoefficient;
     92     }
     93 
     94     // 返回星体数目
     95     Yuint       GetStarsCount() const
     96     {
     97         return m_starsCount;
     98     }
     99 
    100     // 返回星体链表
    101     const Star* GetStarsList() const
    102     {
    103         return m_starsList;
    104     }
    105 
    106     // 返回星体对象
    107     const Star* GetStar(Yuint index) const;
    108 
    109     // 移除星体
    110     bool        RemoveStar(Yuint index);
    111 
    112     // 添加星体
    113     void        AddStar(
    114         const Vector3& vPosition,
    115         const Vector3& vVelocity,
    116         const Yreal fWeight,
    117         Ycolor  color);
    118 
    119 
    120     // 设置星体的重量
    121     void        SetStarWeight(Yuint index, Yreal weight);
    122 
    123     // 设置星体的轨迹是否可见
    124     void        SetStarOrbitVisible(Yuint index, bool visible);
    125 
    126     // 设置星体的颜色
    127     void        SetStarColor(Yuint index, Ycolor color);
    128 
    129     Yuint       GetCurrentOrbit() const
    130     {
    131         return m_currentOrbit;
    132     }
    133 
    134     const Vector3* GetStarPosition(Yuint index) const;
    135 
    136     Yreal       GetStarWeight(Yuint index) const;
    137 
    138     bool        GetStarOrbitVisible(Yuint index) const;
    139 
    140     Ycolor      GetStarColor(Yuint index) const;
    141 
    142     Yuint       GetStarID(Yuint index) const;
    143 
    144     void        GetSpaceCenter(Vector3& vCenter) const;
    145 
    146     void        GetWeightCenter(Vector3& vCenter) const; 
    147 
    148 private:
    149     // 更新星体的位置,速度,动能
    150     void        UpdateStar(Star& star, Yreal t);
    151 
    152     // 根据动能重新计算速度
    153     void        RecalculateStarVelocity(Star& star);
    154 
    155     // 计算每个星体的当前受力
    156     void        CalculateStarsForce();
    157 
    158     // 计算总势能
    159     void        CalculateAmountOfPotentialEnergy();
    160 
    161     // 计算总动能
    162     void        CalculateAmountOfKineticEnergy();
    163 
    164     // 计算总能量
    165     void        CalculateAmountEnergy()
    166     {
    167         CalculateAmountOfPotentialEnergy();
    168         CalculateAmountOfKineticEnergy();
    169         m_fAmountEnergy = m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;
    170     }
    171 
    172 private:
    173     Yreal m_fGravitationCoefficient;    // 引力系数
    174     Yreal m_fAmountOfPotentialEnergy;   // 当前总势能
    175     Yreal m_fAmountOfKineticEnergy;     // 当前总动能
    176     Yreal m_fAmountEnergy;              // 当前能量m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;
    177 
    178     Star* m_starsList;                  // N体链表
    179     Yuint m_starsCount;                 // N体数目
    180 
    181     Yuint m_currentOrbit;
    182 };
    183 
    184 // --------------------------------------------------------------------------------------
    185 
    186 #endif
      1 /****************************************************************
      2 
      3   File name   :  NBodyKernel.cpp
      4   Author      :  叶峰
      5   Version     :  1.0a
      6   Create Date :  2014/09/19
      7   Description :  
      8 
      9 *****************************************************************/
     10 
     11 // --------------------------------------------------------------------------------------
     12 
     13 #include "......WhyYCommon_hWhyMemoryDefines.h"
     14 #include "NBodyKernel.h"
     15 #include <assert.h>
     16 
     17 // --------------------------------------------------------------------------------------
     18 
     19 #define YD_ESP                          0.00001f
     20 
     21 // --------------------------------------------------------------------------------------
     22 
     23 CNBodyKernel::CNBodyKernel()
     24 {
     25     m_fGravitationCoefficient = 100.0f;
     26     m_fAmountOfPotentialEnergy = 0.0f;
     27     m_fAmountOfKineticEnergy = 0.0f;
     28     m_fAmountEnergy = 0.0f;
     29 
     30     m_starsList = NULL;
     31     m_starsCount = 0;
     32 
     33     m_currentOrbit = 0;
     34 } 
     35 
     36 CNBodyKernel::~CNBodyKernel()
     37 {
     38     Clear();
     39 }
     40 
     41 // 清空N体
     42 void        CNBodyKernel::Clear()
     43 {
     44     Star* pDel = m_starsList;
     45     Star* pNext;
     46     while (pDel)
     47     {
     48         pNext = pDel->pNext;
     49         YD_DELETE(pDel);
     50         pDel = pNext;
     51     }
     52 
     53     m_starsList = NULL;
     54     m_starsCount = 0;
     55 }
     56 
     57 // 获得引力系数
     58 void        CNBodyKernel::SetGravitationCoefficient(Yreal c)
     59 {
     60     m_fGravitationCoefficient = c;
     61 
     62     // 计算总能量
     63     CalculateAmountEnergy();
     64 }
     65 
     66 // 计算每个星体的当前受力
     67 void        CNBodyKernel::CalculateStarsForce()
     68 {
     69     // 清空所有引力
     70     Star* pStar = m_starsList;
     71     while (pStar)
     72     {
     73         pStar->vForce = Vector3(0.0f, 0.0f, 0.0f);
     74         pStar = pStar->pNext;
     75     }
     76 
     77     // 计算引力
     78     Star* pCurrent = m_starsList;
     79     Star* pNext;
     80     while (pCurrent)
     81     {
     82         pNext = pCurrent->pNext;
     83         while (pNext)
     84         {
     85             Vector3 vAB = pNext->vPosition - pCurrent->vPosition;
     86             float disSqAB = D3DXVec3LengthSq(&vAB);
     87             if (disSqAB < YD_ESP)
     88             {
     89                 disSqAB = YD_ESP;
     90             }
     91             float disAB = sqrt(disSqAB);
     92             Vector3 vForceDir = vAB/disAB;
     93             Yreal fForce = m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disSqAB;
     94             Vector3 vForce = vForceDir*fForce;
     95 
     96             pCurrent->vForce += vForce;
     97             pNext->vForce -= vForce;
     98 
     99             pNext = pNext->pNext;
    100         }
    101         pCurrent = pCurrent->pNext;
    102     }
    103 }
    104 
    105 void        CNBodyKernel::UpdateStar(Star& star, Yreal t)
    106 {
    107     // 加速度
    108     const Vector3 acc = star.vForce/star.fWeight;
    109 
    110     star.vPosition.x = star.vPosition.x + star.vVelocity.x*t + 0.5f*acc.x*t*t;
    111     star.vPosition.y = star.vPosition.y + star.vVelocity.y*t + 0.5f*acc.y*t*t;
    112     star.vPosition.z = star.vPosition.z + star.vVelocity.z*t + 0.5f*acc.z*t*t;
    113 
    114     star.vVelocity.x += acc.x*t;
    115     star.vVelocity.y += acc.y*t;
    116     star.vVelocity.z += acc.z*t;
    117 
    118     // 重新计算动能
    119     star.CalculateKineticEnergyByVelocity();
    120 }
    121 
    122 void        CNBodyKernel::UpdateNBody(Yreal deltaTime)
    123 {
    124     // 计算每个星体的当前受力
    125     CalculateStarsForce();
    126 
    127     // 更新星体的位置与速度
    128     Star* pStar = m_starsList;
    129     while (pStar)
    130     {
    131         UpdateStar(*pStar, deltaTime);
    132         pStar = pStar->pNext;
    133     }
    134 
    135     // 能量守恒
    136     CalculateAmountOfKineticEnergy();
    137     CalculateAmountOfPotentialEnergy();
    138 
    139     // 理论上的动能
    140     Yreal fKineticEnergy = m_fAmountEnergy - m_fAmountOfPotentialEnergy;
    141     if (fKineticEnergy < 0.0f)
    142     {
    143         fKineticEnergy = 0.0f;
    144     }
    145 
    146     Yreal r = fKineticEnergy/m_fAmountOfKineticEnergy;
    147     pStar = m_starsList;
    148     while (pStar)
    149     {
    150         if (r > YD_ESP)
    151         {
    152             pStar->fKineticEnergy *= r;
    153         }
    154         else
    155         {
    156             pStar->fKineticEnergy = 1.0f;
    157         }
    158         // 根据动能计算速度
    159         pStar->CalculateVelocityByKineticEnergy();
    160         pStar = pStar->pNext;
    161     }
    162 
    163     // 更新轨迹点
    164     pStar = m_starsList;
    165     while (pStar)
    166     {
    167         pStar->verticesOrbit[m_currentOrbit] = pStar->vPosition;
    168         pStar = pStar->pNext;
    169     }
    170     m_currentOrbit++;
    171     if (m_currentOrbit >= YD_ORBIT_VERTICES_COUNT)
    172     {
    173         m_currentOrbit = 0;
    174     }
    175 }
    176 
    177 // 返回星体对象
    178 const Star* CNBodyKernel::GetStar(Yuint index) const
    179 {
    180     if (index >= m_starsCount)
    181     {
    182         return NULL;
    183     }
    184 
    185     const Star* pStar = m_starsList;
    186     while (index && pStar)
    187     {
    188         index--;
    189         pStar = pStar->pNext;
    190     }
    191 
    192     return pStar;
    193 }
    194 
    195 // 移除星体
    196 bool        CNBodyKernel::RemoveStar(Yuint index)
    197 {
    198     if (index >= m_starsCount)
    199     {
    200         return false;
    201     }
    202 
    203     if (index == 0)
    204     {
    205         Star* pStar = m_starsList->pNext;
    206         YD_DELETE(m_starsList);
    207         m_starsList = pStar;
    208     }
    209     else
    210     {
    211         Star* pStar = m_starsList;
    212         Yuint t = index - 1;
    213         while (t && pStar)
    214         {
    215             t--;
    216             pStar = pStar->pNext;
    217         }
    218         assert(pStar);
    219 
    220         Star* pDel = pStar->pNext;
    221         pStar->pNext = pDel->pNext;
    222         YD_DELETE(pDel);
    223     }
    224 
    225     m_starsCount--;
    226 
    227     Yuint id = 0;
    228     Star* pStar = m_starsList;
    229     while (pStar)
    230     {
    231         pStar->id = id;
    232         id++;
    233         pStar = pStar->pNext;
    234     }
    235 
    236     // 重算能量
    237     CalculateAmountEnergy();
    238 
    239     return true;
    240 }
    241 
    242 // 添加星体
    243 void        CNBodyKernel::AddStar(
    244                     const Vector3& vPosition,
    245                     const Vector3& vVelocity,
    246                     const Yreal fWeight,
    247                     Ycolor  color)
    248 {
    249     Star* pNewStar = YD_NEW(Star);
    250     pNewStar->vPosition = vPosition;
    251     pNewStar->vVelocity = vVelocity;
    252     pNewStar->fWeight = fWeight;
    253     pNewStar->fRadiusScale = yf_pow(fWeight, 1.0f/3.0f);
    254     pNewStar->CalculateKineticEnergyByVelocity();
    255     pNewStar->color = color;
    256     pNewStar->id = m_starsCount;
    257     for (Yuint i = 0; i < YD_ORBIT_VERTICES_COUNT; i++)
    258     {
    259         pNewStar->verticesOrbit[i] = vPosition;
    260     }
    261 
    262     if (!m_starsList)
    263     {
    264         m_starsList = pNewStar;
    265     }
    266     else
    267     {
    268         Star* pStar = m_starsList;
    269         while (pStar->pNext)
    270         {
    271             pStar = pStar->pNext;
    272         }
    273         pStar->pNext = pNewStar;
    274     }
    275 
    276     m_starsCount++;
    277 
    278     // 重算能量
    279     CalculateAmountEnergy();
    280 }
    281 
    282 // 设置星体的重量
    283 void        CNBodyKernel::SetStarWeight(Yuint index, Yreal weight)
    284 {
    285     if (index >= m_starsCount)
    286     {
    287         return;
    288     }
    289 
    290     Star* pStar = m_starsList;
    291     while (index && pStar)
    292     {
    293         index--;
    294         pStar = pStar->pNext;
    295     }
    296 
    297     pStar->fWeight = weight;
    298     pStar->fRadiusScale = yf_pow(weight, 1.0f/3.0f);
    299     pStar->CalculateKineticEnergyByVelocity();
    300 
    301     // 重算能量
    302     CalculateAmountEnergy();
    303 }
    304 
    305 // 设置星体的轨迹是否可见
    306 void        CNBodyKernel::SetStarOrbitVisible(Yuint index, bool visible)
    307 {
    308     if (index >= m_starsCount)
    309     {
    310         return;
    311     }
    312 
    313     Star* pStar = m_starsList;
    314     while (index && pStar)
    315     {
    316         index--;
    317         pStar = pStar->pNext;
    318     }
    319 
    320     pStar->bOrbitVisible = visible;
    321 }
    322 
    323 // 设置星体的颜色
    324 void        CNBodyKernel::SetStarColor(Yuint index, Ycolor color)
    325 {
    326     if (index >= m_starsCount)
    327     {
    328         return;
    329     }
    330 
    331     Star* pStar = m_starsList;
    332     while (index && pStar)
    333     {
    334         index--;
    335         pStar = pStar->pNext;
    336     }
    337 
    338     pStar->color = color;
    339 }
    340 
    341 // 计算总动能
    342 void        CNBodyKernel::CalculateAmountOfKineticEnergy()
    343 {
    344     // 动能
    345     m_fAmountOfKineticEnergy = 0.0f;
    346     Star* pStar = m_starsList;
    347     while (pStar)
    348     {
    349         //pStar->CalculateKineticEnergyByVelocity();
    350         m_fAmountOfKineticEnergy += pStar->fKineticEnergy;
    351         pStar = pStar->pNext;
    352     }
    353 }
    354 
    355 // 计算总势能
    356 void       CNBodyKernel::CalculateAmountOfPotentialEnergy()
    357 {
    358     m_fAmountOfPotentialEnergy = 0.0f;
    359 
    360     Star* pCurrent = m_starsList;
    361     Star* pNext;
    362     while (pCurrent)
    363     {
    364         pNext = pCurrent->pNext;
    365         while (pNext)
    366         {
    367             Vector3 vAB = pCurrent->vPosition - pNext->vPosition;
    368             float disAB = D3DXVec3Length(&vAB);
    369             float fPotentialEnergyAB = -m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disAB;
    370             m_fAmountOfPotentialEnergy += fPotentialEnergyAB;
    371             pNext = pNext->pNext;
    372         }
    373         pCurrent = pCurrent->pNext;
    374     }
    375 }
    376 
    377 const Vector3* CNBodyKernel::GetStarPosition(Yuint index) const
    378 {
    379     if (index >= m_starsCount)
    380     {
    381         return 0;
    382     }
    383 
    384     const Star* pStar = m_starsList;
    385     while (index && pStar)
    386     {
    387         index--;
    388         pStar = pStar->pNext;
    389     }
    390 
    391     assert(pStar);
    392     return &pStar->vPosition;
    393 }
    394 
    395 Yreal       CNBodyKernel::GetStarWeight(Yuint index) const
    396 {
    397     if (index >= m_starsCount)
    398     {
    399         return 0.0f;
    400     }
    401 
    402     const Star* pStar = m_starsList;
    403     while (index && pStar)
    404     {
    405         index--;
    406         pStar = pStar->pNext;
    407     }
    408 
    409     assert(pStar);
    410     return pStar->fWeight;
    411 }
    412 
    413 bool        CNBodyKernel::GetStarOrbitVisible(Yuint index) const
    414 {
    415     if (index >= m_starsCount)
    416     {
    417         return false;
    418     }
    419 
    420     const Star* pStar = m_starsList;
    421     while (index && pStar)
    422     {
    423         index--;
    424         pStar = pStar->pNext;
    425     }
    426 
    427     assert(pStar);
    428     return pStar->bOrbitVisible != YD_FALSE;
    429 }
    430 
    431 Ycolor      CNBodyKernel::GetStarColor(Yuint index) const
    432 {
    433     if (index >= m_starsCount)
    434     {
    435         return 0;
    436     }
    437 
    438     const Star* pStar = m_starsList;
    439     while (index && pStar)
    440     {
    441         index--;
    442         pStar = pStar->pNext;
    443     }
    444 
    445     assert(pStar);
    446     return pStar->color;
    447 }
    448 
    449 Yuint      CNBodyKernel::GetStarID(Yuint index) const
    450 {
    451     if (index >= m_starsCount)
    452     {
    453         return 0;
    454     }
    455 
    456     const Star* pStar = m_starsList;
    457     while (index && pStar)
    458     {
    459         index--;
    460         pStar = pStar->pNext;
    461     }
    462 
    463     assert(pStar);
    464     return pStar->id;
    465 }
    466 
    467 void        CNBodyKernel::GetSpaceCenter(Vector3& vCenter) const
    468 {
    469     vCenter = Vector3(0.0f, 0.0f, 0.0f);
    470     if (!m_starsCount)
    471     {
    472         return;
    473     }
    474 
    475     Star* pStar = m_starsList;
    476     while (pStar)
    477     {
    478         vCenter += pStar->vPosition;
    479         pStar = pStar->pNext;
    480     }
    481 
    482     vCenter /= (Yreal)m_starsCount;
    483 }
    484 
    485 void        CNBodyKernel::GetWeightCenter(Vector3& vCenter) const
    486 {
    487     vCenter = Vector3(0.0f, 0.0f, 0.0f);
    488     if (!m_starsCount)
    489     {
    490         return;
    491     }
    492 
    493     Yreal weights = 0.0f;
    494     Star* pStar = m_starsList;
    495     while (pStar)
    496     {
    497         vCenter += pStar->vPosition*pStar->fWeight;
    498         weights += pStar->fWeight;
    499         pStar = pStar->pNext;
    500     }
    501 
    502     vCenter /= weights;
    503 }
    View Code

          程序启动后,会出现4个球体在运动。不知道大家有没有注意到:Win7的开机动画就是四个球体的运动,所以我这程序的默认也是生成四个球体。

    通过UI控件,可以添加删除星体。没有上限,当然星体越多程序越慢。

    可以显示每一个星体的运动轨迹

    用户可以实时修改任意一个星体的质量

    鼠标右键用于控制视角
    键盘U用于开关UI用户界面.
    通过UI用户界面可以设置三个球体的质量,设置万有引力系数,设置天体运行速度,设置球体的显示大小.

    键盘0~9用于开关第N个球体运动轨迹的显示

    键盘G,用于开关三维网格的显示
    键盘C,用于开关坐标轴的显示
    键盘P,用于暂停
    键盘R,用于重置.

    F11 全屏切换

    软件下载地址:http://files.cnblogs.com/WhyEngine/NBody.zip

  • 相关阅读:
    Python入门--14--字典
    Python入门--13--爬虫一
    Python入门--13--递归
    Python入门--12--函数与变量
    Python入门--11--自定义函数
    Python入门--10--序列
    mysql 删除重复记录
    Java 不可编辑的Map
    mysql left join
    mysql 超过5名学生的课
  • 原文地址:https://www.cnblogs.com/WhyEngine/p/3990394.html
Copyright © 2011-2022 走看看