5

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];
    }
}
user3476732
  • 75
  • 7
  • 1
    The calculation of the magnitude of the force looks right. I don't know about the vector. `return fx, fy, fz;` is not going to work. `attraction` returns a float, not 3. Where is `vecLength` defined and why are you defining `av` with a length of 4 instead of 3? – Millie Smith Jan 09 '16 at 03:40
  • i have amended my code but the problem still persists. my test planet does't appear on screen – user3476732 Jan 09 '16 at 03:56
  • my error may lie else where but just not sure about my calculation for direction – user3476732 Jan 09 '16 at 03:59
  • Yeah... your test planet not appearing on screen and your calculation being off are two completely different things... Why don't you print out the result of your force calculation? – Millie Smith Jan 09 '16 at 04:00
  • variable d displays -.1.#ND00 – user3476732 Jan 09 '16 at 04:15
  • 1
    You need to decide what the direction of force is relative to. The force on the centre of each mass is toward the centre of the other mass. So the force on mass #1 is the same magnitude but opposite direction (components the opposite sign) to the force on mass #2. – Peter Jan 09 '16 at 04:18
  • What do you *expect* the force to be? (as a magnitude, not a vector). Can you also post all your numbers? – Millie Smith Jan 09 '16 at 04:19
  • i am trying to calculate the x,y,z of the force being exerted on a given planet, i would need a vector i'd assume – user3476732 Jan 09 '16 at 04:27
  • Yes, but before the vector, you need a *coordinate system* to make the vector meaningful. With gravitational attraction, you won't be bothered with coordinate system transformations, but you will need to fix the *center-of-mass* for each of your planets within the coordinate system in order to compute the relative force vectors between them due to gravity. Once you have your planets fixed at `x,y,z` locations, simple vector subtraction will give you the relative vector between the bodies fixing the direction of your force vector. – David C. Rankin Jan 09 '16 at 06:17
  • See [Is it possible to make realistic n-body solar system simulation](http://stackoverflow.com/a/28020934/2521214) especially the last part about integration precision ... and you should probably port at least to 64 bit `double` instead of float. – Spektre Jan 09 '16 at 10:01

1 Answers1

6

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.

M Oehm
  • 28,726
  • 3
  • 31
  • 42