Gravity Simulation Tutorial

Force fields tutorial
---------------------
Imagine you have a star at position (0, 0) and it produces a gravitational field, that pulls everything around it. I'll describe how to compute those forces.

If you have an object at position (x, y), then the distance from the star is d = sqrt(x^2 + y^2)
And the force field produced by a star is the inverse-square force F = -a/d^2
where a is some constant that determines how strong the gravity is.

So all you have to do is resolve this into components:
forceX = F * dx/d
forceY = F * dy/d
(in this case dx=x and dy=y describes the component displacement from the star)

And integrate these as you do for any other physics:
velocityX += forceX * deltatime
velocityY += forceY * deltatime

X += velocityX * deltatime
Y += velocityY * deltatime

Experiment: Create an array of planets, and give them some random position and velocity, and see them orbit the star. Probably some of the orbits will be elliptical and some will be hyperbolic.

Experiment: It might be interesting to create an exploding planet (that explodes into lots of small planets) and see how the debris field is dragged by the star. You may end up with an asteroid belt of some form.


Repositioning the star
----------------------
Let's say the star is now at position (x0, y0).
The distance to the star from a point (x,y) is now d=sqrt(dx^2 + dy^2) where dx = x-x0 and dy = y-y0

Experiment: Try placing the star at some different points.



Two or more stars
-----------------
You merely have to add the force vectors for each star.
For example, if you have stars at (x0,y0) and (x1,y1) etc, then compute
For star 0 : dx0 = x-x0 ; dy0 = y-y0 ; d0 = sqrt(dx0^2 + dy0^2) ; F0 = -a0/d0^2
For star 1 : dx1 = x-x1 ; dy1 = y-y1 ; d1 = sqrt(dx1^2 + dy1^2) ; F1 = -a1/d1^2
etc
(The different a0 a1 values here represent the different gravitating strengths of each star)

To combine the force vectors, you merely sum the components:
forceX = F0 * (dx0/d0) + F1 * (dx1/d1) + etc
forceY = F0 * (dy0/d0) + F1 * (dy1/d1) + etc

Experiment: Try seeing how planets orbit when there are two or three stars.


Stars attract each other
------------------------
Stars are not stationary objects in space. If two stars are close by, they will attract and move one another. In fact, there is a tension force between every pair of gravitating bodies. A true simulation would have to consider every pairing. With N bodies, that is N*(N+1)/2 pairings.

Instead of thinking of stars as being lone "attractors", think of pairs of bodies instead, and the tension force between each pair.

So if a pair of bodies are positioned at (x0,y0) and (x1,y1) then look at their relationship:
dx = x1-x0 ; dy = y1-y0
d = sqrt(dx^2 + dy^2)
F = -a/d^2

Then apply this tension force "equal and opposite" to the two bodies:

Body 1 feels:
forceX = F * (dx/d) ; forceY = F * (dy/d)

Body 0 feels:
forceX = F * (-dx/d) ; forceY = F * (-dy/d)

Understand the role of the - here is because as x0 and x1 switch roles, dx=x1-x0 becomes -dx=x0-x1


Mass
----

In reality, stars and planets have mass, which has a significant effect on the dynamics, due to the principle of conservation of momentum.

We've actually cheated above where we wrote
velocityX += forceX * deltatime

Because actually,
velocityX += accelerationX * deltatime

Recall Newton's second law: that Force = Mass * Acceleration
To compute the acceleration from a force, we rearrange this to give: Acceleration = Force / Mass


The mathematically correct force between two gravitating bodies is Force = -G*M0*M1/d^2
where G is the 'gravitational constant', M0 and M1 are the respective masses, and d is the distance between them.

This force is a tension force that pulls "equal and opposite" on both gravitating bodies. Differing masses mean they may well "accelerate" at different rates.



Body 0 (with mass M0) feels accelereration = Force / M0 = (-G*M0*M1/d^2)/M0 = -G*M1/d^2
Body 1 (with mass M1) feels accelereration = Force / M1 = (-G*M0*M1/d^2)/M1 = -G*M0/d^2

See that the mass of body 0 plays no role in the acceleration felt by body 0. This means that a big

planet or a small pebble exhibit the same acceleration in the gravitational field of the star.


Putting that all together, we have something like this:

G = 1; //gravitational strength

// For N gravitating bodies (0 upto N-1), consider all pairs:
for(i=0; i<N-1; i++) {
  for(j=i+1; j<N; j++) {

    // consider pair i, j

    dx = StarX[i] - StarX[j];
    dy = StarY[i] - StarY[j];

    d2 = dx*dx + dy*dy; // squared-distance
    d = sqrt(d2); // distance

    F = -G/d2; // tension force is inverse-square

    //apply equal and opposite force
    StarVelocityX[i] += F *(dx/d)*StarMass[j]*deltatime;
    StarVelocityY[i] += F *(dy/d)*StarMass[j]*deltatime;

    StarVelocityX[j] += F *(-dx/d)*StarMass[i]*deltatime;
    StarVelocityY[j] += F *(-dy/d)*StarMass[i]*deltatime;
  }
}
// euler integration of velocities into positions
for(i=0; i<N; i++) {
  StarX[i] += StarVelocityX[i]*deltatime;
  StarY[i] += StarVelocityY[i]*deltatime;
}

There are various ways you can optimise that, factorising deltatime/d into F, is the obvious one.
You may want to check if d=0 so that you don't divide by 0.
You may want to give each stars a radius, and test for collisions.
You may want to try something other than euler integration for higher precision physics: eg, verlet, RK4, or Fehlberg.

Enjoy. Comment. And show me you simulations :)
https://twitter.com/_sorceress