Makie.jl plots

lib/AstrodynamicalSolvers/docs/Project.toml

@@ -2,6 +2,7 @@
 AstrodynamicalCalculations = "c0cf9fb7-f496-4999-a425-c50785d1b88b"
 AstrodynamicalModels = "4282b555-f590-4262-b575-3e516e1493a7"
 AstrodynamicalSolvers = "636ee813-9c9e-4a0a-af07-88b3043dcb77"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"

lib/AstrodynamicalSolvers/docs/src/index.md

@@ -18,6 +18,8 @@ Circular Restricted Three Body Problem dynamics.
 This package contains differential correctors, and helpful wrapper functions, for 
 finding periodic orbits within Circular Restricted Three Body Problem dynamics.
 
+#### `Plots.jl`
+
 ```@example usage
 using AstrodynamicalSolvers
 using AstrodynamicalModels
@@ -45,11 +47,59 @@ end
 plot(planar, extraplanar, layout=(1,2))
 ```
 
+#### `Makie.jl`
+
+```@example
+using AstrodynamicalSolvers
+using AstrodynamicalModels
+using OrdinaryDiffEq
+using CairoMakie
+
+μ = 0.012150584395829193
+
+sol_planar = let
+    ic = halo(μ, 1) # lyapunov (planar) orbit
+    u = Vector(CartesianState(ic))
+    problem = ODEProblem(CR3BFunction(), u, (0, ic.Δt), (μ,))
+    solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14)
+end
+
+sol_extraplanar = let
+    ic = halo(μ, 2; amplitude=0.01) # halo (non-planar) orbit
+    u = Vector(CartesianState(ic))
+    problem = ODEProblem(CR3BFunction(), u, (0, ic.Δt), (μ,))
+    solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14)
+end
+
+fig = Figure(size=(800, 400); fontsize=11)
+
+ax_kwargs_common = (; aspect=:equal, azimuth=-π/3)
+
+ax_left = Axis3(fig[1, 1];
+    title = "Lyapunov Orbit",
+    limits = (0.78, 0.90, -0.09, 0.09, -0.02, 1.04),
+    ax_kwargs_common...,
+)
+ax_right = Axis3(fig[1, 2];
+    title = "Halo Orbit",
+    limits = (1.05, 1.26, -0.1, 0.1, -0.02, 0.01),
+    protrusions = (30, 100, 0, 0),
+    ax_kwargs_common...,
+)
+
+plot!(ax_left, sol_planar; idxs=(:x, :y, :z))
+plot!(ax_right, sol_extraplanar; idxs=(:x, :y, :z))
+
+fig
+```
+
 ### Manifold Computations
 
 Manifold computations, provided by `AstrodynamicalCalculations.jl`, can perturb 
 halo orbits onto their unstable or stable manifolds.
 
+#### `Plots.jl`
+
 ```@example
 using AstrodynamicalSolvers
 using AstrodynamicalCalculations
@@ -125,3 +175,98 @@ scatter!(figure, [1-μ], [0], label="Moon", xlabel="X (Earth-Moon Distance)", yl
 
 figure # hide
 ```
+
+#### `Makie.jl`
+
+```@example
+using AstrodynamicalSolvers
+using AstrodynamicalCalculations
+using AstrodynamicalModels
+using OrdinaryDiffEq
+using LinearAlgebra
+using CairoMakie
+
+μ = 0.012150584395829193
+
+unstable = let
+    ic = halo(μ, 1; amplitude=0.005)
+
+    u = CartesianState(ic)
+    Φ = monodromy(u, μ, ic.Δt, CR3BFunction(stm=true))
+
+    ics = let
+        problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, ic.Δt), (μ,))
+        solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(ic.Δt / 10))
+
+        solution.u
+    end
+
+    perturbations = [
+        diverge(ic[1:6], reshape(ic[7:end], 6, 6), Φ; eps=-1e-7)
+        for ic in ics
+    ]
+
+    problem = EnsembleProblem(
+        ODEProblem(CR3BFunction(), u, (0.0, 2 * ic.Δt), (μ,)),
+        prob_func=(prob, i, repeat) -> remake(prob; u0=perturbations[i]),
+    )
+
+    solution = solve(problem, Vern9(), trajectories=length(perturbations), reltol=1e-14, abstol=1e-14)
+end
+
+stable = let
+    ic = halo(μ, 2; amplitude=0.005)
+
+    u = CartesianState(ic)
+    Φ = monodromy(u, μ, ic.Δt, CR3BFunction(stm=true))
+
+    ics = let
+        problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, ic.Δt), (μ,))
+        solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(ic.Δt / 10))
+
+        solution.u
+    end
+
+    perturbations = [
+        converge(ic[1:6], reshape(ic[7:end], 6, 6), Φ; eps=1e-7)
+        for ic in ics
+    ]
+
+    problem = EnsembleProblem(
+        ODEProblem(CR3BFunction(), u, (0.0, -2.1 * ic.Δt), (μ,)),
+        prob_func=(prob, i, repeat) -> remake(prob; u0=perturbations[i]),
+    )
+
+    solution = solve(problem, Vern9(), trajectories=length(perturbations), reltol=1e-14, abstol=1e-14)
+end
+
+fig = Figure(size=(800, 400), fontsize=20)
+
+ax = Axis(fig[1, 1];
+    xreversed = true,
+    xticks = LinearTicks(5),
+    yticks = LinearTicks(5),
+    aspect = DataAspect(),
+    xlabel = "X (Earth-Moon Distance)",
+    ylabel = "Y (Earth-Moon Distance)",
+    title = "Unstable and Stable Invariant Manifolds",
+    titlesize = 24,
+)
+
+idxs = (:x, :y)
+
+# TODO: replace this manual workaround when
+# https://github.com/SciML/SciMLBase.jl/issues/697#issuecomment-2135801331
+# is addressed
+for (traj, color) in zip(unstable, resample_cmap(:blues, length(unstable)))
+    plot!(ax, traj; idxs, color)
+end
+
+for (traj, color) in zip(stable, resample_cmap(:blues, length(stable)))
+    plot!(ax, traj; idxs, color)
+end
+
+scatter!(ax, [1-μ], [0]; marker='⨯', color=:black, markersize=50, label="Moon")
+
+fig
+```