Solving Gravity

Today, I spent some time working on the realism of our simulation. I wanted to make sure the physics were realistic. I used the equation for centripetal force and set it equal to the equation for gravitational force and solved for the needed tangential velocity for a certain radius and center mass. After trying this out in our current simulation, I found that it wasn't working the way I expected. I created a couple different values for the radius and then experimented with the velocity until I got a perfect orbit (the radius and velocity stayed constant). I found that for a sun with a mass of 100 and a planet at radius 100, the planet's velocity would need to be 100. These units are all arbitrary and our gravitational constant (G) is 1. According to the math, v should = sqrt(GM/r), and the needed velocity should have been 1. This was obviously not the case. I found also that for the same mass, at a radius of 25, the velocity was 50. After some more experimenting, I plotted some of the points in excel (see the graph below) and found that for a constant center mass of 100, the equation for velocity was roughly 10*sqrt(r).

I noticed that the 10 out front was equal to sqrt(M) for M = 100 and I hypothesized that v = sqrt(Mr). I then created a function in Java called orbit() that takes in a mass and radius and returns a velocity. I set this up to make the planet orbit at a variety of different values of M and r and it worked perfectly. This made me think that perhaps I had messed up the force of gravity calculation in the code and the simulation we were seeing was not accurate. At first glance, I couldn't find anything wrong with it but after stepping through with the debugger, I found the problem. Here was the original code:

Vector2 f = r.nor().scl((float) (m1 * m2 / Math.pow(r.len(), 2)));

This is the basic formula for gravity in vector form. The two masses are multiplied together then divided by the magnitude of the radius squared and finally everything is multiplied by the unit vector in the direction of r. The problem here was that calling the function r.nor() didn't simply return the unit vector of r. It also changed the value of the vector object for r in memory. So then later in that same line when we took the magnitude of r with r.len() we were actually taking the magnitude of the unit vector instead, which was 1. This made get rid of the magnitude of the radius entirely from the formula for the force of gravity and it would just equal (Mm)(r-hat). Then when we calculated the velocity for a circular orbit, we would get v=sqrt(Mr) instead of v=sqrt(M/r) like it should have been. I fixed the code by declaring the magnitude and unit vector before plugging them into the formula. I made sure to put the declaration for r_mag before r_hat so that the magnitude of r wouldn't be affected by the normalizing of r. Everything is now working the way it should.

float r_mag = r.len();

Vector2 r_hat = r.nor();

Vector2 f = r_hat.scl((float) (m1 * m2 / Math.pow(r_mag, 2)));