Gravity simulator ..... in Ruby

It seem I have written quite abit of Physics related articles already. That seem to pique my interest mostly in recent times.
I was "brooding" that I have to soon change the theme to Biology or Math for a change... but then a cool idea come to my mind, what if I do something that combines Physics, Math and Programming in a single piece. What I'm talking about is making simple gravity simulator. It is fundamental enough, not very hard to do in simplified form. Allows us to explore basic notions and techniques in a novel way and best of all it is fun.
For the task I picked the Ruby programming language. Ruby is very clean and easy language to read/write and understand. Also I was able to find a nice simple 2D drawing framework, another good coincidence.
Soooo... we are going to write simple GRAVITY SIMULATOR.
If you want to take quick pick before you continue here is the code : gsim.rb
or snatch it from github

The Physics

Lets think, what do we need to know to simulate, for example the Solar system.
First thing that comes to mind is that we have to use in some way the Universal gravitational law, right ? So we will try our intuition...
Universal gravitation law : `F = G * (M*m)/r^2` where : `m` : mass of the first object `M` : mass of the second object `r` : distance between the objects `G` : gravitational constant `F` : exerted force next thing we will use also the Newton second law : `F = m*a` where : `m` : the mass of the object `a` : the acceleration the common thing is the force `F`, so then we can do this : `m*a = - G * (M*m)/r^2` the minus sign is here to represent the idea that the force acts towards the originator of the force i.e. gravitation we can multiply both sides by `m`. If we do this we get : `a = - (G * M) /r^2`
Now that is interesting ... once we know the mass and distance to the remote object we can compute what acceleration it will impart on the source object.
If we have just two objects as you can see we can calculate any situation we may have simply by substituting the values.
But what if we have multiple objects ?
In this case we pick one object and calculate all influences the other objects impart, then move to the next object and do the same.. and then the next .. and the next.. So if we have more than two objects we can't just apply the formula, but have to do this iteratively.
This is where the programming part of our project will come in, but first ... let's finish what we started... our final goal in applying those physical laws is to find what should be the next position to which we have to move the object. We got the acceleration already, so we can find the velocity (which is acceleration over time) and from the velocity we can find the distance (which is velocity over time) the object have to be moved.
Remember we are doing this stepwise, which means that when we calculate velocity we assume acceleration is constant, also we assume that velocity is constant when we compute the change in position.
The step in our case is the time frame we accept that the acceleration and velocity are constant. The smaller it is the more correct values we get the bigger it is the faster the calculation and from there the animation speed. (You could experiment later with time intervals and you will see that sometime you can get drastically different results if you pick longer time intervals).
OK here are our formulas :
`a = v * t` and : `d = v * t` where : `v` : is the velocity `t` : is the time passed `d` : is the travelled distance
OK, now we have all the ingredients , if we know the mass of the objects and the distances between them we can iteratively go trough all and calculate how much any of them have to move. Then we do it for the next time frame and do it again and then again ... ad infintum. There is one more catch still, up to now we assumed the law applied in just one dimension i.e. we can calculate how the objects will move back and forth closer or away from each other as if they live on a single line.

If you want to dig a bit deeper in the Physics of gravity you can click on the [+] sign below. As part of the discussion of Energy I picked Gravity to explain Kinetic and Potential energy :


  Kinetic and potential energyclick to toggle>  [+]
We will use Gravitation to explore the idea of different forms of Energy in Classical mechanics.

To make useful the concept of Energy in Classical mechanics (Remember the idea of Energy comes from Thermodynamics), we would need at least two type of energies, the reason for this is that energy can only be converted from one form to another but it cannot be created or destroyed.
The scientists decided that those two candidates would be called Kinetic (the energy of motion) and Potential (the energy of position/state/"difference of position"). It make sense because Classical mechanics is primary concerned with the description and study of motion of objects.
On the right side in the image you would see all the various equations related to gravitation. I've decided to make this diagram, because most of formulas look very similar and at times it get's too confusing, use them as reference. And don't forget they are for circular orbits.
Now Kinetic energy has to measure the dynamic part of motion.

We know that the change in energy of a system is equal to the Work done on the system. But in Newtonian mechanics Work is equal to the Force applied over Distance. And we also know, force from Newton second law is `F = m*a`. "Mixing" them tighter will give us the following formula for Kinetic energy : (TODO: Add derivation) (1) `E_k = 1/2 *m*v^2` `m` - mass of the object `v` - velocity of the object

Potential/Positional energy on the other hand does not have exact formula. In sense it is whatever we ascribe it to be. How's that ? Let's take for example the Gravitation field, which is a conserved field.
Which in our case means that the Total orbital/gravitational energy `(E)` of an object moving in this field is constant, instead what changes is the ratio between Kinetic and Potential energy i.e. `E = E_k + E_p`. Now we are going somewhere. To figure the exact formulation for Potential energy we would need to use the Universal gravitational law, because as we elaborated earlier change in Energy is the Work-done (and work depends on the gravitational Force).

Everybody knows the G-law : `F = G * (m_1*m_2)/r^2` where : `m_1, m_2` - mass of the two objects `r` - distance between the objects `G` - gravitational constant `F` - exerted force

We will look at the case when only 2 bodies are involved in which `m_1` is much bigger than `m_2` i.e. the influence of `m_2` onto `m_1` is negligible. (Ex: Earth orbiting the Sun, rocket orbiting Earth ....)
So we will say that `m_2` is bound to orbit around `m_1`. And using gravitational law we can calculate what Force will be applied to `m_2` at any point in space we wish. (In fact that is what the idea of gravitational field mean, if you were wondering i.e. we can calculate the Force acted upon any body who happens to be at any point near `m_1`.). Another simplification we will make is that `m_2` orbits in circular orbit.
OK. Let's rehash. If we look at the G-law we will see that the Force acting on `m_2` will be smaller the further the satellite is from the center of the more massive body. Our task is now to figure out the Work we have to do to move object from one orbit to another, this will give us the differences of the energy of the two orbits i.e. the potential energy.
See it was very easy...
Let see how our derivation will look like.
We start with the formula for Work : `W = F * r` Then we use the Work-energy theorem (which states that "change in Kinetic energy is equal to the Work done") : `Delta E_k = W` `Delta E_k - W = 0` And according to Conservation of Energy if there is no outside influence : `Delta E = 0` but : `Delta E = Delta E_k + Delta E_p` so finally : `Delta E_p = - W = - F * r` `W` - work done `F` - force applied `r` - distance over which the Force is applied

Because we need the difference between Work needed to move body from orbit at distance `r_0` to orbit `r_1`, we would use little bit of Calculus. Don't run away, screaming :), let me try to explain it in simple terms.
If we were working with simple difference between two values we would have used simple algebra and subtracted the values calculated between the two points (Normally this is the so called `Delta`-delta, you probably have encountered often in high school physics). But as you may deduced the gravitational influence from point to point is changing gradually, because we can pick any point in between with any precision we like..
So if we calculate the `Delta` we will get just approximation not the real value.
OK you may say, what if we calculate differences every kilometer from one orbit to the other and then sum the list of results. Good point, we will get better approximation than the first case. But if we continue on this path we can go even further and do the calculation every meter, or every centimeter, or even every millimeter.
The smaller the interval the better approximation we will get, but you can see easily it will take us infinity to do all those calculations.
To save us from our misery sir Isaac Newton invented Calculus. What this mean for us is that instead of doing calculation for the infinite steps from orbit one to orbit two, we just apply neat simple formula, which the nice mathematicians cooked for us. And this formula is (definite integral, there is also indefinite integral which is even easier, but remember we are looking for a range between two orbits) :

`int_(x_0)^(x_1) 1/x^2 dx = (- 1/x_1) - (- 1/x_0) = 1/x_0 - 1/x_1 = 1/(x_0 - x_1) `
You can translate it to something like this : split the interval between ZERO and `x_0` to infinite steps and calculate `1/x^2` for all those steps. After you are finished sum the results. Do the same for the interval from ZERO and `x_1`, sum again. Now subtract one sum from the other. Done.
Look at the integral sign it even looks like stylized SUM sign.
Mind you what really happens is not exactly summing, but it is close enough for our current understanding.
So let's apply this formula to our work equation and see what happens.
By the rules of calculus the constants can be moved outside the integration process. So we move `G*m_1*m_2` out. `W = int_(r_0)^(r_1) F\ dr = int_(r_0)^(r_1) (G * m_1* m_2) /r^2\ dr = (G * m_1* m_2) * int_(r_0)^(r_1) 1/r^2\ dr = (G*m_1*m_2) /(r_0 -r_1)` and as we will see in a minute the ZERO point for Potential energy is chosen at infinity i.e. both distances are negative, so `Delta r = (-r_0)-(-r_1) = r_1 - r_0`, but I'm getting ahead of myself. Our solution for Potential energy then will be : `Delta E_p = - W = - (G*m_1*m_2) / (Delta r)` So this then, is our formula for gravitational potential energy : (1) `E_p = - (G*m_1*m_2) / r`
Let me go on a tangent here ... you can skip it if you want.. we already got our formula.
SPECIFIC "thing"
Normally in Physics where mass is involved there is this notion of 'SPECIFIC'-something, which is calculated by dividing the formula we are exploring by the mass. The trick is that "omitting" the mass in this way, we gain a formula for a unit-of-mass i.e. making it independed of the exact object mass. F.e. if we divide the above formula (1) by `m_2`. (And using the symbol `M` instead of `m_1` for the big mass), we get : (2) `E_p = - (G*M)/r` You probably would say, so what's the difference ? The difference is that now in (2) we have formula were we do not care about the exact mass of the satellite. Or alternatively you can think of it like the formula for potential energy per 1 kg i.e. we have invariant way of doing it. So if we need to calculate what potential energy will be for 150 kg for example, we just multiply the result of our calculations by 150 at the end when we are finished. Ok to confuse us even more physicists do not call (2) "specific potential energy", but "gravitaional potential". And the earlier formula (1) is "gravitaional potential energy". Still confused... welcome to the club ;) Can you spot equation (2) in the graphic-diagram above ?

Pfuu... finally we have the formulas for both Kinetic and Potential energy in Gravitational field, but we are not done yet.

Escape velocity

TODO

The Math

What we need to do now is to extend our reasoning so we can handle 2 or more dimensions. We are lucky because there is a mathematical tool to do just that and it is called Vectors. If you are not familiar with the concept open the box below for introduction, otherwise continue ...


  Vectorsclick to toggle>  [+]
What are Vectors and why we need them !?
Good question, vectors seem so deceptively simple when we initially meet them, but then when we dig deeper they tend to become confusing. I think the basic reason for it is that Vectors are introduced in school as pure mathematical "machination", rather than approaching them in more intuitive practical way.
In this article I would like to concentrate more on the practical matter about vectors, which is often unexplored, rather than purely mathematical.
Before Vectors we had Scalars.
Scalars hold single value, but often there is physical quantities that can not be represented by single value or if we can then managing the math becomes tedious to say the least.
Think about velocity, acceleration, position in 2 dimensions or 3 dimensions .....
The fundamental idea behind vectors from practical standpoint is that they can hold/represent multiple values, but still be expressed with single symbol in formulas and calculations.
The big win from this is that we can use ordinary arithmetic and algebra for the most part and still do correct computation for much more complex physical phenomenas. Why throw used and tried practices if we can reuse.
Let me repeat this because I think it is the most important idea behind Vectors.... we can use the same algebraic manipulations we learned in high school to do vector arithmetic's. In most cases we treat the vector as a "scalar-black-box" and when we don't there is pretty good common sense reason for that, mostly expressing new ideas which are not available to us in arithmetic's. So of course if we can't represent those new ideas in scalar arithmetic, then we will invent new conventions.
We represent vectors in Cartesian coordinate system by their x and y coordinates (2D). The assumption is that the vector starts at coordinate [0,0] and end at [x,y] like in the image beside. To calculate the magnitude of a vector we use the Pythagorean theorem.
`v^2 = x^2 + y^2` so : `|vec v| = sqrt(x^2+y^2)` The vertical bars around the vector name signifies magnitude.
Then there is the postulate that a sum of two vectors is a new vector that starts from the beginning of the first vector and ends at the end of the second vector. In our case this is `vec v` in the diagram. Here is something interesting, vectors if transposed are the same, so `vec y = vec y_1`. Then if we apply the rule that we just stated about summing vectors `vec v = vec x + vec y_1 = vec x + vec y`.
That is good pattern but we can do better. When we measure things in the real life we normally say something is so and so meters/ft long. What we are really assuming is that there is somewhere some ideal meter/ft which is "set in stone" with a value of 1 (one) and we use this ideal-meter to measure and compare things.
So then you probably guessed that we will have something similar with Vectors, yep! we have such thing and it is called unit or normalized vector. Now that we have this ideal-measuring-stick device we can say that, vector is equal to the magnitude multiplied by the unit-vector :
`vec x = |vec x| * hat x` the `hat x` thingy here signifies the unit vector and `|x|` is the magnitude i.e. the vector is equal to (magnitude * unit-vector).
The next thing that comes to mind now is that we can represent our original `vec v` as a sum of multiplications of magnitudes and unit vectors, here is how :
Variation 1 `vec v = |vec x| * hat x + |vec y| * hat y` In many books and tutorials instead of putting the arrow over the top over the vector symbol, bold typeface is used, like this : Variation 2 `bbv = |bbx| * hat x + |bby| * hat y` we wont because it is confusing. Other way of representing vectors via unit vectors are the following : Variation 3 `bbv = x*bbi+y*bbj` Example : `bbv_1 + bbv_2 = (x_1 * bbi + y_1*bbj) + (x_2*bbi + y_2*bbj) = (x_1+x_2) * bbi + (y_1 + y_2) * bbj` here `bbi` and `bbj` are the two unit vectors. See this time the "name" tells us which are the unit-vectors, not the boldness or hat. Another way is using pairs of values like this : Variation 4 `bbv = (x,y)` Example : `bbv_1 + bbv_2 = (x_1,y_1) + (x_2,y_2) = (x_1+x_2, y_1+y_2) ` One of the reasons for the confusion with vectors are those different way of representation. I hate it it sucks big time. ( :) can I use "it it" like this ?) Keep in mind that there are variations on the above variations.
Why would we complicate our lives like this ! I'm talking about using unit-vectors not the complications of the representations :).
The benefit is that we can deconstruct the original vector into two scalars. We also can see that unit-vectors associated with every of the scalars are perpendicular to each other. The more general term used is that the vectors are orthogonal to each other i.e. mutually-independent, non-overlapping.
I think now you can appreciate the idea which is : we can represent any vector via 2 mutually independent variables (three for 3D, etc).
Let's list the three important concepts I wanted to convey about vectors :
  • how to calculate magnitude
  • What are unit vectors and what are they good for
  • the concept of orthogonality
As I said I will skip the other common mathematical manipulation like subtraction, multiplication by scalar..etc.. there are plenty of tutorials on Internet about those. What I want to concentrate your attention to are new, as I promised purely Vector operations and their genesis and purpose.
I'm talking about :
  • dot product
  • cross product

TODO:
- dot and cross product

Physics with vectors

As you know or you just figured out using vectors allows us to use the same or almost the same algebra rules we learned in high school, to handle 2D/3D otherwise we have to use trigonometry and separate calculations.
So let's convert our acceleration equation to vector format. For this we preserve the original equation which gave us the magnitude of the Force that was imparted on an object, but we will need to add something else to account for the direction.
Modified acceleration equation : `vec a = ( - (G * M) /r^2) * (vec r/ |r|)` The first part which we know already is the magnitude of the acceleration the second part is a vector of unit length which holds the direction of the calculated force. Vector of unit length in this case is like `+- 1` in arithmetic, i.e. if we multiply by unit-vector we do not change the magnitude but only the direction. If you look carefully it is represented by the vector itself divided by its magnitude... which seems reasonable. Then the other two equations change in the following way : `vec a = vec v * t` and : `vec d = vec v * t` almost the same syntax, of course under the hood vectors holds two or more values. You could also see that the distance became vector too, which is what we would expect, because now we have change in position in more than one axis.
So far so good we are done with the physics and math part. Wasn't hard, was it ? Our next task is to translate this to a working computer program.

The Programming

For our purposes we will use Ruby. It is very easy and straightforward language yet very powerful. I seem to think this is the best introductory language for anyone trying to start with computer programming. By far much nicer, easier and more forgiving than Java and C#, for a first language.
Next because we want to build visual simulator we have to pick some 2D/3D drawing framework. It has to be easy to program and understand, because we are constrained on how long we can make this article.

2D : Green Shoes

When I first started thinking of writing this article I went trough couple of choices, but because either the framework was too low level for our purposes or for lack of documentation I ended up choosing the "Shoes"-framework.
This is very easy to use framework with its own domain specific language seamlessly integrated in Ruby. It also has a very interesting history of development with currently several spin-offs.
The original is called "Red shoes" and is based on GTK2. It is implemented in mix of C and Ruby and programs developed using it should run successfully on Linux, Mac and Windows.
But instead "Red shoes", I went developing this small app using the "Green shoes" variant. (It should work on Red Shoes too). The reason again being that "Green shoes" is much easier to install (so more people will be able to experiment) also Green shoes is written entirely in Ruby.
Here is what you need to do to install it : [Installing] if you have Ruby installed already of course.
And if you are interested in the other flavors of Shoes check this link : [The list of Shoes implementations]
I would also like to thanks the creator of Green shoes ashbb who was very responsive and helped me solve a problem with the planet-trails.
OK we decided it will be Green shoes... let's continue...

The program

The simulator is pretty heavily commented, but let's discuss how did we put our physics in Ruby code. Remember we decided to use Vectors and yes Ruby has them, one less thing to worry about, we just need to "require" them. Let's experiment a little with interactive Ruby console and see how they work :

# start the ruby console
$ irb

#We first need to to load the 'matrix' module, which contains the Vector class :

> require 'matrix'
 => true

#Then let's create two vectors (let's pretend they point to 
# the position of the two gravitationally interacting objects)

> v1 = Vector[5,5]
 => Vector[5, 5]
> v2 = Vector[7,3]
 => Vector[7, 3]

#if that is the case we can find the distance difference between them (remember vector distance):

> r = v1 - v2
 => Vector[-2, 2]

#Of course what we get back is another vector, if we want to get the scalar quantity 
#  i.e the distance we will use the .magnitude() method.

> r.magnitude
 => 2.8284271247461903

#On the other hand if we want to find the direction (described by purely unit vectors) of the force we will .normalize() the distance vector.

> r.normalize
 => Vector[-0.7071067811865475, 0.7071067811865475]

#We can see that if we just do the following experiment with the "unit vectors themselves":

> Vector[-1,1].normalize
 => Vector[-0.7071067811865475, 0.7071067811865475]

#we get the same result...
That is all we need to know to be able to build our program.
The main point again is that we do all the calculations with good old arithmetic methods.
The code is divided in three parts : the front-end program, the Gravity class that does all the physics and the Planet class that holds the planet info, draw the planet and handles the scaling from real world units of meters to pixels.
(The green shoes variant does not support yet translation/scaling, so we have to do it ourselves.. dont worry it is no biggy).
The core of everything is implemented in two routines inside the Gravity class.

You can download the script here : gsim.rb
#This the central piece of the whole simulator.. what it does is calculate
# the velocity change caused by the influence of the remote object

def compute_forces (pa,pb,dt)
  #first calculate the distance between the two objects
  r = pa.pos - pb.pos;#vectors
  #Using Gravitation law and Newton second law, find the acceleration
  pa.acceleration = ((- G*pb.mass) / r.magnitude**2) * r.normalize
  # ...now that we know acceleration, do step-calculation of the velocity
  pa.velocity += pa.acceleration * dt
  # ..... what is left is to calculate the new position.....
end

......

#now that we have the cumulative velocity calculated, find us a new position
def new_position p,dt
  p.backup_pos if @trail and @i % 2 == 0
  p.pos += p.velocity * dt #it happens here
  p.trail if @trail and @i % 2 == 0
  p.draw
end

Did you saw almost verbatim from our previous discussion. Except may be because this is iterative process we accumulate the changes in velocity and postion via this operator +=. You should not expect the velocity and position to grow to infinity, because they are now changing along 2 axes and most of the time because two objects attract each other they will tend to point along the line connecting the two objects in straight line.
OK.. the compute_force() is the meat of the piece.. it does the calculation of the gravitational influence of one body on another. So in the main Shoes program on every iteration of the animation the call is made to the wrapper function called compute_all() which loops over all the objects one by one using compute_force() to do the computation. Once done i.e. the cumulative velocity change was calculated it calls new_postion() which in turn calculates the distance and direction the object has to be moved and then draws the planet in the new position.
One other thing which may be interesting to know is that all the calculations are done in real units like : kg, meters, Newtons...etc, but for animating them we scale to pixels to account for the size of the window surface we have at our disposal to draw. This is done in the Planet class there is two methods for calculating the x() and y() position.
Don forget that x,y are stored as Vector.
Here is how those two methods look like :
#get the distance in meters then scale it to fit into the screen size
# and because Shoes coordinate zero point is the top left corner we have to account
# for that. (Green shoes does not support translate method yet)
def x
  (@pos[0] * @scale) + @center_x
end
def y
  (@pos[1] * @scale) + @center_y
end
The one other small tidbit is calculating the correct scale factor, here is how we do it :
.....
unless @scale
  #calculate the scale factor
  @max_distance = data.map { |h| h[:distance] || 0 }.max
  #half the height should be equal to the max distance
  @scale = @height/(2*@max_distance*1.1)
end
.....
divide the height of the window drawing surface by two. (Normally the height is the smaller than the width that's why I chose it, feel free to modify it as you wish in your experiments).
Then divide the result to the distance of the furthest planet. And finally add 10% fluff.
The simulator does not have any UI, but it is highly configurable via command line arguments, here is the list :
$ ruby gsim.rb -h

 -h : help
 -f : total frames to play (def:1200)
 -p : frames per second, dont push it :)
 -s : scale factor to fit the system on the screen, higher value is closer view
      100 will scale you to the inner Solar system
 -t : time step (def: 10 days).
      Smaller value give more accurate results, but is slower.
 -b : background (def: black)
 -r : leave trail ..... you get nice tail marking the path of the planet.
      But if you push it program may become sluggish.
 -i : adjust the scale factor to see the inner Solar system (scale:120,time-step:1)

Of them the interesting one is probably the -i switch which set the scale factor so you can see only the inner Solar system.

You can download the script here : gsim.rb or github

Further things to think about

Now that we got our code running we can do some experiments.
For starters I created in the data Hash planet called "Havoc" ;), just uncomment it and run the gsim with -i switch.
Or you can uncomment the planet "Postion" this example shows how you can position a planet at specific location via :pos. If you set the :distance instead of :pos the planet is set at position [0,distance]. (The center of the window is the zero point [0,0]).
For example set the Sun position to [-300e9,0] and watch chaos ensues...
As exercise try to extend this small demo to 3D. It should not be very hard because all the calculation as we saw are vector based. See using vectors paid off ;). You just need to extend them to support 3D. Second convert the 3D coordinates into 2D (for drawing) in the Planet class, inside the methods x() and y() ... and that should probably do it I hope.