# Using physics to simulate moon in orbit around planet?

I’m still in “practice” mode since I’m new to Codea and am just exploring what I can do. It’s been great fun so far! I have been looking at the physics functions and, specifically, gravity. Is it possible to apply gravity to an object? I’ve attempted to define a circular physics body and then would like to apply gravity to it, but I’m not sure it’s possible. I envision something along the lines of

``````planet = physics.body(CIRCLE, 100)
planet.type = STATIC
planet.gravity(1000)
``````

Thanks for any tips!

I think the 2d physics thing isn’t really built for this. I guess if you turned off gravity altogether so that CIRCLEs would just sit there, then you could calculate your planetary gravity yourself and use body.applyForce( force ) ( planet.applyForce in your example) to let physics do the rest of the work.

I’m not entirely familiar with applyForce() but my belief is that it applies a force to the specified body. If this is the case, then wouldn’t I have to apply the force to the non-planet objects and essentially “push” them towards the planet?

You are right, I haven’t used applyForce myself either, it just looks right for the task.

Of course notionally you need to calculate the gravity force and apply it to both the satelite and the planet, although the planet should be heavy enough that that is fairly irrelevant…

might be able to do it with distance joint as well, if you just want it to rotate the planet. and then use a force to the moon.

If you’re interested in a non-physics solution, this is the heart of a little orbiting planets routine from the game orb-bits.

``````function Orbiter:draw()
local r
noFill()
stroke(self.clr)
strokeWidth(2)
-- the further from the center, the slower it goes
self.deg = self.deg + (100 / self.radius * -1)
-- don't really have to, but for human readability's sake
if self.deg > 360 then
self.deg = self.deg - 360
end
-- position in orbit
self.x = self.radius * math.cos(r) + (WIDTH / 2)
self.y = self.radius  * math.sin(r) + (HEIGHT / 2)
ellipse(self.x, self.y, self.size)
end
``````

Hi @mcarthey You won’t be able to do what you want until Codea is changed to allow gravity to be towards a point instead of towards the bottom of the screen. To calculate actual orbits, here is an example I wrote that simulates 2 stars orbiting each other. If you comment out the 2 lines of code after the comment “new position of star 2”, it will simulate a moon orbiting a planet. You can change the different values in setup to cause different orbits. I’m not sure if this will interest you, so I just worked on the calculations and not much in the graphics area.

``````
displayMode(FULLSCREEN)
supportedOrientations(PORTRAIT)

function setup()

-- star 1
g_s1 = 400      -- gravity of star 1
x_s1 = 300      -- start x position
y_s1 = 500      -- start y position
dx_s1 = 1.4     -- initial change in x position
dy_s1 = .4      -- initial change in y position

-- star 2
g_s2 = 400      -- gravity of star 2
x_s2 = 350      -- start x position
y_s2 = 600      -- start y position
dx_s2 = -1.4    -- initial change of x position
dy_s2 = -.4     -- initial change in y position

end

function draw()
background(40,40,50)
fill(255)

-- gravity pull of star 2 on star 1
d = g_s2 / math.pow(( (x_s1 - x_s2) * (x_s1 - x_s2) + (y_s1 - y_s2) * (y_s1 - y_s2) ), 1.5 )

-- calculate dx, dy of star 1 towards star 2
dx_s1 = dx_s1 - (( x_s1 - x_s2 ) * d)
dy_s1 = dy_s1 - (( y_s1 - y_s2 ) * d)

-- gravity pull of star 1 on star 2
d = g_s1 / math.pow(( (x_s2 - x_s1) * (x_s2 - x_s1) + (y_s2 - y_s1) * (y_s2 - y_s1) ), 1.5 )

-- calculate dx, dy of star 2 towards star 1
dx_s2 = dx_s2 - (( x_s2 - x_s1 ) * d)
dy_s2 = dy_s2 - (( y_s2 - y_s1 ) * d)

-- new position of star 1
x_s1 = x_s1 + dx_s1
y_s1 = y_s1 + dy_s1

-- new position of star 2
-- comment out the next 2 lines to keep star 2 in a fixed position
x_s2 = x_s2 + dx_s2
y_s2 = y_s2 + dy_s2

-- plot star 1
ellipse(x_s1,y_s1,10)

-- plot star 2
ellipse(x_s2,y_s2,10)
end

``````

It’s looks great and since most of the work involved in any graphics project is the math, I greatly appreciate you sharing this. Thanks!

I wish I could say I understand all the math, but from your comments it appears that the ellipses are actually attracted and circling due to the gravity calculations? Very impressive. Is this simply because their attraction forces are equal or is there more to it?
OK, well, slightly stupid question I know since if it was just gravity they would simply stick to each other. Heh. What is it that keeps them circling?

This is the ‘universal law of gravity’ as was discovered by Newton: a force inversely proportional to the square distance mathematically generates ellipses. The gravity force acts a little bit like a string: if you tie a small stone to a string and spin it around yourself, the stone is following an ellipse (a circle is a particular case of an ellipse) and, once started, the sole force you apply to the stone is pulling towards youself through the string.
they would simply stick to each other.
That would be the case if you and the stone are initially motionless: pulling the rope will make the stone stick to you. The movement is a segment, which is a very flat ellipse.
Heh. What is it that keeps them circling?.
That’s their initial speed: once stone is launched fast enough with a speed not towards you, pulling the rope only curves the direction of the stone. Of course on earth the normal gravity and the air friction makes the movement more complex to sustain, so the example of the rope has some limits, but that’s the idea.

Hi @mcarthey You’re correct about the gravity calculations. I’m constantly calculating the speed in 2 directions. One direction is star 1 towards star 2 due to gravity, and the other is the forward direction of star 1 around star 2. Think of it as throwing a baseball level to the ground. Gravity pulls the ball down while the speed of the ball moves it forward. The harder you throw the ball the farther it goes before it hits the ground. If you could throw it fast enough and ignoring air resistance, the curvature of the earth would drop as fast as the ball was falling so the ball wouldn’t hit the ground. I’m not sure if these calculations are what you’re after because they’re unstable due to the precision of the Codea math. This is just a simple simulation, so I wasn’t concerned about the precision. Here is another simple variation of the program. It’s a moon orbiting a planet that allows you to increase or decrease the forward speed of the moon. Each time you tap “faster” or “slower”, the speed changes. You have to tap each time for a change, not just tap and hold. You can change the orbit size by increasing or decreasing the speed.

``````
displayMode(FULLSCREEN)
supportedOrientations(PORTRAIT)

function setup()
inc=1
-- moon
x_s1 = 375      -- start x position
y_s1 = 400      -- start y position
dx_s1 = 1       -- initial change in x position
dy_s1 = 0       -- initial change in y position

-- planet
g_s2 = 400      -- gravity of star 2
x_s2 = 375      -- x position
y_s2 = 600      -- y position
end

function touched(t)
if t.state==BEGAN then
if (t.y-950)^2/50^2+(t.x-600)^2/100^2 <= 1 then
inc=1.05
colorF=true
end
if (t.y-850)^2/50^2+(t.x-600)^2/100^2 <= 1 then
inc=.95
colorS=true
end
end
if t.state==ENDED then
colorF=false
colorS=false
end
end

function draw()
background(40,40,50)
speed()

dx_s1 = dx_s1*inc
dy_s1 = dy_s1*inc
inc=1

d = g_s2 / math.pow(( (x_s1 - x_s2) * (x_s1 - x_s2) + (y_s1 - y_s2) * (y_s1 - y_s2) ), 1.5 )
dx_s1 = dx_s1 - (( x_s1 - x_s2 ) * d)    -- calc dx,dy of moon
dy_s1 = dy_s1 - (( y_s1 - y_s2 ) * d)

x_s1 = x_s1 + dx_s1  -- new position of moon
y_s1 = y_s1 + dy_s1

fill(255)
ellipse(x_s1,y_s1,5)    -- plot moon
ellipse(x_s2,y_s2,20)   -- plot planet
end

function speed()
fill(255)
if colorF then
fill(255,0,0)
end
ellipse(600,950,200,50)

fill(255)
if colorS then
fill(255,0,0)
end
ellipse(600,850,200,50)

fill(255, 0, 0, 255)
text("Faster",600,950)
text("Slower",600,850)
end

``````

Thanks @Jmv38 and @dave1707 for the detailed explanations. I think I’m going to need a second cup of coffee.

Coffee is a brewed beverage with a distinct aroma and flavor prepared from the roasted seeds of the Coffea plant. The seeds are found in coffee “cherries”,…
No, just joking!

I took dave1707’s solution and refactored it a bit. You can now have as many celestial bodies as you like. Also, touching within 100 pixels of a body will select it, and then you can hit slower/faster to adjust it’s speed by 10%

Initial setup is 2 planets and a moon, playing with faster/slower is a good way to lose a body… You can add more planets/etc by adding bodies[x] = lines in setup.

``````displayMode(FULLSCREEN)
--supportedOrientations(PORTRAIT)

function setup()

bodies = {}

--initial planets
bodies[1] = { m = 400, x= 300, y= 500, dx = 1.4, dy = .4 }
bodies[2] = { m = 400, x= 300, y= 600, dx = -1.4, dy = -.4 }
bodies[3] = { m = 4, x = 300, y= 480, dx = -3.1, dy = .4 }

selectedBody = 1
end

function touched(touch)
if touch.state == BEGAN then
if touch.x > WIDTH - 100 and touch.y < HEIGHT - 85 and touch.y > HEIGHT - 125 then
--touched slower
bodies[selectedBody].dx = bodies[selectedBody].dx * .9
bodies[selectedBody].dy = bodies[selectedBody].dy * .9
elseif touch.x > WIDTH - 100 and touch.y < HEIGHT - 140 and touch.y > HEIGHT - 180 then
--touched faster
bodies[selectedBody].dx = bodies[selectedBody].dx * 1.1
bodies[selectedBody].dy = bodies[selectedBody].dy * 1.1
else
--change slected planet if touch is with 100 of it
closest = nil
dist = 100
for k,v in pairs(bodies) do
calcDist = math.sqrt((touch.x - v.x)^2+(touch.y - v.y)^2)
if  calcDist < dist then
dist = calcDist
closest = k
end
end
if closest ~= nil then
selectedBody = closest
end
end
end
end

function draw()
background(40,40,50)
strokeWidth(0)
fill(255)

applyGravity()

moveBodies()

-- plot bodies
for k,v in pairs(bodies) do
if selectedBody == k then
fill(255, 9, 0, 255)
else
fill(255,255,255,255)
end
ellipse(v.x, v.y, v.m/40+2)
end

drawControl()
end

function drawControl()
fill(255,255,255,255)
text(string.format("Body: %d", selectedBody), WIDTH - 50, HEIGHT - 50)
fill(150,150,150,255)
rect(WIDTH - 100, HEIGHT - 125, 100, 40)
rect(WIDTH - 100, HEIGHT - 180, 100, 40)
fill(255,255,255,255)
text("slower", WIDTH - 50, HEIGHT - 105)
text("faster", WIDTH - 50, HEIGHT - 160)
v = math.sqrt(bodies[selectedBody].dx^2+bodies[selectedBody].dy^2)
text(string.format("Velocity: %f", v), WIDTH - 100, HEIGHT - 200)
end

function applyGravity()
for k,v in pairs(bodies) do
for k2,v2 in pairs(bodies) do
if k2 > k then
--gravity factor - this is distance facotring for gravity and trigonometry
d = 1 / math.pow(( (v.x - v2.x)^2 + (v.y - v2.y)^2), 1.5)

v.dx = v.dx - (v.x - v2.x) * d * v2.m
v.dy = v.dy - (v.y - v2.y) * d * v2.m

v2.dx = v2.dx - (v2.x - v.x) * d * v.m
v2.dy = v2.dy - (v2.y - v.y) * d * v.m
end
end
end
end

function moveBodies()
for k,v in pairs(bodies) do
v.x = v.x + v.dx
v.y = v.y + v.dy
end
end
``````

Thanks everyone. I am quickly learning lots of things, not the least of which is that I have a lot of learning to do. The way you guys whip out this complicated code is very impressive and I greatly appreciate it. I have some ideas for apps and will be working to include much of what I’m learning here. Great stuff, and thanks again!

One last one. This is tweaked a little to allow for scaling by adjusting the gravitational constant, and some fudging of scaling for the ellipses… mainly so that I could model a somewhat stable solar system.

I’ve left both planetary with a variable called setup - if you change it to 2 it’ll be the 2 planet 1 moon as previous

Also, if you change leaveTrails to true it will leave trails as the bodies move which will allow you to fine tune your orbits perhaps…

``````
displayMode(FULLSCREEN)
--supportedOrientations(PORTRAIT)

function setup()
leaveTrails = false
if leaveTrails then
backingMode(RETAINED)
end
bodies = {}

setup = 1

if setup == 1 then
--solar system
--gravitation constant for scaling
G = .005
--initial planets
--sun
bodies[1] = { m= 40000, x = WIDTH/2, y= HEIGHT/2, dx=0, dy=0 }
--planets
bodies[2] = { m= 15, x = WIDTH/2, y= HEIGHT/2 - 65, dx= 1.8, dy = 0 }
bodies[3] = { m= 20, x = WIDTH/2, y= HEIGHT/2 - 100, dx=1.5, dy=0 }
bodies[4] = { m= 30, x = WIDTH/2, y= HEIGHT/2 - 150, dx=1.3, dy=0 }
bodies[5] = { m= 40, x = WIDTH/2, y= HEIGHT/2 - 180, dx=1.2, dy= 0}
bodies[6] = { m= 60, x = WIDTH/2, y= HEIGHT/2 - 300, dx=0.85, dy= 0 }
--moons
bodies[7] = { m= 0.1, x = WIDTH/2, y= HEIGHT/2 - 155, dx= 1.47, dy= 0 }
end

if setup == 2 then
--2 planets and a moon
G = 1
--initial planets
bodies[1] = { m = 400, x= 300, y= 500, dx = 1.4, dy = .4 }
bodies[2] = { m = 400, x= 300, y= 600, dx = -1.4, dy = -.4 }
bodies[3] = { m = 4, x = 300, y= 480, dx = -3.1, dy = .4 }
end

selectedBody = 1
end

function touched(touch)
if touch.state == BEGAN then
if touch.x > WIDTH - 100 and touch.y < HEIGHT - 85 and touch.y > HEIGHT - 125 then
--touched slower
bodies[selectedBody].dx = bodies[selectedBody].dx * .9
bodies[selectedBody].dy = bodies[selectedBody].dy * .9
elseif touch.x > WIDTH - 100 and touch.y < HEIGHT - 140 and touch.y > HEIGHT - 180 then
--touched faster
bodies[selectedBody].dx = bodies[selectedBody].dx * 1.1
bodies[selectedBody].dy = bodies[selectedBody].dy * 1.1
else
--change slected planet if touch is with 100 of it
closest = nil
dist = 100
for k,v in pairs(bodies) do
calcDist = math.sqrt((touch.x - v.x)^2+(touch.y - v.y)^2)
if  calcDist < dist then
dist = calcDist
closest = k
end
end
if closest ~= nil then
selectedBody = closest
end
end
end
end

function draw()
if leaveTrails == false then
background(40,40,50)
end
strokeWidth(0)
fill(255)

applyGravity()

moveBodies()

-- plot bodies
for k,v in pairs(bodies) do
if selectedBody == k then
fill(255, 9, 0, 255)
else
fill(255,255,255,255)
end
ellipse(v.x, v.y, math.pow(v.m, .2)*2 + 2)
end

drawControl()
end

function drawControl()
fill(255,255,255,255)
text(string.format("Body: %d", selectedBody), WIDTH - 50, HEIGHT - 50)
fill(150,150,150,255)
rect(WIDTH - 100, HEIGHT - 125, 100, 40)
rect(WIDTH - 100, HEIGHT - 180, 100, 40)
fill(255,255,255,255)
text("slower", WIDTH - 50, HEIGHT - 105)
text("faster", WIDTH - 50, HEIGHT - 160)
v = math.sqrt(bodies[selectedBody].dx^2+bodies[selectedBody].dy^2)
text(string.format("Velocity: %f", v), WIDTH - 100, HEIGHT - 200)
end

function applyGravity()
for k,v in pairs(bodies) do
for k2,v2 in pairs(bodies) do
if k2 > k then
--gravity factor - this is distance facotring for gravity and trigonometry
d = 1 / math.pow(( (v.x - v2.x)^2 + (v.y - v2.y)^2), 1.5)

v.dx = v.dx - (v.x - v2.x) * d * v2.m * G
v.dy = v.dy - (v.y - v2.y) * d * v2.m * G

v2.dx = v2.dx - (v2.x - v.x) * d * v.m * G
v2.dy = v2.dy - (v2.y - v.y) * d * v.m * G
end
end
end
end

function moveBodies()
for k,v in pairs(bodies) do
v.x = v.x + v.dx
v.y = v.y + v.dy
end
end
``````

@spacemonkey Nice change. I was able to get the moon to alternate orbits from one planet to the other a few times before flying off the screen. I found it hard to select the moon, so I changed your code so if I touched the screen away from the buttons, it cycled thru each object to select it. That was easier than trying to tap on a moving object.

Not sure I recommend using `setup` as a variable name …

Is there some sort of ‘Physics Equations for Programmers Who Slept Through Class’ book? Or maybe a site that would summarize this sort of thing? The logic of the code isn’t terribly complex but it’s very intimidating to go through this sort of thing from scratch without any understanding of where the formulas come from.

The challenge is not helped because myself and dave1707 have done some optimising, compressing multiple equations into 1 to reduce the amount of calculation being done.

Below is a deoptimised piece of code, I took out the touch and side display stuff and elaborated with comments on the applyGravity section.

Hope this helps

``````

displayMode(FULLSCREEN)
--supportedOrientations(PORTRAIT)

function setup()
leaveTrails = false
if leaveTrails then
backingMode(RETAINED)
end
bodies = {}

planetarySetup = 1

if planetarySetup == 1 then
--solar system
--gravitation constant for scaling
G = .005
--initial planets
--sun
bodies[1] = { m= 40000, x = WIDTH/2, y= HEIGHT/2, dx=0, dy=0 }
--planets
bodies[2] = { m= 15, x = WIDTH/2, y= HEIGHT/2 - 65, dx= 1.8, dy = 0 }
bodies[3] = { m= 20, x = WIDTH/2, y= HEIGHT/2 - 100, dx=1.5, dy=0 }
bodies[4] = { m= 30, x = WIDTH/2, y= HEIGHT/2 - 150, dx=1.3, dy=0 }
bodies[5] = { m= 40, x = WIDTH/2, y= HEIGHT/2 - 180, dx=1.2, dy= 0}
bodies[6] = { m= 60, x = WIDTH/2, y= HEIGHT/2 - 300, dx=0.85, dy= 0 }
--moons
bodies[7] = { m= 0.1, x = WIDTH/2, y= HEIGHT/2 - 155, dx= 1.465, dy= 0 }
end

if planetarySetup == 2 then
--2 planets and a moon
G = 1
--initial planets
bodies[1] = { m = 400, x= 300, y= 500, dx = 1.4, dy = .4 }
bodies[2] = { m = 400, x= 300, y= 600, dx = -1.4, dy = -.4 }
bodies[3] = { m = 4, x = 300, y= 480, dx = -3.1, dy = .4 }
end

selectedBody = 1
end

function draw()
if leaveTrails == false then
background(40,40,50)
end
strokeWidth(0)
fill(255)

applyGravity()

moveBodies()

-- plot bodies
for k,v in pairs(bodies) do
if selectedBody == k then
fill(255, 9, 0, 255)
else
fill(255,255,255,255)
end
ellipse(v.x, v.y, math.pow(v.m, .2)*2 + 2)
end
end

function applyGravity()
--loop through the planets
for k,v in pairs(bodies) do
--loop through them again, we will apply the gravity from all planets in loop two to the current planet for loop one
for k2,v2 in pairs(bodies) do
--don't apply gravity to itself
if k2 ~= k then
--determine the distance between the planets (pythagorus theorem)
distance = math.sqrt((v.x - v2.x)^2 + (v.y - v2.y)^2)

--determine force of gravity (newton F = G*(m1*m2)/r^2
gravityForce = G*(v.m*v2.m)/(distance^2)

--determine how much acceleration this will have on planet 1 (F = m*a -> a = F/m)
acceleration = gravityForce/v.m

--now adjust velocity, but to convert this into the right vector use trigonometry
--sohcahtoa -> for x is cah (cosine = adjacent/hypotenuse)
-- for y is soh (sin = opposit/hypotenuse)
xCos = (v.x - v2.x)/distance
ySin = (v.y - v2.y)/distance

--apply the acceleration to the current vector adjusted by the factor for direction
--reverse trigonometry -> x (adjacent) = hypotenuse * Cos
-- y (opposite) = hypotenuse * Sin
v.dx = v.dx - (acceleration * xCos)
v.dy = v.dy - (acceleration * ySin)
end
end
end
end

function moveBodies()
for k,v in pairs(bodies) do
v.x = v.x + v.dx
v.y = v.y + v.dy
end
end
``````