计算行星之间的吸引力方向

Question

我正在尝试使用C在opengl中建立太阳系的模拟.我似乎无法弄清楚如何编写用于计算另一个行星施加在行星上的力的方向的代码.这就是我到目前为止:

void attraction(planet self, planet other, float *direction)
{
   float d = vecDistance(self.p, other.p);

   //calc the force of attraction
   float f = G * ((self.mass * other.mass) / (pow(d, 2)));

   //calc the direction of the force
   float vectorDist[3];
   vecSub(self.p, other.p, vectorDist);

   direction[0] = f*vectorDist[0] / d;
   direction[1] = f*vectorDist[1] / d;
   direction[2] = f*vectorDist[2] / d;
}



float vecDistance( float *pV1, float *pV2 )
{
    float fLen=0.0f;

    if(pV1 && pV2)
    {
        float av[3];

        vecSub(pV2, pV1, av);
        fLen=vecLength(av);
    }

    return fLen;
}


void vecSub( float *pV0, float *pV1, float *pVRes )
{
    if(pV0 && pV1 && pVRes)
    {
        pVRes[0]=pV0[0]-pV1[0];
        pVRes[1]=pV0[1]-pV1[1];
        pVRes[2]=pV0[2]-pV1[2];
    }
}

Answer 1

没有输入数据很难知道,但我猜你正在遇到浮点溢出,因为你操作的值太大了.

例如,如果您的单位在SI单位系统中,那么:

venus.mass = 4.87e+24;    // kg
earth.mass = 5.98e+24;    // kg

并且等式中的分子变为:

self.mass * other.mass == 2.91047e+49;

单精度浮点数不能大于3.4e + 38,因此质量积被视为无穷大.

当然高指数在某种程度上相互抵消,因为距离也很大(大约1e + 10)并且引力常数很小,但如果上述产品是一个中间结果,所有依赖结果也是错误的.(-1.0#IND意味着不确定并且对应于NaN,而不是数字,无效的浮点数.)

有几种方法可以解决这个问题:

• 使用可以在浮点状态下安全平方的值.例如,如果使用地球质量对质量进行标准化,则得到的数字大约为1.0,这对于进行计算应该是安全的.同样,一个好的距离单位可能是天文单位.

• 重新排列表达式,以免中间表达式溢出.例如,代替(m1 * m2) / (d*d),写(m1 / d) * (m2 / d).

• 使用double而不是浮动.可表示的最大值double约为.1.8E + 308.请注意,这个问题不会消失duoble,只是为您提供更多的操作空间.举个例子,你应该很好.