AstrodynamicalSolvers.jl
Common solvers within orbital mechanics and astrodynamics.
Installation
pkg> add AstrodynamicalSolvers
Getting Started
This package currently provides periodic orbit, and manifold computations within Circular Restricted Three Body Problem dynamics.
Periodic Orbits
This package contains differential correctors, and helpful wrapper functions, for finding periodic orbits within Circular Restricted Three Body Problem dynamics.
using AstrodynamicalSolvers
using AstrodynamicalModels
using OrdinaryDiffEq
using Plots
μ = 0.012150584395829193
planar = let
u, T = halo(μ, 1) # lyapunov (planar) orbit
u = [u.x, 0, 0, 0, u.ẏ, 0]
problem = ODEProblem(CR3BFunction(), u, (0, T), (μ,))
solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14)
plot(solution, idxs=(:x,:y,:z), title = "Lyapunov Orbit", label=:none, size=(1600,900), dpi=400, aspect_ratio=1)
end
extraplanar = let
u, T = halo(μ, 2; amplitude=0.01) # halo (non-planar) orbit
u = [u.x, 0, u.z, 0, u.ẏ, 0]
problem = ODEProblem(CR3BFunction(), u, (0, T), (μ,))
solution = solve(problem, Vern9(), reltol=1e-14, abstol=1e-14)
plot(solution, idxs=(:x,:y,:z), title = "Halo Orbit", label=:none, size=(1600,900), dpi=400, aspect_ratio=1)
end
plot(planar, extraplanar, layout=(1,2))
Manifold Computations
Manifold computations, provided by AstrodynamicalCalculations.jl
, can perturb halo orbits onto their unstable or stable manifolds.
using AstrodynamicalSolvers
using AstrodynamicalCalculations
using AstrodynamicalModels
using OrdinaryDiffEq
using LinearAlgebra
using Plots
μ = 0.012150584395829193
unstable = let
u, T = halo(μ, 1; amplitude=0.005)
u = [u.x, 0, u.z, 0, u.ẏ, 0]
Φ = monodromy(u, μ, T)
ics = let
problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, T), (μ,))
solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(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, 2T), (μ,)),
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
u, T = halo(μ, 2; amplitude=0.005)
u = [u.x, 0, u.z, 0, u.ẏ, 0]
Φ = monodromy(u, μ, T)
ics = let
problem = ODEProblem(CR3BFunction(stm=true), vcat(u, vec(I(6))), (0, T), (μ,))
solution = solve(problem, Vern9(), reltol=1e-12, abstol=1e-12, saveat=(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.1T), (μ,)),
prob_func=(prob, i, repeat) -> remake(prob; u0=perturbations[i]),
)
solution = solve(problem, Vern9(), trajectories=length(perturbations), reltol=1e-14, abstol=1e-14)
end
figure = plot(;
aspect_ratio = 1.0,
background = :transparent,
grid = true,
title = "Unstable and Stable Invariant Manifolds",
size = (1600,900),
dpi = 400,
)
plot!(figure, unstable, idxs=(:x, :y), aspect_ratio=1, label=:none, palette=:blues)
plot!(figure, stable, idxs=(:x, :y), aspect_ratio=1, label=:none, palette=:blues)
scatter!(figure, [1-μ], [0], label="Moon", xlabel="X (Earth-Moon Distance)", ylabel="Y (Earth-Moon Distance)", marker=:x, color=:black, markersize=10,)