In this post we are going to do some basic astrophysics. We are going to simulate the movements of the Moon around the Earth using Julia.
Lets warm up with some imports:
using JLD using Plots using LinearAlgebra
Now we will write two functions, one that given an array of all positions, velocities and masses returns the velocities at that point, so, trivially:
velocities(p, v, m) = v
and another which calculates the acceleration given the same inputs. For that we will need to iterate through all point masses and calculate the acceleration between those bodies, summing everything for each body.
function accelerations(p, v, m) local bodies = length(p) map(function (i) local sum = [0; 0] for j in 1:bodies if i != j local d = p[j]-p[i] sum += normalize(d)*(6.67408e-11*m[j]/norm(d)^2) end end sum end, 1:bodies) end
Do notice the gravitational constant tucked in there!
Now we create our integrator, in this case the
function integrate_rk4(p, v, m, Δt::Real) kp₁ = velocities(p, v, m) kv₁ = accelerations(p, v, m) kp₂ = velocities(p + kp₁*Δt/2, v + kv₁*Δt/2, m) kv₂ = accelerations(p + kp₁*Δt/2, v + kv₁*Δt/2, m) kp₃ = velocities(p + kp₂*Δt/2, v + kv₂*Δt/2, m) kv₃ = accelerations(p + kp₂*Δt/2, v + kv₂*Δt/2, m) kp₄ = velocities(p + kp₃*Δt, v + kv₃*Δt, m) kv₄ = accelerations(p + kp₃*Δt, v + kv₃*Δt, m) p+Δt/6*(kp₁+2*kp₂+2*kp₃+kp₄), v+Δt/6*(kv₁+2*kv₂+2*kv₃+kv₄), m end
Notice how each step uses the values from the previous step. Finally, lets define the initial masses, positions and velocities
# Earth at origin p = [ [0; 0], [0; 384402e3 ] ] # Earth still v = [ [0; 0], [1.022e3; 0] ] # Masses m = [ 5.972e24 , 7.348e22 ]
Finally, lets just call our integrator over and over, with time step of one hour, for 26.8 days, which is the time the Moon takes for one full revolution around the Earth.
orbits_x =  orbits_y =  for t ∈ 1:24*26.8 global p, v, m, orbits_x, orbits_y p, v, m = integrate_rk4(p, v, m, 3600) for i ∈ 1:length(p) push!(orbits_x, p[i]) push!(orbits_y, p[i]) end end plot( orbits_x, orbits_y, seriestype = :scatter, title = "Moon around the Earth", label = "Celestial Bodies Positions", aspect_ratio=:equal )
And behold, we have simulated, with physically accurate constants, the movements of the moon: