Celestial mechanics

This is a very simple Sun Earth Moon simulation, which can easily be extended with more bodies.

The tricky part was to get the numbers right to have the moon rotate around the earth.

https://gist.github.com/1396195

There was an embarrassing bug (forgot the r-square in the gravitation force). This is fixed here: https://gist.github.com/1397257

That’s a really cool simulation - though it’s quite hard to see it all moving at once given the distance of the sun.

It would be kind of cool to see the trajectory of an asteroid passing through.

Showing trajectories, or at least trails of them, would be nice. I need to think about it.

In the meantime, here is a write-up about this simulation, it may contain information of interest for people starting with Codea.

@brab: Your write-up works very well as a tutorial (read: ambrosia for us neophytes). Would you consider adding it to the wiki please? Or adding a link to the page on your own site?

Thanks!

I added a link on the GettingStarted page. Don’t hesitate to let me know if you think it’s not appropriate.

Added under User Contirbutions as well https://bitbucket.org/TwoLivesLeft/codea/wiki/UserContrib

I just saw this - very nice. In college we had to do this as an exercise in a Physics class, but no graphics, just listing coordinates through time. Most people did it in Fortran(!) but I used APL, and it was very concise. The prof didn’t know APL and had to take it on faith that it actually worked.

I did something similar here:
http://twolivesleft.com/Codea/Talk/discussion/comment/17036#Comment_17036

Here is something I wrote awhile back, but I don’t know if I posted it or not. I couldn’t find it in the forum. This is based on code I posted previously back then. This shows 2 stars orbiting each other with a 3rd star orbiting the 2. This will continue for quite awhile before they start to break down and one of them is thrown out. You can move them up or down on the screen to keep them in view. A double tap on the screen will turn the trail mode on or off.


displayMode(FULLSCREEN)
supportedOrientations(PORTRAIT)

function setup()
    show=0
    yy=0
    fill(255)
    
    -- star 1
    gs1 = 400      -- gravity of star 1
    xs1 = 300      -- start x position
    ys1 = 500      -- start y position
    dxs1 = 0       -- initial change in x position
    dys1 = 2       -- initial change in y position
   
    -- star 2 
    gs2 = 400      -- gravity of star 2
    xs2 = 350      -- start x position
    ys2 = 500      -- start y position
    dxs2 = 0       -- initial change of x position
    dys2 = -2      -- initial change in y position
    
    -- star 3
    gs3 = 400      -- gravity of star 3
    xs3 = 600      -- start x position
    ys3 = 500      -- start y position
    dxs3 = 0       -- initial change of x position
    dys3 = 1.45    -- initial change in y position
end

function touched(t)
    if t.state==BEGAN and t.tapCount==2 then
        if show==1 then
            show=0
        else
            show=2
        end
    end
    if t.state==MOVING then
        yy=yy-t.deltaY
        background()
    end
end

function draw()
    if show==2 then
        backingMode(RETAINED)
        size=2
        show=1
    elseif show==0 then
        background()
        size=8
    end
        
    -- gravity star 2  to star 1
    d = gs2 / math.pow(( (xs1 - xs2) * (xs1 - xs2) + (ys1 - ys2) * (ys1 - ys2) ), 1.5 )
    dxs1 = dxs1 - (( xs1 - xs2 ) * d)
    dys1 = dys1 - (( ys1 - ys2 ) * d)

    -- gravity star 1  to star 2
    d = gs1 / math.pow(( (xs2 - xs1) * (xs2 - xs1) + (ys2 - ys1) * (ys2 - ys1) ), 1.5 )
    dxs2 = dxs2 - (( xs2 - xs1 ) * d)
    dys2 = dys2 - (( ys2 - ys1 ) * d)
    
    -- gravity star 2  to star 3
    d = gs2 / math.pow(( (xs3 - xs2) * (xs3 - xs2) + (ys3 - ys2) * (ys3 - ys2) ), 1.5 )
    dxs3 = dxs3 - (( xs3 - xs2 ) * d)
    dys3 = dys3 - (( ys3 - ys2 ) * d)
    
    -- gravity star 3  to star 2
    d = gs3 / math.pow(( (xs2 - xs3) * (xs2 - xs3) + (ys2 - ys3) * (ys2 - ys3) ), 1.5 )
    dxs2 = dxs2 - (( xs2 - xs3 ) * d)
    dys2 = dys2 - (( ys2 - ys3 ) * d)
    
    -- gravity star 1  to star 3
    d = gs1 / math.pow(( (xs3 - xs1) * (xs3 - xs1) + (ys3 - ys1) * (ys3 - ys1) ), 1.5 )
    dxs3 = dxs3 - (( xs3 - xs1 ) * d)
    dys3 = dys3 - (( ys3 - ys1 ) * d)
    
    -- gravity star3  to star 1
    d = gs3 / math.pow(( (xs1 - xs3) * (xs1 - xs3) + (ys1 - ys3) * (ys1 - ys3) ), 1.5 )
    dxs1 = dxs1 - (( xs1 - xs3 ) * d)
    dys1 = dys1 - (( ys1 - ys3 ) * d)
    
    xs1 = xs1 + dxs1
    ys1 = ys1 + dys1
    xs2 = xs2 + dxs2
    ys2 = ys2 + dys2
    xs3 = xs3 + dxs3
    ys3 = ys3 + dys3
       
    ellipse(xs1,ys1-yy,size)
    ellipse(xs2,ys2-yy,size)
    ellipse(xs3,ys3-yy,size)
end