public
anonymous / Orbit integrator.lua
Last active

  • Download Gist
Orbit integrator.lua
Lua
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
-- Notice: this is the simplest, but worst
-- possible integration routine.
 
sqrt = math.sqrt
 
AU = 1.5e11
G = 6.67e-11
Msun = 2e30
day = 86400
 
-- initial velocity
v0 = sqrt(G*Msun/AU)
-- initial position vector
x = { AU, 0 }
-- initial velocity vector
v = { 0, v0 }
 
t = 0
tout = 0
dt = 864
 
-- loop for 2 yrs
while (t < 2 * 365 * day) do
-- distance from center
r = sqrt(x[1]^2+x[2]^2)
vold = v
-- add acceleration to velocity vector
v[1] = v[1] + dt * (- G*Msun/r^2) * x[1]/r
v[2] = v[2] + dt * (- G*Msun/r^2) * x[2]/r
 
-- add velocity to position vector
x[1] = x[1] + dt * vold[1]
x[2] = x[2] + dt * vold[2]
 
t = t + dt
tout = tout + dt
-- every day, print out time and coordinates
if (tout > day) then
print(t, x[1], x[2], v[1], v[2])
tout = 0.
end
end

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.