When KSP was published, planetary mechanics got quite some attention. Coding this on your own is not hard at all and requires just some basic mathematical knowledge. During this tutorial you will learn everything you need to know. A c# console project accompanies this tutorial and can be downloaded from google drive (zip, 40kB). It contains all the code described here and you may use, adapt, modify and use it in your game.

Equations of motion

Let's start simple: Imagine a particle which can move only in one dimension (so to speak: on a line). To describe the motion of this particle all you need are some variables. The first one is, obviously the position x of the particle. That would be all we need, if we just have a look at it at one single moment in time. As we are interested in the so called trajectory of the particle, we need some more variables. One of them is obviously the time t. The other variable required is the velocity v of the particle. The position and velocity would be all we need, if the particle is just moving without any interactions. If we want to have the particle (satellite) be attracted by other objects (planets) we need one more variable which is the acceleration a. In Physics you can show that those three variables are enough to describe the movement of the particle completely deterministic.

Let's have a look at those variables and the equations, which connect them on physical level! As you know, velocity is the change of the position with time. Think about driving a car: your speed (the velocity) is given in kilometers per hour (for the metric unit system) which is exactly what I just wrote: The change of the position in a given time. When you drive 60km/h for one hour, your position will be different by 60 km after that hour. To put that into an equation, we can easily write:

v = delta_x / delta_t

Note that the delta corresponds to a difference. Think about your car at two different points in time: t1 and t2. delta_t is just the difference between them: delta_t = t2 - t1 The same holds for delta_x.

The acceleration is the change of the velocity with time. Think about ads for cars: "Accelerates from 0 to 100 in 15 seconds." 0 and 100 are known to be velocities. So it can be written:

a = delta_v / delta_t

For our one dimensional case all of those variables are numbers (double), but we can easily see, that all the relations are still true if we exchange the scalar values by vectors (of dimension two or three), as the different axis (x, y, z) can be calculated independently.

Setting up those equations

In the final code, we will be interested in the position at a given time. Therefore some refinements of the equations are required. We know the acceleration. from that we can compute the velocity:

delta_v = a * delta_t

as this is just the change in velocity the real velocity will be

v_new = v_old + delta_v

or as it will be written later:

v += a * delta_t

The same is true for the position and the velocity.

delta_x = v * delta_t
x_new = x_old + delta_x
x += v*delta_t

Those equations now tell us, how a new state of the system will look like, based on the old state of the system. Mathematically this is called "integrating the equations of motion".

Newton's axioms

For particles to interact(planets should attract each other) we need to know, how the forces translate into those equations of motion, we just talked about. This physics building block can be derived from Newton's axiom. One of them states that Force equals mass times acceleration. If we put that into an equation it can be written as:

F = m * a.

OK, so what are those F and m values. F is the type of the force we are looking at. If you want to see satellites move around planets, you need the law of gravity. If you want to have a look at electrons and protons, you need Coulomb's law. All we should care about now is that F tells you "what" is acting on the particle and only depends on the position of the particles. The mass m is just the mass of something you can pick up. The more mass it has, the harder it is to make it move. If you want to lift a light feather (low mass), that is easy. If you want to lift a large heavy crate, that's hard. The mass is sort of the "resistance" of the object against being moved.

As we know the Forces of the particles, we can now use that to calculate the acceleration of the particle.

a = F/m

The Gravitational law

So how does gravity work? I can not give you an explanation to that on a fundamental level, but for what is required here you just need to know that two masses (m1 and m2) attract each other by the following law:

F_G = G * (m1*m2)/(r**2)

G is called the gravitational constant and r is the distance between the two masses (their centers). This distance can be calculated by Pythagoras' theorem in any dimension.

Show me some code already!

So let's put all that physics knowledge to some code!

We need a PhysObject class, which will store those variables and do the computations for us. The constructor just initializes the values. The addForce function will let a force act on the object. The most important part is the update function. Let's have a closer look at that! The update function puts the equations mentioned above into order.

public void update ( double dt )
{
    vx += ax * dt;
    vy += ay * dt;

    x += vx * dt;
    y += vy * dt;

    ax = ay = 0;
}

Those are exactly the equations explained above. At the end of the function, we reset the accelerations. This is required, since the particles have moved. As the Forces depend on the positions, we need to calculate them in each step.

So in the main function, we need multiple particles, so we create a list of PhysObjects. Then we can add as many particles with any starting conditions we like.

System.Collections.Generic.List<PhysObject> pobjs = new List<PhysObject>();
{
    PhysObject p1 = new PhysObject();
    p1.x = 0;
    p1.y = 1.4;
    p1.vx = 1.2;
    pobjs.Add(p1);

    PhysObject p2 = new PhysObject();
    p2.x = 0.2;
    p2.y = 0.2;
    p2.m = 100000;
    p2.vx = 0;
    pobjs.Add(p2);
}

Then we need an update loop which calculates as many integration steps as we want. In every step we need to determine the forces which are acting on each particle. This is done by looping over all the other particles and calculating the pairwise forces. The Force calculation depends on the distance, which we can easily calculate. After that we normalize the vector from object 1 to object 2 and use those values for the gravitational law. This yields the force. As we know "actio = reactio" from Newton, we can use this to apply the same force on both particles, just with different signs. This helps to reduce computational cost, as the second iteration loop can start at m+1. All previous particle pairs have already been treated completely.

for (int m = 0; m != pobjs.Count; ++m)
{
    for (int n = m+1; n != pobjs.Count; ++n)    // note that all particles <m have already been treated.
    {
        PhysObject p1 = pobjs[m];
        PhysObject p2 = pobjs[n];   // get the objects

        double dx = p2.x - p1.x;    // distance in x and y
        double dy = p2.y - p1.y;

        double r = Math.Sqrt(dx * dx + dy * dy);    // euclidean distance

        dx /= r;    //normalized vector
        dy /= r;

        double c = p1.m * p1.m / r / r; // law of gravity is \propto m1*m2/r**2
                                                    // G is set to 1 without loss of generality.

        p1.addForce( dx * c,  dy * c);  // adding forces symmetrically. 
        p2.addForce(-dx * c, -dy * c);
    }
}

After we have added all the forces, we can now call update. We cannot do the update earlier, as each particle has to know the forces which are acting on it.

// now the forces are calculated. Do the integration step now.
for (int m = 0; m != pobjs.Count; ++m)
{
    pobjs[m].update(0.01);
}

A word about numerics

What is written in the code is the most basic numerical integrator, which is called after the mathematician Leonhard Euler. There are multiple versions of that (implicit, explicit) and many other sort of integrators (Midpoint, Runge Kutta (of arbitrary order), Leap Frog, Velocity Verlet, ... ) and all of them have their benefits and drawbacks. The benefit of this explicit Euler method is that it's simple and that it is computationally very efficient. Unfortunately it is not stable at all. For a large system (many particles, and many in the sense of ten thousands), the largest computational effort will be determining the forces, as each body interacts with all other bodies. This is called an N-squared loop. There are ways to boil that down, but they are out of the scope of this tutorial.

Another important point is the time step deltat. The integrator works well if deltat is small. If it is too large, strange things might happen. One of the drawbacks of the Euler method is, that it requires quite a small time step to work. So if strange things happen with your physics, try reducing the time step first.