标签：style blog http color io 使用 ar 2014 div
这里将以万有引力和势能动能守恒定律为基础,实现行星运动轨迹.然后再假设有两个固定的恒星,让行星在这两个恒星的力场中运动(这是三体问题的一种噢).前面我写过关于混沌曲线的文章:混沌数学及其软件模拟.这类混沌曲线的本质是一个导数方程,即我不知道这条曲线是什么样子,也不知道这条曲线最终通往何处去,我只知道,曲线上的任意一点的切线方向,从而得到它下一点的位置.从而得到整条曲线上的顶点位置.学过物理的人都知道,太阳系中行星的运动轨迹大致是一个椭圆,我们可以知道每个行星的椭圆方程.记得我刚学图形学时,就看过一个程序是太阳行星运动.但它不是以万有引力为基础的,只是让球体绕着固定的椭圆轨迹旋转.
本来我只打算使用万有引力定律,没有考虑势能动能守恒定律.但程序运行后发现,由于没有使用微积分,即使采样间隔时间再小,误差也很大.其实目前的算法误差也不小,呵呵,模拟一下吧,不要太计较.
先帖下基类代码,这代码与混沌数学及其软件模拟中的很相似,即是一个导数方程,用于由当前点计算下一点.
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; // 引力势能
};
这是我写的测试程序,写它是为下一步的三体模拟软件做准备.
标签：style blog http color io 使用 ar 2014 div
原文地址：http://www.cnblogs.com/WhyEngine/p/3977764.html
