Makie Visualization for Planetary Simulation for High School Students

Hi,

I am going to teach some simple orbital dynamics to high school students in summer.
The students are going to implement the gravitational forces and some sort of time integration (simple explicit or implicit euler) by themselves.
Also, we want to include some rockets designs with simplified thrust models and weight loss. So that is mostly their learning part.

Is there a easy way to use Makie as some 3D viewer where one can maybe change the perspective and navigate in space?

Frankly I am a bit overwhelmed about all the options so I’d be very happy about any hints to visualize this.

Thank you :)),

Felix

1 Like

See GitHub - HolyLab/FlyThroughPaths.jl: Create "fly-throughs" in 3d visualization, you could either just use it directly or see how it manipulates the camera to do it yourself!

You can also set camera settings in LScene’s camera I think, via cameracontrols(scene) (or cameracontrols(lscene.scene)).

Also see the controls here:

1 Like

You can also try to make a small interactive plot with GLMakie like this.
Note that I just made a quick example so it is not consistent with astrodynamics.

using GLMakie

sun = Point3f(0.0,0.0,0.0)
planenet1(t) = Point3f(1.0*cos(0.02 *t), 1.2*sin(0.02*t), 0.0)
planenet2(t) = Point3f(2.3*cos(0.01 *t + 1), 2.1*sin(0.01*t + 1 ), 0.0)

let 
    fig = Figure()
    ax = Axis3(fig[1,1], 
        limits=((-3.0,3.0),(-3.0,3.0),nothing))

    sg = SliderGrid(fig[2, 1],(label = "t", range = 0:0.1:1000, startvalue = 0.0))

    p1 = Observable([planenet1(0)])
    p2 = Observable([planenet2(0)])
    p1_track = Observable(planenet1.(0:0))
    p2_track = Observable(planenet2.(0:0))

    scatter!(ax, sun, color=:orange, markersize=30)
    scatter!(ax, p1, color=:green)
    scatter!(ax, p2, color=:blue)
    
    lines!(ax, p1_track, color=:green)
    lines!(ax, p2_track, color=:blue)

    on(sg.sliders[1].value) do t
        p1[] = [planenet1(t)]
        p2[] = [planenet2(t)]

        p1_track[] = planenet1.(0:0.1:t)
        p2_track[] = planenet2.(0:0.1:t)
    end
    
    fig
end

5 Likes

This looks promising!

I give it a try :)!

Note that you can rotate the view with a left mouse click and drag on the figure. Use right click and drag to pan the view.

1 Like

Oh, it doesn’t seem to work with Pluto.
The plot window is entirely static and not interactive.

But just adding using WGLMakie fixed it!

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

2 Likes

Awesome! If you happen to want some more realism, then you can use mesh with a translation instead of scatter, and meshes can have textures. See Beautiful Makie (you can ignore the SSAO stuff). This was just hacking the function you wrote to use a planet_colors dictionary to get a color from a planet name, and then using mesh with a sphere of size 1e10 to render. Unfortunately it’s a bit complex to do this in pixel space, but if that’s something you’d like to do then we can certainly pull it off. You will need using GeometryBasics, FileIO for this. Makie already loads those so there’s no extra dependencies / cost.

Just replace your plotting loop with:


    for (planet_name, data) in planets
        mp = GLMakie.mesh!(ax, uv_normal_mesh(Sphere(Point3f(0.0, 0.0, 0.0), 1.0e10)); color = get(planet_colors, planet_name, :gray))
        on(mp, data[:observable]) do position
            translate!(mp, position)
        end
        GLMakie.lines!(ax, data[:track])
    end

where planet_colors is some kwarg to the function. Then,

planet_colors = Dict(
    "mars" => FileIO.load(download("https://www.solarsystemscope.com/textures/download/2k_mars.jpg")), 
    "sun" => FileIO.load(download("https://www.solarsystemscope.com/textures/download/2k_sun.jpg")), 
    "earth" => FileIO.load(download("https://www.solarsystemscope.com/textures/download/2k_earth_daymap.jpg")), 
    "moon" => :gray,
    "jupiter" => FileIO.load(download("https://www.solarsystemscope.com/textures/download/2k_jupiter.jpg")), 
    "mercury" => FileIO.load(download("https://www.solarsystemscope.com/textures/download/2k_mercury.jpg")), 
    "saturn" => FileIO.load(download("https://www.solarsystemscope.com/textures/download/2k_saturn.jpg")),
    "venus" => FileIO.load(download("https://www.solarsystemscope.com/textures/download/2k_venus_surface.jpg")),
)
1 Like