Intersection of trajectories

intersections.jl

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