這里將以萬有引力和勢能動能守恆定律為基礎,實現行星運動軌跡.然后再假設有兩個固定的恆星,讓行星在這兩個恆星的力場中運動(這是三體問題的一種噢).前面我寫過關於混沌曲線的文章:混沌數學及其軟件模擬.這類混沌曲線的本質是一個導數方程,即我不知道這條曲線是什么樣子,也不知道這條曲線最終通往何處去,我只知道,曲線上的任意一點的切線方向,從而得到它下一點的位置.從而得到整條曲線上的頂點位置.學過物理的人都知道,太陽系中行星的運動軌跡大致是一個橢圓,我們可以知道每個行星的橢圓方程.記得我剛學圖形學時,就看過一個程序是太陽行星運動.但它不是以萬有引力為基礎的,只是讓球體繞着固定的橢圓軌跡旋轉.
本來我只打算使用萬有引力定律,沒有考慮勢能動能守恆定律.但程序運行后發現,由於沒有使用微積分,即使采樣間隔時間再小,誤差也很大.其實目前的算法誤差也不小,呵呵,模擬一下吧,不要太計較.
先帖下基類代碼,這代碼與混沌數學及其軟件模擬中的很相似,即是一個導數方程,用於由當前點計算下一點.
class DifferentiationFunction { public: virtual void Differentiate(float x, float y, float z, float t, float& outX, float& outY, float& outZ) = NULL; virtual float GetExtendT() const = NULL; virtual float GetStartX() const { return 0.0f; } virtual float GetStartY() const { return 0.0f; } virtual float GetStartZ() const { return 0.0f; } };
行星方程:
// 行星的導數曲線
class Planet : public DifferentiationFunction { public: Planet() { m_star_weight = 1.0f; m_planet_x = 5.0f; m_planet_y = 8.0f; m_planet_z = 1.0f; m_planet_weight = 0.1f; m_planet_speed_x = 4.0f; m_planet_speed_y = 0.0f; m_planet_speed_z = 0.0f; m_g = 100.0f; m_ek = 0.5f*m_planet_weight*(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); // 1/2*m*v*v
float r = sqrt(m_planet_x*m_planet_x + m_planet_y*m_planet_y + m_planet_z*m_planet_z); m_ep = -/*0.5f**/m_g*m_star_weight*m_planet_weight/r; m_e = m_ek + m_ep; } void Differentiate(float x, float y, float z, float t, float& outX, float& outY, float& outZ) { t = t*10.0f; float sqd = x*x + y*y + z*z; float d = sqrt(sqd); float a = m_g*m_star_weight/sqd; float ax = -a*x/d; float ay = -a*y/d; float az = -a*z/d; outX = x + m_planet_speed_x*t + 0.5f*ax*t*t; outY = y + m_planet_speed_y*t + 0.5f*ay*t*t; outZ = z + m_planet_speed_z*t + 0.5f*az*t*t; m_planet_speed_x += ax*t; m_planet_speed_y += ay*t; m_planet_speed_z += az*t; float r = sqrt(outX*outX + outY*outY + outZ*outZ); m_ep = -/*0.5f**/m_g*m_star_weight*m_planet_weight/r; m_ek = m_e - m_ep; if (m_ek < 0.0f) { m_ek = 0.0f; } float v = sqrt(2*m_ek/m_planet_weight); float w = sqrt(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); m_planet_speed_x *= v/w; m_planet_speed_y *= v/w; m_planet_speed_z *= v/w; } float GetExtendT() const { return 20.0f; } float GetStartX() const { return m_planet_x; } float GetStartY() const { return m_planet_y; } float GetStartZ() const { return m_planet_z; } public: float m_star_weight; float m_planet_x; float m_planet_y; float m_planet_z; float m_planet_weight; float m_planet_speed_x; float m_planet_speed_y; float m_planet_speed_z; float m_g; // 萬有引力系數
float m_e; float m_ek; // 動能
float m_ep; // 引力勢能
};
行星在這兩個恆星的力場中運動(三體問題)
// 三體問題的導數曲線
class ThreeBody : public DifferentiationFunction { public: ThreeBody() { m_star1_x = -10.0f; m_star1_y = 0.0f; m_star1_z = 0.0f; m_star1_weight = 1.0f; m_star2_x = 10.0f; m_star2_y = 0.0f; m_star2_z = 0.0f; m_star2_weight = 1.0f; m_planet_x = 5.0f; m_planet_y = 5.0f; m_planet_z = 0.1f; m_planet_weight = 0.1f; m_planet_speed_x = 0.0f; m_planet_speed_y = 2.0f; m_planet_speed_z = 0.0f; m_g = 50.0f; m_ek = 0.5f*m_planet_weight*(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); // 1/2*m*v*v
float d1x = m_star1_x - m_planet_x; float d1y = m_star1_y - m_planet_y; float d1z = m_star1_z - m_planet_z; float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z; float d1 = sqrt(sqd1); m_ep1 = -m_g*m_star1_weight*m_planet_weight/d1; float d2x = m_star2_x - m_planet_x; float d2y = m_star2_y - m_planet_y; float d2z = m_star2_z - m_planet_z; float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z; float d2 = sqrt(sqd2); m_ep2 = -m_g*m_star2_weight*m_planet_weight/d2; m_e = m_ek + m_ep1 + m_ep2; } void Differentiate(float x, float y, float z, float t, float& outX, float& outY, float& outZ) { t = t*20.0f; float d1x = m_star1_x - x; float d1y = m_star1_y - y; float d1z = m_star1_z - z; float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z; float d1 = sqrt(sqd1); float d2x = m_star2_x - x; float d2y = m_star2_y - y; float d2z = m_star2_z - z; float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z; float d2 = sqrt(sqd2); float a1 = m_g*m_star1_weight/sqd1; float a1x = a1*d1x/d1; float a1y = a1*d1y/d1; float a1z = a1*d1z/d1; float a2 = m_g*m_star2_weight/sqd2; float a2x = a2*d2x/d2; float a2y = a2*d2y/d2; float a2z = a2*d2z/d2; outX = x + m_planet_speed_x*t + 0.5f*(a1x + a2x)*t*t; outY = y + m_planet_speed_y*t + 0.5f*(a1y + a2y)*t*t; outZ = z + m_planet_speed_z*t + 0.5f*(a1z + a2z)*t*t; m_planet_speed_x += (a1x + a2x)*t; m_planet_speed_y += (a1y + a2y)*t; m_planet_speed_z += (a1z + a2z)*t; { float d1x = m_star1_x - outX; float d1y = m_star1_y - outY; float d1z = m_star1_z - outZ; float sqd1 = d1x*d1x + d1y*d1y + d1z*d1z; float d1 = sqrt(sqd1); m_ep1 = -m_g*m_star1_weight*m_planet_weight/d1; float d2x = m_star2_x - outX; float d2y = m_star2_y - outY; float d2z = m_star2_z - outZ; float sqd2 = d2x*d2x + d2y*d2y + d2z*d2z; float d2 = sqrt(sqd2); m_ep2 = -m_g*m_star2_weight*m_planet_weight/d2; m_ek = m_e - m_ep1 - m_ep2; if (m_ek < 0.0f) { m_ek = 0.0f; } float v = sqrt(2*m_ek/m_planet_weight); float w = sqrt(m_planet_speed_x*m_planet_speed_x + m_planet_speed_y*m_planet_speed_y + m_planet_speed_z*m_planet_speed_z); m_planet_speed_x *= v/w; m_planet_speed_y *= v/w; m_planet_speed_z *= v/w; } } float GetExtendT() const { return 20.0f; } float GetStartX() const { return m_planet_x; } float GetStartY() const { return m_planet_y; } float GetStartZ() const { return m_planet_z; } public: float m_star1_x; float m_star1_y; float m_star1_z; float m_star1_weight; float m_star2_x; float m_star2_y; float m_star2_z; float m_star2_weight; float m_planet_x; float m_planet_y; float m_planet_z; float m_planet_weight; float m_planet_speed_x; float m_planet_speed_y; float m_planet_speed_z; float m_g; // 萬有引力系數
float m_e; float m_ek; // 動能
float m_ep1; // 引力勢能
float m_ep2; // 引力勢能
};
這是我寫的測試程序,寫它是為下一步的三體模擬軟件做准備.下篇文章更精彩:三體運動的程序模擬