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
    r = math.rad(self.deg)
    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.
That being said, the answer to your questions are:
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. :wink:

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! :wink:

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)
                
                --adjust v velocity
                v.dx = v.dx - (v.x - v2.x) * d * v2.m
                v.dy = v.dy - (v.y - v2.y) * d * v2.m
                
                --adjust v2 velocity
                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. :slight_smile: 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)
                
                --adjust v velocity
                v.dx = v.dx - (v.x - v2.x) * d * v2.m * G
                v.dy = v.dy - (v.y - v2.y) * d * v2.m * G
                
                --adjust v2 velocity
                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