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

calculate direction of force of attraction between planets

我正在尝试使用 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];
    }
}

没有输入数据很难知道,但我猜你 运行 进入 floating-point 溢出,因为你操作的值太大了。

例如,如果您的单位采用国际单位制,则:

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

等式中的分子变为:

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

Single-precision floating-point numbers不能大于3.4e+38,所以质量乘积被视为无穷大。

当然高指数相互抵消了一些,因为距离也大(1e+10量级),引力常数小,但如果上面的乘积是中间结果,所有依赖结果也是错误的。 (-1.0#IND 表示 不确定 并且对应于 NaN 不是数字 ,无效的 floating-point 数字.)

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

  • 使用可以在 floating-point 范围内安全平方的值。例如,如果你用地球的质量标准化质量,你会得到大约 1.0 的数字,这应该可以安全地进行计算。同样,一个好的距离单位可能是 astronomical unit.

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

  • 使用 double 而不是浮点数。最大可表示 double 约为。 1.8e+308。请注意,duoble 并没有解决这个问题,它只是为您提供了更多操作余量。不过,对于你的例子,你应该是好的。