Last active
September 27, 2021 03:23
-
-
Save mforets/a87d294b41e850ee8929144ad2b0c1b9 to your computer and use it in GitHub Desktop.
Intersection of trajectories
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using Plots | |
using Unitful | |
using DifferentialEquations | |
using GeneralAstrodynamics | |
using LazySets | |
# Find a periodic orbit near Earth (orbit, period) | |
Oₑ, Pₑ = halo(SunEarth; Az=100_000u"km", L=2) | |
# Find a periodic orbit near Jupiter (orbit, period) | |
Oⱼ, Pⱼ = halo(SunJupiter; Az=300_000u"km", L=1) | |
# Compute an unstable manifold departing Earth | |
Mₑ = manifold(Oₑ, Pₑ; trajectories=50, saveat=0.01, duration=3Pₑ, eps=1e-6, Trajectory=Val{false}, direction=Val{:unstable}) | |
# Compute a stable manifold arriving at Jupiter | |
Mⱼ = manifold(Oⱼ, Pⱼ; trajectories=50, saveat=0.01, duration=3Pⱼ, eps=-1e-6, Trajectory=Val{false}, direction=Val{:stable}) | |
# Assume there's some function that takes a state vector, | |
# and transforms it to a common reference frame | |
# magic_transformation(state::AbstractVector) = ... # another state vector | |
# What's the intersection of Mₑ and Mⱼ? | |
plot(Mₑ, vars=:xy) | |
plot!(Mⱼ; vars=:xy, palette=:greens) | |
xcoords_manifold(M) = reduce(vcat, [[M[i][j].x for j in 1:length(M[i])] for i in 1:length(M)]) | |
ycoords_manifold(M) = reduce(vcat, [[M[i][j].y for j in 1:length(M[i])] for i in 1:length(M)]) | |
points_manifold(M) = [Singleton([x, y]) for (x, y) in zip(xcoords_manifold(M), ycoords_manifold(M))] | |
function group_manifold_points(M, ndivx=20, ndivy=20, dom=nothing) | |
points = points_manifold(M) | |
if isnothing(dom) | |
dom = box_approximation(UnionSetArray(points)) | |
end | |
part = split(dom, [ndivx, ndivy]) | |
idx_dom = [] | |
for D in part | |
idx = findall(p -> p ⊆ D, points) | |
isempty(idx) && continue | |
push!(idx_dom, idx) | |
end | |
UnionSetArray([ConvexHullArray(points[ii]) for ii in idx_dom]) | |
end | |
Xₑ = group_manifold_points(Mₑ, 40, 40); | |
Xⱼ = group_manifold_points(Mⱼ, 40, 40); | |
plot(Xₑ, c=:lightblue, alpha=.5) | |
plot!(Mₑ, vars=:xy, c=:blue) | |
plot(Xⱼ, c=:lightgreen, alpha=.5) | |
plot!(Mⱼ, vars=:xy, c=:green) | |
# bounding box | |
O = Hyperrectangle(low=[0.9, -0.1], high=[1.1, 0.1]) | |
Xₑ = group_manifold_points(Mₑ, 40, 40, O); | |
Xⱼ = group_manifold_points(Mⱼ, 20, 20, O); | |
idxⱼ = findall(D -> !isdisjoint(O, D), Xⱼ.array) | |
idxₑ = findall(D -> !isdisjoint(O, D), Xₑ.array); | |
# finds intersection (naive way) | |
out = [] | |
for i in idxₑ | |
for j in idxⱼ | |
a = overapproximate(Xₑ.array[i], HPolygon, 1e-3) | |
b = overapproximate(Xⱼ.array[j], HPolygon, 1e-3) | |
R = intersection(a, b) | |
if !isempty(R) | |
push!(out, R) | |
end | |
end | |
end | |
Xint = identity.(out); | |
plot(Xⱼ, c=:lightgreen) | |
plot!(Mⱼ, vars=:xy, c=:green) | |
plot!(Mₑ, vars=:xy, c=:blue) | |
plot!(Xₑ, c=:lightblue) | |
plot!(Xint, lc=:magenta, c=:white, lw=2.0) | |
xlims!(0.9, 1.1); ylims!(-0.1, 0.1) | |
@show sum(area, Xint) | |
sum(area, Xint) = 0.0005889003496004155 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment