Orbit integrator.lua

-- 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