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))
Example block output

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,)
Example block output