计算行星间引力方向。

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];
    }
}

- user3476732

---

力的大小计算看起来没问题，但向量我不确定。return fx, fy, fz; 是行不通的。attraction 返回一个浮点数，而不是3个值。vecLength 定义在哪里？为什么你要用长度为4而不是3来定义 av？ - Millie Smith

我已经修改了代码，但问题仍然存在。我的测试星球没有出现在屏幕上。 - user3476732

我的错误可能在其他地方，但是我对于我的方向计算不确定。 - user3476732

是的...你的测试星球没有出现在屏幕上，而你的计算错误则是两件完全不同的事情...为什么不将你的力量计算结果打印出来呢？ - Millie Smith

你需要决定相对于什么方向施加力。每个质量中心上的力都指向另一个质量的中心。因此，质量#1上的力与质量#2上的力具有相同的大小但相反的方向（分量为相反符号）。 - Peter

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而不是float。最大可表示的double约为1.8e+308。请注意，这个问题并没有通过duoble消失，只是给您更多的操作余地。对于您的示例，应该很好。

- M Oehm