Felipe Tavares' Avatar

Integrators: Predicting the Movement of the Moon

February 25, '21

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 RK4 integrator:

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][1])
        push!(orbits_y, p[i][2])
    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:

The Moon taking a full revolution around the Earth