Calculate the direction of gravity between the planets

Question

I am trying to build a simulation of the solar system in opengl using C. I cannot figure out how to write code to calculate the direction of the force acting on another planet. This is what I still have:

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 is hard to know without input, but I think you are working with floating point overflow, because the values ​​you are working with are too large.

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; 

The single precision of floating point numbers cannot be greater than 3.4e + 38, so the mass product is considered to be infinity.

Of course, the high exponents partially cancel each other, since the distances are also large (of the order of 1e + 10), and the gravitational constant is small, but if the above product is the final result, all the dependent results are wrong, too. ( -1.0#IND means indefinte and matches NaN , not a number, an invalid floating point number.)

There are several ways to fix this:

• Use values ​​that can be safely squared in floating point mode. For example, if you normalize the masses with the mass of the earth, you get numbers around 1.0, which should be safe for calculations. Similarly, an astronomical unit can be a good unit for distances.

• Reorder the expression so that you do not overflow intermediate 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 does not disappear with duoble , it just gives you more options for work. For example, you must be kind.