Just for the sake of completeness.
With simple forward Euler (I know there is better ways, this is educational) we can achieve reasonable results:
# ╔═╡ d2c56190-fde5-11ef-01aa-e9d57cbfb029
using GLMakie, Makie, WGLMakie,LinearAlgebra
# ╔═╡ c388060c-3a0d-43a7-9718-88aeee6358b9
function Fg(r1, r2, m₁, m₂)
G = 6.67430 * 10^(-11)
F = G * m₁ * m₂ * (r2 .- r1) / sqrt(sum(abs2, r2 .- r1))^3
return F
end
# ╔═╡ 2f0d4d40-722a-421b-aaf2-ce6eab2c484c
struct Planet{T}
name::String
mass::T
position::Vector{T}
velocity::Vector{T}
end
# ╔═╡ 8c975dc2-c66c-4fd0-8b20-9b7947d98bb4
begin
earth = Planet("earth", 5.9722e24, [0.0, 152.10e9, 0.0], [29.29e3, 0.0, 0.0])
sun = Planet("sun", 1.9885e30, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0])
mars = Planet("mars", 6.4171e23, [0.0, 206.650e9 * cosd(1.850), 206.650e9 * sind(1.850)],
[26.50e3 * cosd(1.850), 0.0, 26.50e3 * sind(1.850)])
jupiter = Planet("jupiter", 1.8982e27, [0.0, 816.62e9 * cosd(1.303), 816.62e9 * sind(1.303)],
[13.72e3 * cosd(1.303), 0.0, 13.72e3 * sind(1.303)])
venus = Planet("venus", 4.8675e24, [0.0, 108.94e9 * cosd(3.394), 108.94e9 * sind(3.394)],
[35.02e3 * cosd(3.394), 0.0, 35.02e3 * sind(3.394)])
mercury = Planet("mercury", 3.3011e23, [0.0, 69.816e9 * cosd(7.004), 69.816e9 * sind(7.004)],
[38.86e3 * cosd(7.004), 0.0, 38.86e3 * sind(7.004)])
saturn = Planet("saturn", 5.6834e26, [0.0, 1_514.50e9 * cosd(2.485), 1_514.50e9 * sind(2.485)],
[10.18e3 * cosd(2.485), 0.0, 10.18e3 * sind(2.485)])
end
# ╔═╡ 8f7d9716-401f-4e6e-80d2-43e4e13a8b34
function simulate(planets::Vector{Planet{T}}, times) where T
# initialize logs
logs = Dict{String, Vector{Planet{T}}}()
for p in planets
logs[p.name] = [p]
end
Δt = times[2] - times[1]
for t in times
for p1name in [p.name for p in planets]
p1 = logs[p1name][end]
F = [0.0, 0.0, 0.0]
for p2name in [p.name for p in planets]
if p1name != p2name
p2 = logs[p2name][end]
F .+= Fg(p1.position, p2.position, p1.mass, p2.mass)
end
end
v = p1.velocity + F / p1.mass * Δt
r = p1.position + v * Δt
push!(logs[p1.name], Planet(p1.name, p1.mass, r, v))
end
end
return logs
end
# ╔═╡ c30c82c3-dc3d-4b67-aa18-ad0857e6652f
times = range(0, 86_400 * 10_000, 10_000)
# ╔═╡ a46a7d23-5e0d-4d6e-a3c9-607e6003f444
@time logs = simulate([earth, sun, mars, mercury, saturn, jupiter], times)
# ╔═╡ aa768373-e687-4aec-808a-b7d716f2bcc8
plot([l.position[1] for l in logs["earth"]], [l.position[2] for l in logs["earth"]])
# ╔═╡ aadafd7e-d148-488d-a22c-cbca61111a3b
# Define a function to create planet observables and plots
function plot_planet_trajectory(logs, planet_names, times_ints)
# Initialize the plot
fig = Figure()
ax = Axis3(fig[1, 1], limits=((-1000e9, 1000e9), (-1000e9, 1000e9), (-200e9, 200e9)))
sg = GLMakie.SliderGrid(fig[2, 1], (label="t", range=times_ints, startvalue=10))
# Initialize observables and data for each planet
planets = Dict()
for planet_name in planet_names
m_log = [Point3f(logs[planet_name][i].position...) for i in times_ints]
planets[planet_name] = Dict(
:m_log => m_log,
:observable => Observable(m_log[1]),
:track => Observable([m_log[1]])
)
end
# Plot the initial positions and tracks for each planet
#colors = Dict("earth" => :blue, "sun" => :yellow, "mars" => :red, "moon" => :gray)
for (planet_name, data) in planets
GLMakie.scatter!(ax, data[:observable])
GLMakie.lines!(ax, data[:track])
end
# Update function when slider is moved
on(sg.sliders[1].value) do t
for (planet_name, data) in planets
data[:observable][] = data[:m_log][t]
data[:track][] = @view data[:m_log][begin:t]
end
end
return fig
end
# ╔═╡ 35a54dcb-4204-49f6-abc2-9145ceadebe4
@time plot_planet_trajectory(logs, ["earth", "sun", "mars", "saturn", "mercury", "jupiter"], 1:length(times))
