calculate direction of force of attraction between planets

Question

I am trying to build a simulation of a solar system in opengl using C. I can't seem to figure out how to write the code for calculating the direction of the force exerted upon a planet by another. This what I have so far:

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

It's hard to know without input data, but I guess you are running into floating-point overflow, because the values you operate with are too big.

For example if your units are in the SI unit system, then:

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

and the numerator in your equation becomes:

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

Single-precision floating-point numbers cannot be greater than 3.4e+38, so the mass product is treated as infinity.

Of course the high exponents cancel each other out somewhat, because the distances are also big (on the order of 1e+10) and the gravitational constant is small, but if the above product is an intemediary result, all dependent results are wrong, too. (The -1.0#IND means indefinte and corresponds to NaN, not a number, an invalid floating-point number.)

There are several ways to fix this:

• Use values that can be squared safely in the floating-point regime. For example if you normalise the masses with the mass of earth, you get numbers around 1.0, which should be safe to do calculations on. Likewise, a good unit for distances might be the astronomical unit.

• Rearrange the expression so that you don't get overflow of intermediary expressions. For example instead of (m1 * m2) / (d*d), write (m1 / d) * (m2 / d).

• Use double instead of float. The maximum representable double is approx. 1.8e+308. Note that this problem doesn't go away with duoble, it just gives you more headroom to operate. For your example, you should be good, though.