N體運動的程序模擬

      這依然是與《三體》有關的一篇文章。空間中三個星體在萬有引力作用下的運動被稱之為三體問題，參見我的上一篇文章：三體運動的程序模擬。而這一節，對三體問題進行了擴展，實現了空間中N個星體在萬有引力下的運動模擬。

程序中使用了兩個物理定律：

（1）萬有引力定律

這是牛頓發現的:任意兩個質點有通過連心線方向上的力相互吸引。該引力的大小與它們的質量乘積成正比，與它們距離的平方成反比，與兩物體的化學本質或物理狀態以及中介物質無關。

以數學表示為：

其中：

• : 兩個物體之間的引力
• : 萬有引力常數
• : 物體1的質量
• : 物體2的質量
• : 兩個物體之間的距離

依照國際單位制，F的單位為牛頓(N)，m₁和m₂的單位為千克(kg)，r 的單位為米(m)，常數G近似地等於6.67 × 10⁻¹¹ N m² kg⁻²（牛頓米的平方每千克的平方）。當然程序中G不可能設置為這么小的一個數,用戶可以自己設置萬有引力常數,以實現合理的運動.

(2)勢能與動能守恆定律

引力位能或引力勢能是指物體因為大質量物體的萬有引力而具有的位能，其大小與其到大質量的距離有關。

其中G為萬有引力常數，M、m分別為兩物體質量，r為兩者距離。

動能，簡單的說就是指物體因運動而具有的能量。數值上等於（1/2）Mv^2. 動能是能量的一種，它的國際單位制下單位是焦耳(J)，簡稱焦。

勢能與動能守恆即在整個運動過程中,其總能量是不變的,即N體的動能與勢能之和為固定值.

     本來我只打算使用萬有引力定律,沒有考慮勢能動能守恆定律.但程序運行后發現,由於沒有使用微積分,即使采樣間隔時間再小,誤差也很大.這讓我想起當年學大學物理的時候,大一上學期開高數,然后下學期開了物理.現在對大學物理講的東西,只記得是用高數中的微積分解物理問題.可那時,我高數學得就是個渣,所以物理也沒辦法學好.當然就算物理學好了,現在也早忘光了.沒有用微積分,這個程序中算法誤差不小,模擬一下吧,不要太計較.

      核心代碼也發一下,希望有興趣的人可以對其做出優化:

/****************************************************************

  File name   :  NBodyKernel.h
  Author      :  葉峰
  Version     :  1.0a
  Create Date :  2014/09/19   
  Description :  

*****************************************************************/

// --------------------------------------------------------------------------------------

#ifndef _NBodyKernel_H_
#define _NBodyKernel_H_

// INCLUDES -----------------------------------------------------------------------------

#include "..\..\..\Why\YCommon_h\WhyMath.h"
#include "..\..\..\Why\YInterface\Others\YdD3D9Math.h"

// --------------------------------------------------------------------------------------

#define YD_ORBIT_VERTICES_COUNT         1024

// --------------------------------------------------------------------------------------

class Star
{
public:
    Vector3 vPosition;          // 位置
    Vector3 vVelocity;          // 速度
    Vector3 vForce;             // 受力
    Yreal   fWeight;            // 質量
    Yreal   fRadiusScale;       // 縮放
    Yreal   fKineticEnergy;     // 動能
    Ycolor  color;              // 顏色
    Ybool   bOrbitVisible;      // 軌跡可見
    Yuint   id;
    Vector3 verticesOrbit[YD_ORBIT_VERTICES_COUNT];
    Star*   pNext;

    Star()
    {
        fWeight = 0.0f;
        fRadiusScale = 1.0f;  
        fKineticEnergy = 0.0f;
        color = 0xffffffff;
        bOrbitVisible = YD_FALSE;
        id = -1;
        memset(verticesOrbit, 0, sizeof(verticesOrbit));
        pNext = NULL;
    }

    // 根據速度計算動能
    void CalculateKineticEnergyByVelocity()
    {
        Yreal sqv = D3DXVec3LengthSq(&vVelocity);
        fKineticEnergy = 0.5f*fWeight*sqv;
    }

    // 根據動能計算速度
    void CalculateVelocityByKineticEnergy()
    {
        float v = sqrt(2*fKineticEnergy/fWeight);
        float w = D3DXVec3Length(&vVelocity);
        vVelocity *= v/w;
    };
};

// --------------------------------------------------------------------------------------

class CNBodyKernel
{
public:
    CNBodyKernel();

    ~CNBodyKernel();

    // 更新N體
    void        UpdateNBody(Yreal deltaTime);

    // 清空N體
    void        Clear();

    // 設置引力系數
    void        SetGravitationCoefficient(Yreal c);

    // 獲得引力系數
    Yreal       GetGravitationCoefficient() const
    {
        return m_fGravitationCoefficient;
    }

    // 返回星體數目
    Yuint       GetStarsCount() const
    {
        return m_starsCount;
    }

    // 返回星體鏈表
    const Star* GetStarsList() const
    {
        return m_starsList;
    }

    // 返回星體對象
    const Star* GetStar(Yuint index) const;

    // 移除星體
    bool        RemoveStar(Yuint index);

    // 添加星體
    void        AddStar(
        const Vector3& vPosition,
        const Vector3& vVelocity,
        const Yreal fWeight,
        Ycolor  color);


    // 設置星體的重量
    void        SetStarWeight(Yuint index, Yreal weight);

    // 設置星體的軌跡是否可見
    void        SetStarOrbitVisible(Yuint index, bool visible);

    // 設置星體的顏色
    void        SetStarColor(Yuint index, Ycolor color);

    Yuint       GetCurrentOrbit() const
    {
        return m_currentOrbit;
    }

    const Vector3* GetStarPosition(Yuint index) const;

    Yreal       GetStarWeight(Yuint index) const;

    bool        GetStarOrbitVisible(Yuint index) const;

    Ycolor      GetStarColor(Yuint index) const;

    Yuint       GetStarID(Yuint index) const;

    void        GetSpaceCenter(Vector3& vCenter) const;

    void        GetWeightCenter(Vector3& vCenter) const; 

private:
    // 更新星體的位置,速度,動能
    void        UpdateStar(Star& star, Yreal t);

    // 根據動能重新計算速度
    void        RecalculateStarVelocity(Star& star);

    // 計算每個星體的當前受力
    void        CalculateStarsForce();

    // 計算總勢能
    void        CalculateAmountOfPotentialEnergy();

    // 計算總動能
    void        CalculateAmountOfKineticEnergy();

    // 計算總能量
    void        CalculateAmountEnergy()
    {
        CalculateAmountOfPotentialEnergy();
        CalculateAmountOfKineticEnergy();
        m_fAmountEnergy = m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;
    }

private:
    Yreal m_fGravitationCoefficient;    // 引力系數
    Yreal m_fAmountOfPotentialEnergy;   // 當前總勢能
    Yreal m_fAmountOfKineticEnergy;     // 當前總動能
    Yreal m_fAmountEnergy;              // 當前能量m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;

    Star* m_starsList;                  // N體鏈表
    Yuint m_starsCount;                 // N體數目

    Yuint m_currentOrbit;
};

// --------------------------------------------------------------------------------------

#endif

/****************************************************************

  File name   :  NBodyKernel.cpp
  Author      :  葉峰
  Version     :  1.0a
  Create Date :  2014/09/19
  Description :  

*****************************************************************/

// --------------------------------------------------------------------------------------

#include "..\..\..\Why\YCommon_h\WhyMemoryDefines.h"
#include "NBodyKernel.h"
#include <assert.h>

// --------------------------------------------------------------------------------------

#define YD_ESP                          0.00001f

// --------------------------------------------------------------------------------------

CNBodyKernel::CNBodyKernel()
{
    m_fGravitationCoefficient = 100.0f;
    m_fAmountOfPotentialEnergy = 0.0f;
    m_fAmountOfKineticEnergy = 0.0f;
    m_fAmountEnergy = 0.0f;

    m_starsList = NULL;
    m_starsCount = 0;

    m_currentOrbit = 0;
} 

CNBodyKernel::~CNBodyKernel()
{
    Clear();
}

// 清空N體
void        CNBodyKernel::Clear()
{
    Star* pDel = m_starsList;
    Star* pNext;
    while (pDel)
    {
        pNext = pDel->pNext;
        YD_DELETE(pDel);
        pDel = pNext;
    }

    m_starsList = NULL;
    m_starsCount = 0;
}

// 獲得引力系數
void        CNBodyKernel::SetGravitationCoefficient(Yreal c)
{
    m_fGravitationCoefficient = c;

    // 計算總能量
    CalculateAmountEnergy();
}

// 計算每個星體的當前受力
void        CNBodyKernel::CalculateStarsForce()
{
    // 清空所有引力
    Star* pStar = m_starsList;
    while (pStar)
    {
        pStar->vForce = Vector3(0.0f, 0.0f, 0.0f);
        pStar = pStar->pNext;
    }

    // 計算引力
    Star* pCurrent = m_starsList;
    Star* pNext;
    while (pCurrent)
    {
        pNext = pCurrent->pNext;
        while (pNext)
        {
            Vector3 vAB = pNext->vPosition - pCurrent->vPosition;
            float disSqAB = D3DXVec3LengthSq(&vAB);
            if (disSqAB < YD_ESP)
            {
                disSqAB = YD_ESP;
            }
            float disAB = sqrt(disSqAB);
            Vector3 vForceDir = vAB/disAB;
            Yreal fForce = m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disSqAB;
            Vector3 vForce = vForceDir*fForce;

            pCurrent->vForce += vForce;
            pNext->vForce -= vForce;

            pNext = pNext->pNext;
        }
        pCurrent = pCurrent->pNext;
    }
}

void        CNBodyKernel::UpdateStar(Star& star, Yreal t)
{
    // 加速度
    const Vector3 acc = star.vForce/star.fWeight;

    star.vPosition.x = star.vPosition.x + star.vVelocity.x*t + 0.5f*acc.x*t*t;
    star.vPosition.y = star.vPosition.y + star.vVelocity.y*t + 0.5f*acc.y*t*t;
    star.vPosition.z = star.vPosition.z + star.vVelocity.z*t + 0.5f*acc.z*t*t;

    star.vVelocity.x += acc.x*t;
    star.vVelocity.y += acc.y*t;
    star.vVelocity.z += acc.z*t;

    // 重新計算動能
    star.CalculateKineticEnergyByVelocity();
}

void        CNBodyKernel::UpdateNBody(Yreal deltaTime)
{
    // 計算每個星體的當前受力
    CalculateStarsForce();

    // 更新星體的位置與速度
    Star* pStar = m_starsList;
    while (pStar)
    {
        UpdateStar(*pStar, deltaTime);
        pStar = pStar->pNext;
    }

    // 能量守恆
    CalculateAmountOfKineticEnergy();
    CalculateAmountOfPotentialEnergy();

    // 理論上的動能
    Yreal fKineticEnergy = m_fAmountEnergy - m_fAmountOfPotentialEnergy;
    if (fKineticEnergy < 0.0f)
    {
        fKineticEnergy = 0.0f;
    }

    Yreal r = fKineticEnergy/m_fAmountOfKineticEnergy;
    pStar = m_starsList;
    while (pStar)
    {
        if (r > YD_ESP)
        {
            pStar->fKineticEnergy *= r;
        }
        else
        {
            pStar->fKineticEnergy = 1.0f;
        }
        // 根據動能計算速度
        pStar->CalculateVelocityByKineticEnergy();
        pStar = pStar->pNext;
    }

    // 更新軌跡點
    pStar = m_starsList;
    while (pStar)
    {
        pStar->verticesOrbit[m_currentOrbit] = pStar->vPosition;
        pStar = pStar->pNext;
    }
    m_currentOrbit++;
    if (m_currentOrbit >= YD_ORBIT_VERTICES_COUNT)
    {
        m_currentOrbit = 0;
    }
}

// 返回星體對象
const Star* CNBodyKernel::GetStar(Yuint index) const
{
    if (index >= m_starsCount)
    {
        return NULL;
    }

    const Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    return pStar;
}

// 移除星體
bool        CNBodyKernel::RemoveStar(Yuint index)
{
    if (index >= m_starsCount)
    {
        return false;
    }

    if (index == 0)
    {
        Star* pStar = m_starsList->pNext;
        YD_DELETE(m_starsList);
        m_starsList = pStar;
    }
    else
    {
        Star* pStar = m_starsList;
        Yuint t = index - 1;
        while (t && pStar)
        {
            t--;
            pStar = pStar->pNext;
        }
        assert(pStar);

        Star* pDel = pStar->pNext;
        pStar->pNext = pDel->pNext;
        YD_DELETE(pDel);
    }

    m_starsCount--;

    Yuint id = 0;
    Star* pStar = m_starsList;
    while (pStar)
    {
        pStar->id = id;
        id++;
        pStar = pStar->pNext;
    }

    // 重算能量
    CalculateAmountEnergy();

    return true;
}

// 添加星體
void        CNBodyKernel::AddStar(
                    const Vector3& vPosition,
                    const Vector3& vVelocity,
                    const Yreal fWeight,
                    Ycolor  color)
{
    Star* pNewStar = YD_NEW(Star);
    pNewStar->vPosition = vPosition;
    pNewStar->vVelocity = vVelocity;
    pNewStar->fWeight = fWeight;
    pNewStar->fRadiusScale = yf_pow(fWeight, 1.0f/3.0f);
    pNewStar->CalculateKineticEnergyByVelocity();
    pNewStar->color = color;
    pNewStar->id = m_starsCount;
    for (Yuint i = 0; i < YD_ORBIT_VERTICES_COUNT; i++)
    {
        pNewStar->verticesOrbit[i] = vPosition;
    }

    if (!m_starsList)
    {
        m_starsList = pNewStar;
    }
    else
    {
        Star* pStar = m_starsList;
        while (pStar->pNext)
        {
            pStar = pStar->pNext;
        }
        pStar->pNext = pNewStar;
    }

    m_starsCount++;

    // 重算能量
    CalculateAmountEnergy();
}

// 設置星體的重量
void        CNBodyKernel::SetStarWeight(Yuint index, Yreal weight)
{
    if (index >= m_starsCount)
    {
        return;
    }

    Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    pStar->fWeight = weight;
    pStar->fRadiusScale = yf_pow(weight, 1.0f/3.0f);
    pStar->CalculateKineticEnergyByVelocity();

    // 重算能量
    CalculateAmountEnergy();
}

// 設置星體的軌跡是否可見
void        CNBodyKernel::SetStarOrbitVisible(Yuint index, bool visible)
{
    if (index >= m_starsCount)
    {
        return;
    }

    Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    pStar->bOrbitVisible = visible;
}

// 設置星體的顏色
void        CNBodyKernel::SetStarColor(Yuint index, Ycolor color)
{
    if (index >= m_starsCount)
    {
        return;
    }

    Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    pStar->color = color;
}

// 計算總動能
void        CNBodyKernel::CalculateAmountOfKineticEnergy()
{
    // 動能
    m_fAmountOfKineticEnergy = 0.0f;
    Star* pStar = m_starsList;
    while (pStar)
    {
        //pStar->CalculateKineticEnergyByVelocity();
        m_fAmountOfKineticEnergy += pStar->fKineticEnergy;
        pStar = pStar->pNext;
    }
}

// 計算總勢能
void       CNBodyKernel::CalculateAmountOfPotentialEnergy()
{
    m_fAmountOfPotentialEnergy = 0.0f;

    Star* pCurrent = m_starsList;
    Star* pNext;
    while (pCurrent)
    {
        pNext = pCurrent->pNext;
        while (pNext)
        {
            Vector3 vAB = pCurrent->vPosition - pNext->vPosition;
            float disAB = D3DXVec3Length(&vAB);
            float fPotentialEnergyAB = -m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disAB;
            m_fAmountOfPotentialEnergy += fPotentialEnergyAB;
            pNext = pNext->pNext;
        }
        pCurrent = pCurrent->pNext;
    }
}

const Vector3* CNBodyKernel::GetStarPosition(Yuint index) const
{
    if (index >= m_starsCount)
    {
        return 0;
    }

    const Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    assert(pStar);
    return &pStar->vPosition;
}

Yreal       CNBodyKernel::GetStarWeight(Yuint index) const
{
    if (index >= m_starsCount)
    {
        return 0.0f;
    }

    const Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    assert(pStar);
    return pStar->fWeight;
}

bool        CNBodyKernel::GetStarOrbitVisible(Yuint index) const
{
    if (index >= m_starsCount)
    {
        return false;
    }

    const Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    assert(pStar);
    return pStar->bOrbitVisible != YD_FALSE;
}

Ycolor      CNBodyKernel::GetStarColor(Yuint index) const
{
    if (index >= m_starsCount)
    {
        return 0;
    }

    const Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    assert(pStar);
    return pStar->color;
}

Yuint      CNBodyKernel::GetStarID(Yuint index) const
{
    if (index >= m_starsCount)
    {
        return 0;
    }

    const Star* pStar = m_starsList;
    while (index && pStar)
    {
        index--;
        pStar = pStar->pNext;
    }

    assert(pStar);
    return pStar->id;
}

void        CNBodyKernel::GetSpaceCenter(Vector3& vCenter) const
{
    vCenter = Vector3(0.0f, 0.0f, 0.0f);
    if (!m_starsCount)
    {
        return;
    }

    Star* pStar = m_starsList;
    while (pStar)
    {
        vCenter += pStar->vPosition;
        pStar = pStar->pNext;
    }

    vCenter /= (Yreal)m_starsCount;
}

void        CNBodyKernel::GetWeightCenter(Vector3& vCenter) const
{
    vCenter = Vector3(0.0f, 0.0f, 0.0f);
    if (!m_starsCount)
    {
        return;
    }

    Yreal weights = 0.0f;
    Star* pStar = m_starsList;
    while (pStar)
    {
        vCenter += pStar->vPosition*pStar->fWeight;
        weights += pStar->fWeight;
        pStar = pStar->pNext;
    }

    vCenter /= weights;
}

      程序啟動后,會出現4個球體在運動。不知道大家有沒有注意到：Win7的開機動畫就是四個球體的運動，所以我這程序的默認也是生成四個球體。

通過UI控件，可以添加刪除星體。沒有上限，當然星體越多程序越慢。

可以顯示每一個星體的運動軌跡

 

用戶可以實時修改任意一個星體的質量

 

鼠標右鍵用於控制視角
鍵盤U用於開關UI用戶界面.
通過UI用戶界面可以設置三個球體的質量,設置萬有引力系數,設置天體運行速度,設置球體的顯示大小.

鍵盤0~9用於開關第N個球體運動軌跡的顯示

鍵盤G,用於開關三維網格的顯示
鍵盤C,用於開關坐標軸的顯示
鍵盤P,用於暫停
鍵盤R,用於重置.

F11 全屏切換

軟件下載地址:http://files.cnblogs.com/WhyEngine/NBody.zip