Sun, Earth, and Moon in Codea
I spent a few hours last week-end coding a very simple simulation of
the gravitational interaction of the Sun, Earth, and Moon system. Even
if it’s very basic, I was so proud of it I posted a link to it on
Twitter, which lead to the great honor or having
Dr Drang himself ask for more
details. Here they are, along with the code. (You can find the whole
code as a gist, if you want to play
with it. The code below has additional line breaks to make it look not to bad.)
The simulation is written with Codea,
a programming environment on the iPad. Codea is currently named Codify
but this will change when Apple approves the latest version (soon,
hopefully).
The programming language Codea uses is Lua,
with a few extensions: a class system, a drawing loop very similar to
the one of Processing, and ways to interact
with the iPad (touches, accelerometer, and so on). I did not know Lua
before playing with Codea, so please excuse the ugly code that is
coming.
The code
The goal of the simulation is to have several “bodies” which interact
through gravity only. I don’t assume much knowledge of Codea, but if
things are not clear there are some great
Codea tutorials online.
Let’s start with the class describing these bodies.
Body = class()
Our Body object creation method takes the initial position, speed,
body mass, and color, and sets some instance variables.
function Body:init(x,y,vx,vy,m,c)
-- you can accept and set parameters here
self.pos=vec2(x,y)
self.speed=vec2(vx,vy)
self.mass=m
self.color=c
end
The vec2 object is a point. It has several interesting methods
including one that computes the distance between two points. We make
use of it in the next method, where v is a vec2.
function Body:dist(v)
return self.pos:dist(v)
end
We now get to the first method which is about physics:
Newton’s second law of motion
“F = m . a” (force is equal to mass times acceleration, where
acceleration is the rate at which speed changes:
a = (new_speed - old_speed)/time).
This method applies a force f for a given time t on the body,
which results on the body’s speed changing.
function Body:applyForce(f,t)
self.speed = self.speed + t * f / self.mass
end
Since speed is the rate at which position changes, we can describe how
the position evolves during some time t.
function Body:updatePos(t)
self.pos = self.pos + t * self.speed
end
The last two methods are about the graphical interface. The first one
returns the on-screen position of the body. The parameter z is the
zoom factor (how many distance units are in a pixel), and c is the
position that is at the center of the screen.
function Body:screenPosition(z,c)
return vec2((self.pos.x-c.x)*z + WIDTH/2,
(self.pos.y-c.y)*z + HEIGHT/2)
end
Finally, we draw a body as a small circle of the body’s color at the
body’s screen position.
function Body:draw(z,c)
local position = self:screenPosition(z,c)
fill(self.color)
ellipse(position.x,position.y,10)
end
Let’s now move to the next class, which will be at the heart of the
simulation.
Bodies = class()
An instance of this class is an array of the bodies we want to
simulate.
function Bodies:init()
self.blist={}
end
function Bodies:addBody(b)
table.insert(self.blist,b)
end
function Bodies:bodyAtIndex(i)
return self.blist[i]
end
The central method of the whole program is the next one. Given a body
(identified by its position in the array), compute the force the other
bodies exert on it. We accumulate this force in force, and heavily
rely on the vector operations (addition or subtraction of vectors,
multiplication and subtraction by a scalar) to compute it.
More precisely, we enumerate the bodies of the simulation (in the
for loop), and for each one, called b, we compute the force it
exerts on body. This force is a vector from body to b whose
length is (according to Newton’s law of universal gravitation)
G * m1 * m2 / d^2, where m1 and m2 are the masses of the bodies,
d is the distance between the bodies, and G is the gravitational
constant.
The code is almost exactly this formula, except that there is no
G. As I’ll explain later, we choose our units such that G is equal
to 1. The formula (b.pos - body.pos)/d is the vector from body to
b whose size is 1. We multiply it by the masses of both bodies,
and divide by the square of the distance.
When we’re done accumulating the force, we simply return it.
function Bodies:computeForce(bnum)
local body=self.blist[bnum]
local force=vec2(0,0)
local d
for i,b in ipairs(self.blist) do
if i ~= bnum then
d = body:dist(b.pos)
if d ~= 0 then
force = force +
((b.pos - body.pos)/d) *
body.mass * b.mass / (d*d)
else
print("null")
end
end
end
return force
end
We could optimize the previous method in two ways. We can delay
the multiplication by body.mass until the end (since it can be
factored in each individual force). We can even forget completely
about body.mass and change Body:applyForce so that it does not
divide the force by body.mass (as the multiplication and division
cancel each other). As it makes the physics less clear in both cases,
I decided to avoid these optimizations.
The next method is the simulation main loop: for each body, compute
the force and apply it to the body. (The astute reader will notice
that applying the force changes the speed of the body, not its
position. The position will remain unchanged for the computation of
other forces and the speed is not used to compute the forces, so the
order in which computations occur does not matter.)
A second loop simulates the bodies moving, changing their position
according to their speed.
function Bodies:simulate(time)
local force
for i,b in ipairs(self.blist) do
force = self:computeForce(i,b)
b:applyForce(force,time)
end
for i,b in ipairs(self.blist) do
b:updatePos(time)
end
end
The last two methods are for the graphical interface. The first one
searches whether there is a body close (less than 20 pixels) to a
given screen point.
function Bodies:bodyAtPosition(touch,zoom,center)
for i,b in ipairs(self.blist) do
if touch:distSqr(b:screenPosition(zoom,center))
< 400 then
return i
end
end
return 0
end
The last function draws everything.
function Bodies:draw(z,c)
for i,b in ipairs(self.blist) do
b:draw(z,c)
end
end
We can now turn to the Main class, which will let us look deeper
into Codea itself.
We first set some parameters. They appear as sliders on the side of
the screen. The first one controls the speed of the simulation (how
much time passes for each simulation step), and the second how much
zoomed in we are. (I should use the usual zooming gesture for that,
but I haven’t taken the time to do so.)
parameter("Speed",0,1000,10)
parameter("Zoom",0,1000,10)
We next set some constants. I’ll talk more about this below, but in a
nutshell we are choosing the mass of the Earth as mass unit, and 1000
km as distance unit.
SunMassInEarth = 330000
EarthMassInMoon = 81
EarthMoon = 384 -- thousands of kilometers
EarthSun = 149600 -- thousands of kilometers
The setup() function is called once when the program is launched. We
use it to add three bodies (the Sun, the Earth, the Moon) at some
initial positions (so that distances between them are as expected) and
some initial velocity I tweaked to get the orbits to look right (more
about this below).
-- Use this function to perform your initial setup
function setup()
bl=Bodies()
bl:addBody(Body(0,0,0,0,SunMassInEarth,
color(255,255,0,255)))
bl:addBody(Body(0,EarthSun,1.5,0,1,
color(0,255,255,255)))
bl:addBody(Body(EarthMoon,EarthSun,1.5,-0.05,
1/EarthMassInMoon,
color(255,255,255,255)))
center = vec2(0,0)
zfactor = 10000
centeredOnBody = 0
print("Tap on a body to track it")
print("Tap anywhere else to stop tracking it")
print("Drag to move to view point")
end
The touched(touch) function is called every time a touch event
occurs. If the touch has just began, then we check if there is a body
near. If so, we remember to track this body (to set the center of the
display on this body). If the touch is MOVING, it means that the
user is dragging his finger, and we move the center of the display
accordingly.
function touched(touch)
if touch.state == BEGAN then
centeredOnBody = bl:bodyAtPosition(
vec2(touch.x,touch.y),
Zoom/zfactor,center)
elseif touch.state == MOVING then
center.x = center.x +
(lastx - touch.x)*zfactor/Zoom
center.y = center.y +
(lasty - touch.y)*zfactor/Zoom
centeredOnBody = 0
end
lastx = touch.x
lasty = touch.y
end
Finally, the draw() function is automatically called every frame
(this is the main loop of the program). There we erase everything, and
if the simulation is not paused, we make it progress for Speed
amount of time. If we’re tracking a body, we set the center of the
display to this body. Finally, we redraw everything.
-- This function gets called once every frame
function draw()
background(0)
if Speed > 0 then bl:simulate(Speed) end
if (centeredOnBody ~= 0) then
local b = (bl:bodyAtIndex(centeredOnBody))
center.x = b.pos.x
center.y = b.pos.y
end
bl:draw(Zoom/zfactor,center)
end
I could say “that’s it”, but there is one more thing…
On units of measurement
This simulation tries to be realistic. But except for the mass and
distance of the bodies, there is no magic number. In particular, I
chose G to be equal to 1 instead of it’s real value (6.673
10^-11). How can this work?
It works because I never specify anywhere what the time unit is, and
mass, distance, and time are
independant base units. I
thus have 3 degrees of freedom to use, and by choosing a unit of mass
(weight of the Earth), a unit of distance (1000 km), and the constant
G (1), I fix the time unit to something we can compute. This is how
I did it using Soulver, keeping
probably too many digits…
| # |
Soulver Input |
Result |
|---|
| 1 |
In m3 kg-1 s-2, G: 6.673 × 10-11 |
6.673×10-11 |
| 2 |
In kg, Earth mass (EM): 5.974 × 1024 |
5.974×1024 |
| 3 |
In m3 EM-1 s-2, G: 6.673×1011 × 5.974×1024 |
3.987×1014 |
| 4 |
In (1000km)3 EM-1 s-2, G: 3.987×1014/((106)3) |
0.0003986 |
| 5 |
thus our time unit is: √(1/0.0003986) |
50.08 |
| 6 |
|
|
| 7 |
Earth speed: 2 × π × 149,6 × 106 × 103/(365 × 24 × 3600) |
29806 |
| 8 |
in our units: 29806 × 50.08/106 |
1.493 |
| 9 |
Moon speed: 2 × π × 384 × 103 × 103/(28 × 24 × 3600) |
997.3 |
| 10 |
In our units: 997.3 × 50.08/106 |
0.04995 |
As G is expressed in “something per kg”, to change it into
“something per Earth mass”, we multiply it by the mass of the
Earth (line 3). To go from meters to thousands of kilometers, we
divide the result by one million thrice (since G is in meters cube,
line 4). To conclude, let’s call foo our time units, which is a
number of seconds. To convert G in this unit, we need to multiply it
by foo twice (or by foo squared), and the result would be 1 (which
is what we chose for G):
line4 * foo^2 = 1
foo^2 = 1 / line4
foo = √(1 / line4)
which is what we compute on line 5. Our time unit is thus 50
seconds. If we call our simulation with a speed of 10, for each frame
we simulate 500 seconds.
The last few lines are checks that everything is consistent. In line
7, I compute the Earth orbital speed in meters per seconds. It’s
simply going around one orbit 2 × π × 149,6 × 10^6 × 10^3 (Earth is
150 million of kilometers away from the sun) in one year. When
converting this to our units (thousands of kilometers and fifties of
seconds), we get 1.49, which is quite close to the 1.5 I tweaked for
initial Earth velocity to get a circular-looking orbit. Similarly, we
can compute the orbital speed of the Moon, and find it’s consistent as
well.
Final words, or tldr
There are many things I want to add to this, such as trails for
trajectories, zooming using gestures, working in 3d (both for the
simulation, easy, and for the interface, much harder)… But with this
very simple program, I was able to get my son interested in celestial
mechanics and in the mind-blowing distances in the solar system.
This post was way longer than I expected, and probably too
technical. If there is one thing to take away, try the following:
Codea is a really cool app, and it’s quite easy and to play with
physical simulations in it. And physics, even as simple as this, is
fun!
