Extraction and Visualization of Poincaré Map Topology for Spacecraft Trajectory Design
EExtraction and Visualization of Poincar ´e Map Topology forSpacecraft Trajectory Design
Xavier Tricoche, Wayne Schlei, and Kathleen C. Howell
Earth Moon1 2 3 4 5 6
Fig. 1. An overview of the 1450 periodic orbits (most of which previously unknown) extracted by the proposed method in the Poincar´emap of the Earth-Moon three-body system at an energy level that enables spacecraft motion across the entire phase space. A sampleof orbits shows the corresponding trajectories around the Earth and the Moon. Our method allows spacecraft trajectory designers tovisually find and evaluate natural connections between these orbits as they seek energy-efficient paths that satisfy mission constraints.
Abstract —Mission designers must study many dynamical models to plan a low-cost spacecraft trajectory that satisfies missionconstraints. They routinely use Poincar´e maps to search for a suitable path through the interconnected web of periodic orbits andinvariant manifolds found in multi-body gravitational systems. This paper is concerned with the extraction and interactive visualexploration of this structural landscape to assist spacecraft trajectory planning. We propose algorithmic solutions that address thespecific challenges posed by the characterization of the topology in astrodynamics problems and allow for an effective visual analysisof the resulting information. This visualization framework is applied to the circular restricted three-body problem (CR3BP), where itreveals novel periodic orbits with their relevant invariant manifolds in a suitable format for interactive transfer selection. Representativedesign problems illustrate how spacecraft path planners can leverage our topology visualization to fully exploit the natural dynamicspathways for energy-efficient trajectory designs.
Index Terms —Trajectory planning and design, Poincar´e map, dynamical systems, topology extraction, invariant manifolds, chaos,visual analysis
NTRODUCTION
The design of a spacecraft trajectory is key to the success of any spacemission. The spacecraft path must deliver the scientific objectives ofthe mission under the constraints imposed by the laws of physics anda targeted mission price. Although many factors influence the mission • Xavier Tricoche is an Associate Professor of Computer Science at PurdueUniversity, e-mail: [email protected]. • Wayne Schlei is a Mission Design Engineer, JHU Applied Physics Lab,e-mail: [email protected]. • Kathleen Howell is the Hsu-Lo Professor of Aeronautics and Astronauticsat Purdue University, email: [email protected] received xx xxx. 201x; accepted xx xxx. 201x. Date of Publicationxx xxx. 201x; date of current version xx xxx. 201x. For information onobtaining reprints of this article, please send e-mail to: [email protected] Object Identifier: xx.xxxx/TVCG.201x.xxxxxxx cost, the price is strongly driven by the spacecraft mass. Beside thespacecraft operational equipment ( e.g. , antenna, engines, solar arrays,etc...), this mass is made up of payload and propellant . The payload isthe collection of scientific instruments required to acquire the missiondata. Course corrections or maneuvers are accomplished by perform-ing a change in spacecraft velocity (called ∆ V ) while expelling propel-lant. Although propellant is necessary to perform maneuvers, payloadmass is typically favored over propellant mass. Such a trade-off pro-duces more scientific return for the mission and may reduce its cost.The role of a spacecraft trajectory designer is therefore to devise apathway that minimizes the amount of propellant required to transportthe vehicle to mission objectives.An essential aspect of the path planner’s work consists in findingsuitable trajectory candidates within relevant multi-body gravitationalmodels. Specifically, designers strive to exploit the dynamics that nat-urally exists in the model to control the spacecraft’s path while requir-ing little or no maneuvers. Widely used in practice, the circular re- a r X i v : . [ n li n . C D ] S e p tricted three-body problem (CR3BP) offers a versatile model in whichthe gravitational field is defined by two massive bodies (the Sun anda planet or a planet and a moon). The CR3BP exhibits a rich dynam-ics, and lends itself to a realistic simulation of the spacecraft behaviorduring a planned mission. Since the corresponding design space, the phase space , is high-dimensional, the dynamics is preferrably stud-ied on a properly chosen planar surface utilizing the Poincar´e map (orfirst-return mapping). The Poincar´e map offers a phase-space snap-shot of all crossing trajectories while reducing the dimension of theproblem. Fig. 1 shows a Poincar´e map of the Earth-Moon CR3BP inwhich periodic trajectories have been highlighted. Yet, even in thismore tractable setup, the analysis of the CR3BP remains a difficulttask. Knowing the topological skeleton of the Poincar´e map, whichindicates how orbital structures are naturally interconnected, wouldsupply spacecraft path-planning with a plethora of options and path-ways that save on propellant usage. Unfortunately, the automatic ex-traction of the CR3BP’s Poincar´e map topology is challenging due tonumerical sensitivities, the extensive presence of chaotic dynamics,and violated assumptions in the construction of the map. As a result,the topological insight that has been used so far to design low-costtransfers remains very limited. Instead, spacecraft trajectory designersusually probe Poincar´e maps in basic puncture plots, which they thenuse to initiate the search for relevant connections through numericaloptimization, a process that can be ineffective and tedious.In this context, this paper makes the following contributions. First,we present a number of significant technical solutions to the extractionof the Poincar´e map topology that decisively improve upon existingmethods in the planar CR3BP. Second, we propose an interactive vi-sualization that delivers this topological information in an actionablevisual form to the space trajectory designer. Specifically, we record thehierarchical make-up of the separatrices of the topology, which allowsdesigners to instantaneously navigate these manifolds up and downstream at no additional computational cost as they explore the visual-ization. We demonstrate the effectiveness of our solution in the contextof representative design scenarios, in which our topology visualizationenables the rapid definition of free-flowing connections between dis-tant spatial locations.
ELATED W ORK
The three-body problem defines a vector field and the topologicalstructure of the CR3BP that we aim to visualize in a Poincar´e mapis the discrete signature of a three-dimensional vector field topol-ogy. There exists an abundant literature in the visualization commu-nity focused on the use of topological methods to represent vectorfields [12, 20, 10]. In the discrete realm, L¨offelmann et al. intro-duced multiple methods for the visualization of maps [17, 15, 14, 16],thereby emphasizing the continuous structures that control the dynam-ics. This work followed the seminal contributions of Abraham andShaw to the expressive visualization of dynamical systems [1]. Peik-ert and Sadlo applied a Poincar´e map approach to the visualization ofvortex rings in a flow recirculation bubble [18, 19]. Their work in-troduced a method to find the location of center-type fixed points andrevealed the tangles formed by the separatrices of the topology. Tric-oche et al. proposed to compute the finite-time Lyapunov exponent(FTLE) [9] of the Poincar´e map to obtain an image-based visualiza-tion of separatrices in area-preserving maps [30]. Sanderson et al. pre-sented a method to extract the geometry of the island chains present inPoincar´e maps [23] from the successive returns of randomly selectedstreamlines. They also described a method to infer the location of fixedpoints from that geometry. Tricoche et al. later proposed a method toextract the fixed points of area-preserving maps along with their as-sociated separatrices [29] and applied it to magnetic confinement ina fusion reactor. Sagrista et al. introduced an extension of FTLE toinertial dynamics and applied it to the analysis of this dynamic in high-dimensional phase space [22].The visualization of the Poincar´e map has also been considered inthe context of astrodynamics and space mission design applications.Short et al. applied the concept of Lagrangian coherent structures [9]to the visualization of the Poincar´e map in the CR3BP and other multi- body systems [26, 27]. Haapala et al. presented a glyph-based solutionfor the visualization of four-dimensional Poincar´e maps in the study ofthe spatial three-body problem [8]. Davis et al. proposed a visualiza-tion method for multidimensional Poincar´e maps to identify designoptions around Saturn’s poles [4]. Most closely related to the presentwork, Schlei et al. [25] presented several modifications to the fixedpoint extraction method of Tricoche et al. [29] to improve its results inmulti-body gravitational environments. We discuss their work in moredetail in the presentation of our algorithm in Sect. 4.Unfortunately, none of the techniques presented to date are ableto capture the complex topology of the CR3BP, nor do they supportthe kind of visual analysis of the topology that spacecraft trajectorydesigners need.
PPLICATION B ACKGROUND
We provide in this section a brief introduction to space trajectory de-sign, which provides the application context and the motivation for thework presented in the remainder of this paper. We start with a shortpresentation of the circular restricted three-body problem in its planarform before discussing its analysis through Poincar´e maps.
The gravitational model considered in this paper describes the mo-tion of a spacecraft under the influence of two celestial bodies thatform an orbital system ( e.g. , Earth and Moon or Sun and Earth). Themotion of a spacecraft under the influence of the combined grav-itational field is then modeled as the circular restricted three-bodyproblem (CR3BP). Assume a pair of gravitating bodies, P and P with corresponding masses m > m and associated gravity param-eter µ = m / ( m + m ) , revolve around their common barycenter B incircular orbits (see Fig. 2, left). We further assume that the spacecraftresides in the orbital plane (that is z ≡ state vector x = [ x , y , ˙ x , ˙ y ] T (with ˙ a denoting the time derivative dadt ). Coordinates are expressed ina rotating reference frame with the origin at the barycenter, the axis ˆ x aligned with the −−→ P P line and ˆ y is aligned with the velocity vector of P with respect to P . B
This section reviews basic theoretical results about area-preservingmaps and Poincar´e maps that apply to the PCR3BP. Additional detailscan be found in classical references [7, 13].A dynamical system associated with a vector field v defines a flowmap x f with ˙ x f = v ( x f ) such that x f ( t , t , x ) describes the transportfrom an initial state x at time t to its later state at time t .Let Σ represent a plane that is transverse to the flow and let x be an initial state on Σ . The Poincar´e map, or the first-return map P ( x ) : = x (cid:55)→ P Σ ( x ) corresponds to the first crossing of Σ by thetrajectory starting at x . Multiple iterations of the Poincar´e map arethen computed by compounding the first return map, e.g. , P p ( x ) = P Σ ( P Σ ( ... P Σ ( x ))) for p returns. Both the initial state and first returnto Σ are shown on the yellow plane in Fig. 3(a). Later sections refer tothe vector-valued (displacement) mapping ∆∆∆ p associated with a givenperiod p , which is defined as: ∆∆∆ p ( x ) = P p ( x ) − x . (4)Three dynamic behaviors co-exist on a Poincar´e map of the pla-nar PCR3BP: periodicity , quasi-periodicity , and chaos . Visible inFig. 3(a), a periodic state, x ∗ , returns to the same state through thePoincar´e map, i.e., ∆∆∆ p ( x ∗ ) = p represents the number of iter-ations required for a p -periodic trajectory to complete an orbit. These p distinct returns are called fixed points of the Poincar´e map and, ow-ing to the aera-preserving nature of the PCR3BP, they can either be centers or saddles . Refer to Fig. 3(b). So-called stable and unsta-ble manifolds (known as separatrices in vector field topology) emergefrom the saddle fixed points and flow into and out of the periodic or-bits, respectively. A fundamental feature of Poincar´e map topology is The plane is transverse if the vector field is nowhere tangent to Σ the existence of connections between saddle fixed points, in which sta-ble and unstable manifolds intersect an infinite number of times, cre-ating tangles as seen in Fig. 3(b). Connections between saddle fixedpoints of the same periodic orbit are called homoclinic connections ,whereas those that join saddle fixed points of two distinct orbits arecalled heteroclinic connections . y
It can be shown that, for a given value of the Jacobi constant C , thedynamics of the PCR3BP is confined to a torus [28], in which the mo-tion is then characterized by the so-called winding number w = ω ω ,where ω and ω are the poloidal (measured in small circle) and toroidal (measured in large circle) rotation frequencies, respectively.The winding number permits us to classify trajectories: rational num-bers w = qp , p , q ∈ N ∗ , p and q mutually prime correspond to periodicorbits. In this case, q corresponds to the number of poloidal rotationsperformed during p toroidal rotations and p is the period of the fixedpoint. In contrast, quasi-periodic trajectories possess irrational wind-ing numbers: a quasi-periodic orbit will never trace exactly the samepath along the torus, and in the case of chaotic trajectories, the windingnumber is undefined.In the vicinity of a fixed point of period p , the spatial derivative ofthe flow map x f after p returns is known as the monodromy matrix M . The eigenvalues of M determine the type of a fixed point. Sincethe flow map is area-preserving, these eigenvalues come in reciprocalpairs. The distinction between stable, center, and unstable types isdone as follows. | λ i | <
1: stable , | λ i | =
1: center , | λ i | >
1: unstable . (5)In the case of a saddle-type fixed point, one eigenvalue has a stabletype while the other eigenvalue has an unstable type. The correspond-ing eigenvectors ( e i with i = S , U for stable and unstable, respectively)indicate the tangent of its invariant manifolds W i .An alternative stability classification is possible through a stabilityindex , ν SI , defined as ν SI = ( tr ( M ) − ) . (6)Orbits are unstable when | ν SI | > . OINCAR ´ E M AP T OPOLOGY E XTRACTION IN THE
PCR3BP
As noted previously, no existing method can successfully extract thetopology of the Poincar´e map in the PCR3BP. Nonetheless, our workis closely related to the work of Schlei et al. [25]. In this section, wefirst summarize the main points of this existing approach and discussits shortcomings. We then present our algorithmic contributions. Notethat in the following we refer to a maximum period p max to consider.A maximum period is needed because the PCR3BP possesses a fractal topology that contains an infinite number of periodic orbits and asso-ciated invariant manifolds for arbitrarily high values of the period. Amaximum period is therefore a means to bound the scope of the topol-ogy extraction. From a practical standpoint, a maximum period alsoamounts to an upper bound on the time and distance that a spacecraftwould need to travel to return to the same point. .1 Existing Approach and Limitations Schlei et al. proposed a method to extract fixed points and separatricesof the Poincar´e map topology [25], which extends to the PCR3BP amethod proposed by Tricoche et al. [29]. Its main input parametersare p max : maximum period of the Poincar´e map to consider and G :regular grid that covers the considered domain of the Poincar´e sectionin the ( x , ˙ x ) , y = S t a b l e S t a b l e Unstable Δ ¯ p ( x ) P ¯ p ( x ) x Unstable discontinuity
A. B.C.D.
Fig. 4. Poincar´e index computation for ∆∆∆ ¯ p . A. index computed alongthe edges of a cell containing a saddle fixed point of period ¯ p . B. thecorresponding ∆∆∆ ¯ p values complete 1 clockwise rotation (index -1). C.Nonlinearity in the rotation of ∆∆∆ ¯ p is resolved via adaptive sampling. D.Adaptive sampling fails to resolve orientation change when a transver-sality violation creates a discontinuity of ∆∆∆ ¯ p . Map sampling:
Evaluate the Poincar´e map, through numerical inte-gration of the equation of motion (Equation 2), at each vertex of thegrid G , for a number of iterations k ≥ p max . For each computed or-bit, record a vector W = ( w xy , w x ˙ x , w ˙ x , y ) of 3 winding numbers (seeSect. 3.2) in the ( x , ˙ x ) , ( x , y ) , and ( ˙ x , y ) planes, respectively. For eachwinding number w , compute a best rational approximation w ≈ q ∗ / p ∗ such that p ∗ ≤ p max and store p ∗ . In general, the computed orbit willbe quasi-periodic (see Sect. 3.2) and w will be an irrational number. Index computation:
For each cell of the grid G , consider the set ofperiods p ∗ recorded at its 4 vertices. For each such period ¯ p , computevia adaptive sampling the Poincar´e index of the displacement map ∆∆∆ ¯ p (see Equation 4), i.e., the number of counterclockwise rotations of ∆∆∆ ¯ p , along the cell edges. Refer to Fig. 4 for an illustration. Assuminga fine enough sampling , the values of the index are expected to beeither − + Fixed point extraction:
If the index is non-zero, find the location ofthe fixed point, which must satisfy ∆∆∆ ¯ p = , using a Newton method.If a saddle-type fixed point is found, compute the eigenvectors ofthe monodromy matrix M (the spatial derivative of the flow map,Sect. 3.2). Separatrices construction:
For each saddle-type fixed point, com-pute the stable and unstable manifolds ( i.e., separatrices) associatedwith its eigenvectors using an iterative shooting method [5]. Theexisting method [25] suffers from significant limitations that make itunsuitable for mission design applications. Indeed, it can only find avery limited number of periodic orbits and corresponding fixed points,and does so at great computational expense. Specifically, the resultsreported in the paper [25] for the Earth-Moon system correspond toaround 50 periodic orbits detected, a result that was only possiblethanks to the tedious repeated application of the method in 100 smallsubdomains, which was required to use the index to detect fixed pointsin cells. Similarly, the manifold construction method used previouslywas only successful in a small number of cases ( e.g. , those related toLyapunov orbits, see Sect. 3.2 and Fig. 2) because it is oblivious to the The cell must contain no more than one fixed point of the same period forthe detection criterion to be valid. fundamental issue of transversality violation that is discussed below.In contrast, our solution can reliably and more efficiently characterizethe topological structure of the PCR3BP to support an effective visualanalysis. For comparison, our method found 1450 periodic orbits inthe same system by inspecting 2 large domains and 2 densely popu-lated smaller domains around the Moon. Similarly, the manifold ex-traction algorithm presented hereafter succeeded in constructing mostof the invariant manifolds corresponding to periodic orbits with lim-ited instability ( | ν SI | < , see Equation 6). Indeed, spacecrafts canmaintain their position even on unstable orbits with manageable con-trol effort, up to a certain value of the stability index. As previouslystated, these topological objects are extremely valuable for our appli-cation as they can be used as anchor points and pathways through thedynamics of the PCR3BP to chart a trajectory that satisfies missionobjectives. To permit a reliable detection of fixed points, the Poincar´e index of thedisplacement map ∆∆∆ p (Equation 4) must be evaluated around areas ofthe map that are small enough to contain at most a single fixed pointfor the considered period [28]. Hence, a very high resolution sam-pling yielding tiny cells is typically desirable. However, this approach,which was used in prior work [25], is computationally prohibitive anda more subtle data-driven sampling is needed.Our key observation is that the winding numbers associated witheach trajectory (Sect. 3.2) are locally smoothly varying characteristicparameters within regions of regular dynamics. Therefore, the valuesof all the winding number vectors W = ( w xy , w x ˙ x , w ˙ xy ) (Sect. 4.1) mea-sured in and around a cell should lie within a small interval. Hence,we apply this criterion in each cell to the W values measured at thevertices of the cell, as well as to the trajectories that intersected the interior of the cell during the Poincar´e map sampling. If this criterionis not satisfied, the cell is marked for subdivision.Practically, we use a quadtree-like data structure that records wind-ing numbers both at the sampling vertices and inside the cells. Sinceall Poincar´e section crossings of a given trajectory share the samewinding numbers set W as the initial vertex, each one of the p in-termediate returns is assigned the same W values. These values areadded to the cell containing the return and are then tested as part ofthe subdivision criterion discussed previously. When required, cellsare subdivided with internal data assigned to the corresponding quad-rant within the original cell. A user-specified maximum depth levelparameter, d max , is employed to represent the total number of subdi-vision layers allowed. Cells are also subdivided if any corner has anundefined winding number, which occurs in chaotic regions.Fig. 5 shows the adaptive resolution mesh produced in the domain D EM for a maximum refinement depth d =
3. The initial grid is shownin thick gray lines. The highest resolution is achieved in chaotic re-gions where the dynamics is most complex. This result is excel-lent from an astrodynamics perspective since the saddles embeddedin chaos offer the most versatile transfer opportunities. Regions ofregular dynamics, in contrast, are coarsely resolved, as expected.
An issue that can prevent the computation of the Poincar´e index inprior work [25] is the presence of discontinuities of the ∆∆∆ mappingalong cell edges. Two properties of the PCR3BP can explain this be-havior: highly sensitive dynamics and transversality violation of theflow map for the chosen section Σ . Indeed, the the ability of thePoincar´e map to capture the dynamics of the system is premised onthe assumption that the flow is transverse to the chosen section in theconsidered region. Transversality violations are typically the result ofone of two specific trajectory events. First, trajectories that are tangentto the section along their path generate discontinuities in ∆∆∆ , see Fig. 6.A second event is an intersection by the trajectory of a primary bodyin CR3BP (see Fig. 6). ig. 5. Adaptive cell subdivision based on the winding number set W applied to the domain D EM with parameters C = . and d max = .Fig. 6. Transversality violation types in the PCR3BP. Left: section tan-gency, Right: singularity intersection. We perform the evaluation of the Poincar´e index in cells wheretransversality violations occur by considering the behavior of ∆∆∆ (Equa-tion 4) in the limit approaching a transversality violation. A disconti-nuity of ∆∆∆ at some location g on the closed curve Γ ( i.e., the boundaryof a cell) for period p creates a discontinuity in the angle coordinate of ∆∆∆ at g (Fig. 4, D.). However, since the limits of the value ∆∆∆ ( g ) exist onboth sides of g , the Poincar´e index can be expressed as the summationof improper integrals κ = π (cid:73) Γ d α ( ∆∆∆ ) = π (cid:18) (cid:90) g γ d α ( ∆∆∆ ) + (cid:90) γ g d α ( ∆∆∆ ) (cid:19) , (7)where γ is a starting point along Γ ( γ (cid:54) = g ) and α ( ∆∆∆ ) is the angleformed by ∆∆∆ with the x axis.The edge sampling approach used in prior work [25] is thereforemodified to detect transversality violations and ensure fine samplingof the edge around the corresponding discontinuity. The success of a numerical search depends heavily on the quality ofthe initial guess. Our solution starts by sampling the displacement map ∆∆∆ at a set of regularly distributed positions within the cell with nonzeroindex. Practically, if the considered variable is ζ = x − s with s repre-senting the (unknown) saddle-type fixed point location and x = ( x , y ) is the dynamic state of the system (see Equation 2), then a quadratic model of the sampled dynamics is formed as ˙ ζ = A s ζ + ζ T Q ζ . (8)Note, A s is a 2 × Q is a 2 × × Q = ζ = , which yields an initial guess for aniterative numerical search. The Newton search and multiple shooting solutions used in priorwork to find fixed points often fail. To remedy this situation, we adopta multipronged approach. Starting with a single shooting method forfixed point refinement, we switch to a multiple shooting method ifthis first attempt fails. If the second solution fails as well, we applya refined precision Newton method, which is significantly more com-putationally expensive than the previous two but possesses strongerconvergence properties. If a saddle type is identified, the construction of the stable and un-stable manifolds constitutes the last step of our topology extractionand we derive eigenvectors and stability index (Equation 6) from M .Prior work [29, 25] constructs invariant manifolds through a series oftwo-point boundary problems that aim to ensure smoothness and finesampling of the manifold [5]. Unfortunately, this solution does nothandle the issue raised by transversality violations in the PCR3BP (seeSect. 4.3), nor does it provide any guidance to accommodate the nu-merical challenges associated with this particular system. We describein the following our improvements of this method. We then presenta data structure that allows the mission designer to instantly navigatethe manifolds up- and downstream as demonstrated in Sect. 5. Manifold Extraction with Curve-Refinement.
Consider two adja-cent positions φ i and φ i + that form a segment w = φ i φ i + on the man-ifold. We further assume that both positions are close enough suchthat linear interpolation between these two positions yields positionsthat are themselves on the manifold, which is a valid approximationas long as the manifold does not come too close to another saddlefixed point. Applying the Poincar´e map P ¯ p ( ¯ p is the period of the cor-responding saddle fixed point) to any such intermediate position willresult in a new position further downstream on the manifold. Referto Fig. 7 (top). The algorithm used in prior work [25] performs an Fig. 7. Top: Schematic of a 1D invariant manifold curve on the Poincar´esection. Bottom: Generating new downstream manifold points and seg-ments through a transversality violation. adaptive sampling of the segment w by the Poincar´e map, controlledby curve quality checks, to construct the next segment = φ i + φ i + onthe manifold and ultimately extract the entire manifold [5]. Our so-lution follows the same approach while simultaneously checking fortransversality violations.The heuristics for detecting transversality violations duringPoincar´e index evaluation are reapplied alongside the curve-refinement criteria. If a downstream transversality violation is detectedbetween consecutive segment samples, the segment is bisected on thatinterval. Subdivision continues until the distance between consecutivepoints reaches a user-prescribed minimal distance ( u min , which cor-responds to a relative distance ∆∆∆ τ min on the segment). An exampleis depicted in Fig. 7 (bottom) where downstream mappings are color-coded by their initial position on the active segment: a downstreamtransversality violation exists between φ i and the midpoint x mw , andsubdivision localizes the separation condition when the parameter dif-ferential is below ∆∆∆ τ min . Multiple shooting methods introduce intermediate steps to the solution. topping Criteria.
We found two criteria primarily effective in con-trolling the useful downstream length of an invariant manifold as itstarts to form tangles in the vicinity of another saddle point. Referto Fig. 3(b). The first stopping criteria tracks a practical measure forspacecraft trajectory planning. The manifold construction algorithmcreates a parent-child link between upstream and downstream seg-ments. Our algorithm caps the manifold progression by stopping whenthe depth of the resulting tree reaches a maximum depth d w , max = P p ( x ) ). The second stopping condition observes simultaneous ad-vection of manifolds from the same periodic orbit for the detection ofso-called homoclinic connections . Screening Computations
Given the high computational cost of man-ifold construction, we perform several tests beforehand to prevent un-necessary computations. First, we pre-screen for potentially impracti-cal structures for spaceflight, i.e., periodic orbits with exceptionallyhigh instabilities. A threshold cutoff is established on stability in-dex magnitude at | ν SI | > (refer to Equation 6), whereby a lower | ν SI | cutoff can be used to further reduce the overall workload of ourmethod.A second important observation is that a lower bound is necessaryfor u min (which regulates upstream manifold segment subdivisions) forrealistic spaceflight. The accuracy limitations of sensors and enginesprevent an exact knowledge of the spacecraft position and speed andmake maneuvers impractical below a certain ∆ V magnitude. Space-craft state determination outside of low Earth orbits is limited to anaccuracy of 3 km for position and 0.1 mm/s for velocity based onmeasurement error of standard capabilities [31]. Practically, we re-quire u min values above 2 × − (nondimensional map displacement)for the Earth-Moon system, which is equivalent to 2.05 cm/s for ve-locity and 7.69 km for position. Note that the value of u min is differentacross CR3BP systems. Suggested values are provided in appendix. Manifold Arc Extraction
Implementing manifolds within trajectorydesign applications requires the ability to select an invariant mani-fold state from the Poincar´e map visualization and reconstitute thearc. Simply propagating a user-selected invariant manifold state up-stream ( i.e., forward-time for W S and reverse-time for W U ) often failsbecause of the flow-dividing nature of the invariant manifolds, fromwhich neighboring trajectories diverge exponentially.Upstream manifold arc reconstruction with our method circumventsthe challenges associated with propagating manifold arcs upstreamwith a very effective data structure format resulting from the mani-fold advection procedure. Our manifold extraction process advectsstates sampled from an upstream manifold segment to create a newgroup of downstream segments through the Poincar´e map; this pro-cedure naturally forms a manifold segment tree . To illustrate, samplemanifold segments near the origin location of the advection procedureare arranged as a staircase schematic indicating depth levels ( d w ) as inFig. 8; each step down symbolizes the downstream progression to thenext group of segments after p map iterates. To reconstruct the entire Fig. 8. As the construction of the invariant manifold progresses viacurve-refinement, the spawning of new manifold segments generatesa tree structure that can be employed for accessing data. upstream trajectory from the fixed point to an arbitrary user-selectedpoint, the downstream mappings recorded during manifold construc-tion supply a reliable framework to return the intermediate trajectorystates between the parent and child manifold segments. Practically,the local (linear) coordinates of a state on a child segment are mappedupstream to a state with the same local coordinates on the parent seg-ment, which eventually leads to the immediate vicinity of the saddle fixed point from which the manifold originated.
OPOLOGY V ISUALIZATION
We present in this section the visualizations of the PCR3BP that weremade possible by our topology extraction method and illustrate thevariety of structures that can be observed in this context.
We visualize in the following the fixed points that our method wasable to discover. Among these fixed points are many that are eitherchallenging to extract with existing means ( e.g. , they require veryhigh resolution sampling in a specific region and lengthy numericalsearch), or are presented here for the first time. It is important to notethat the topological complexity (including the number of fixed points )of the PCR3BP varies significantly with the considered energy level,as quantified by the Jacobi constant C (Equation 3). Specifically, highvalues of C correspond to low energy and exhibit a simpler topologythan lower C values, i.e., higher energy levels. Both cases are consid-ered next to illustrate this effect. Fig. 9. Fixed points and selected periodic orbits found in the Earth-Moonsystem at C = . . Practically, in the EM system, C = . inwhich it is impossible to reach the Moon from the Earth as they be-long to disconnected regions of the dynamics. In this case, our visu-alization of the identified fixed points, shown in Fig. 9, reveals sad-dles and centers grouped in island chains (repeated patterns similarto what is shown in Fig. 3(b)), as well as isolated saddle-type fixedpoints embedded in chaotic regions. In this and following images,the fixed points belonging to the same periodic orbit are assigned thesame color. Nonetheless, the large number of periodic orbits makesit challenging to visually identify all the fixed points associated withany given orbit. It is instead more insightful to analyze the geometry ofspecific periodic orbits while highlighting their associated fixed points.Their visualization, displayed in the xy plane in Fig. 9, shows each se-lected periodic orbit as a closed curve, sometime with confoundinggeometric complexity. The background image is formed by orbitalconvolution, an extension of line integral convolution method [2] tomaps [29, 25]. At this energy level, transversality violations are rarewithin the considered region and the fixed point extraction is relativelystraightforward with little need for cell subdivision.A more complex case is considered next, in which a higher energylevel is selected where chaos exists throughout the domain of interest.At C = .
96 the Earth-Moon system permits trajectories everywherein the xy plane, and in particular between the Earth and the Moon.Broader sampling parameters are applied over a larger analysis do-main (see Trial 1 in Table 1). We find fixed points throughout D EM but only a small number in close proximity of the Moon. Increasingthe sampling resolution near the Moon allows us to uncover a large We only consider fixed points of period p < p max This energy level is still much higher than typical levels in low Earth orbit,for which C ∈ [ , ] umber of additional fixed points in that region. Refer to Trials 2, 3,and 4 in Table 1. As shown in Fig. 1, many fixed points were extractedfor C = .
96 in the EM system, reaching a total of 1450 distinct peri-odic orbits. A closeup view centered on L reveals the density of fixedpoints in that region. Trial C Domain ( x , ˙ x ) (nondim) Resolution l min p max [ . , . ] × [ − . , . ] ×
16 8 × −
121 2.96 [ . , . ] × [ − . , . ] ×
16 8 × −
122 2.96 [ . , . ] × [ − . , . ] × × −
123 2.96 [ . , . ] × [ − . , . ] × × −
64 2.96 [ . , . ] × [ − . , . ] } × × − Table 1. Parameters used in the Earth-Moon system.
Among the fixed points found at C = .
96, many novel saddle-typeperiodic orbits were identified. As shown in Fig. 1, several periodicorbits are commonly known such as Orbit 4 (the L Lyapunov), Orbit3 (stable 3:2 resonant orbit, which means that it completes 3 revolu-tions around the Earth while the Moon completes two), and Orbit 5(the p = e.g. , Orbit1). Yet others like Orbit 6 visit all the aforementioned regions, makingsuch orbits compelling for transfer design. Though the analysis is onlyperformed within the primary analysis domain D EM on the Σ : y = L , L , and L vicinities.Clearly, our results offer a vivid dynamical understanding of this par-ticular system. With fixed points extracted, the Poincar´e map topology structure isfurther characterized by means of our invariant manifold extractionalgorithm, which is first demonstrated in the Earth-Moon system at C = .
2. The large-scale topology extraction result appears in Fig. 10with unstable manifolds ( W U ) and stable manifolds ( W S ) colored inred and blue, respectively. At C = .
2, invariant manifolds are ex-tracted throughout the chaotic areas, thoroughly filling in the phasespace areas between quasi-periodic islands. Our algorithm capturessaddle-center island chains except on some islands near the Moon.Difficulties near the Moon can be explained by numerical sensitivityand numerical error build-up during integration as trajectories pass ex-ceptionally close to the singularity multiple times before completingthe p -th iterate. As with fixed point extraction, constructing invariantmanifolds for the Earth-Moon system at C = . C = .
96 in the EM system, our manifold construction pro-duces a visualization of both stable (blue) and unstable (red) manifoldsfor the periodic orbits shown in Fig. 11. We apply stability filtering( | ν SI | ≤ PPLICATION TO T RAJECTORY D ESIGN
We present in the following the insights gained from the application ofour topology visualization technique to two representative trajectorydesign problems. As shown below, the topological scaffold formedby periodic orbits and invariant manifolds, and further enhanced by (a) Primary analysis domain D EM (b) Zoom-in on indicated domain Fig. 10. The Poincar´e map topology skeleton (unstable manifolds in redand stable in blue) computed with the manifold extraction algorithm inthe Earth-Moon system within the domain D EM at C = . . interactive design support, broadens the design possibilities availableto domain experts, delivering new options and the ability to quicklyexamine trade-off decisions. The manifold arc selection solution presented in Sect. 4.6 was createdto allow for the interactive visual inspection of transfer solutions alongmanifolds. Specifically, this selection mechanism allows the user, forany arbitrary location selected interactively along a manifold, to in-stantly determine the precise path (including entry point) that a space-craft would need to follow to arrive at that location. In particular, givena position where two manifolds intersect in the visualization, the usercan readily query the path, upward and downward, that will be trav-eled from that point with no or minimal correction need. This makesit possible to use the visualization of the topology to quickly examinetransfer possibilities.To illustrate this process, we first consider two simple transfers,namely from Earth proximity to the L Lyapunov orbit (which sur-rounds L , see Fig. 2) on one hand, and to the distant retrograde orbit(DRO) of period p =
3, on the other hand. Note that both targetedorbits are unstable, i.e., their corresponding fixed points are saddles.We start by visualizing their invariant manifolds (see Fig. 12(a)). Notethat the fragmentation of the manifold is due to transversality viola-tions. The stable and unstable manifold pairs are colored with blackand crimson for the L Lyapunov and with blue and red for the p = stable manifold of the targeted orbit since thedynamics along that manifold is known to converge toward the orbit,by definition (see Fig. 3(b)).Using this filtered visualization of the topology, the user can eas-ily select entry points on a L Lyapunov stable manifold (black) and ig. 11. The Poincar´e map topology skeleton (unstable manifolds in redand stable manifolds in blue) computed in the Earth-Moon system at C = . . on a p = in Fig. 12(c)) initially perform an elliptical orbit around theEarth before progressively shifting towards an asymptotic approachof the targeted orbits. An important practical question from a designstandpoint is the duration (or time of flight ) of each transfer. The blackarc enters the L Lyapunov orbit after 38.10 days whereas the blue arcarrives in the p = i.e., maneuver-free) capture trajectories.As second example, we now consider an interesting transfer startingat the p = L and L vicinities while also closely ap-proaching the Moon several times. We refer to this orbit as Orbit O ∗ in the following. The visualization of the invariant manifold curvesof that orbit (displayed in Fig. 13 with indigo and tan colors) revealsto the user that a lot of locations naturally flow into this orbit with thelarge dispersion of stable manifolds. A transfer from the p = A popular topic in the astrodynamics community is the design of alow-cost pathway to an orbit around a moon of Saturn called Ence-ladus. A tantalizing scenario in that regard would be a fully ballis-tic capture, namely a solution in which gravity alone would depositthe spacecraft into a close orbit around Enceladus. Fortunately, ourtopology visualization makes it easy to construct such a path, evenwithout detailed knowledge of the Saturn-Enceladus (
SEnc ) systembeforehand.Practically, we select an energy level C = . p = SEnc system. The corresponding visualization is An inertial frame is a reference frame, e.g. , defined with respect to fixeddistant stars, in which the rotating frame rotates. (a) Manifolds of the L Lyapunov and p = Fig. 12. Invariant manifolds for the L Lyapunov and p = DRO saddle-type orbits extracted in the Earth-Moon system ( C = . ). Selected sta-ble manifold arcs are displayed in the rotating (b) and inertial (c) frames. shown in Fig. 14(a). As in our first example (Sect. 6.1), we are hereagain seeking a transfer to a stable manifold of the DRO to achieveguaranteed convergence. This manifold, shown in light blue, exhibitsmany transfer options in the form of intersections with the dark blueunstable manifolds of the many saddle-type periodic orbits that pop-ulate that region. It can further be seen that these manifold segmentspopulate not only the immediate vicinity of Enceladus but also a re-gion situated on the opposite side of Saturn, shown in Fig. 14(b) overan orbital convolution texture [29] for context.A possible solution consists in selecting the position shown inFig. 14(b) where our topology visualization indicates the presence ofan intersection between the unstable manifold of another periodic or-bit and a segment of the DRO stable manifold on the opposite sideof Saturn [24]. Note that the dark blue unstable manifolds shown inFig. 14(a) densely cover this region as well but are omitted in the sec-ond image for clarity. The trajectory associated with this selection isshown in Fig. 14(c). It can be seen that it performs one and a halfrevolutions around Saturn before a close passage of Enceladus thatprecedes the expected asymptotic approach of the DRO. Note that thisparticular solution is only one out many alternative options that ourtopology visualization automatically reveals, here again affording thedesigner precious flexibility in their selection of a suitable path. ISCUSSION AND C ONCLUSION
We have presented a Poincar´e map topology visualization method thatallows space trajectory designers to interactively discover novel nat-ural transfer options in the circular restricted three-body problem asthey piece together the path of a spacecraft in a planned mission. Tomake this visualization possible, we have devised a technical solutionfor the extraction of the Poincar´e map topology in the CR3BP thatallows us to identify thousands of fixed points and associated mani-folds, the overwhelming majority of which were previously unknown.We have shown how a detailed visualization of the topological skele-ton reveals crucial information to domain experts as they seek low- a) Manifolds for the p = O ∗ (b) H c : p = O ∗ ( ∆ t = .
22 days)
Fig. 13. Invariant manifolds for the p = DRO and Orbit O ∗ extracted byour method and a selected heteroclinic connection in the Earth-Moonsystem ( C = . ). propellant trajectories. In particular, we have presented representativedesign scenarios in which the knowledge of the heteroclinic connec-tions that exist between the invariant manifolds of unstable periodicorbits provided the domain expert with the essential insight neededto construct suitable transfer solutions. Our method, by displaying theconnectivity of orbital structures, offers spacecraft trajectory designersa broad range of interactive options without external computation. Ourcollaborators’ experimentation with this approach suggests that a de-signer could employ automated topological skeletons as a visual inputcatalog from which pathways can be interactively selected, visualized,and assessed.Our approach and its current implementation have some limitationsthat we wish to address. First, the computational cost of our topol-ogy extraction method is significant. While our fixed point extractiontechnique is essentially an embarassingly parallel task that can be com-pleted within a few minutes, the individual manifold construction pro-cedure is an expensive operation that can take many hours to finish ifthe numerical parameters that control its behavior are not chosen prop-erly. We provide in appendix the parameter values that we have foundto strike a good balance between accuracy and computation time. Amore principled way to choose these parameters, or alternatively, away to predict their impact on the compute time would greatly bene-fit this work. A second important issue concerns the accuracy of ourevaluation of the Poincar´e map. By definition, each iteration of thismap requires the integration of a strongly nonlinear differential sys-tem, which is an inherent source of error. It would be important to (a) Poincar´e section in Enceladus area(b) Poincar´e section near L (c) Selected capture arc ( p = W S ) Fig. 14. A plausible region for ballistic capture about Enceladus indi-cated on a Poincar´e section with a selected example p = DRO mani-fold arc in the
SEnc system at C = . . quantify the uncertainty that it introduces in our characterization of thetopology and to convey it visually to the designer. Another outstand-ing problem is the visual complexity of the topology and the difficultyto effectively manage it to allow the spacecraft trajectory designer toquickly access the relevant information. We have so far dealt with thisissue by applying, after the fact, stability thresholds to fixed points andmanifolds but more sophisticated and discriminative filtering criteriawould bring a major improvement. In fact, it would be most desirableto be able to apply this kind of filtering in a pre-processing step in or-der to limit the scope of the computation to the most relevant featuresof the topology for the considered problem. A CKNOWLEDGMENTS
This research was supported in part by NSF CAREER Award R (cid:13) ) for implementation assis-tance with the visualization tools employed in this work. EFERENCES [1] R. Abraham and C. Shaw.
Dynamics – The Geometry of Behavior .Addison-Wesley, 1992.[2] B. Cabral and C. Leedom. Imaging vector fields using line integral con-volution. In
Proceedings of SIGGRAPH 93 , pages 263–270. ACM SIG-GRAPH, 1993.[3] J. M. Danby.
Fundamentals of Celestial Mechanics . Willmann-Bell, Inc.,Richmond, Virginia, 2nd edition, 1992.[4] D. C. Davis, S. M. Phillips, and B. P. McCarthy. Trajectory design forsaturnian ocean worlds orbiters using multidimensional poincar´e maps.
Acta Astronautica , 143:16–28, 2018.[5] J. England, B. Krauskopf, and H. Osinga. Computing one-dimensionalglobal manifolds of poincar´e maps by continuation.
SIAM Journal ofApplied Dynamical Systems , 4(4):1008–1041, 2005.[6] G. G´omez, W. S. Koon, M. Lo, J. E. Marsden, J. Masdemont, and S. D.Ross. Connecting orbits and invariant manifolds in the spatial restrictedthree-body problem.
Nonlinearity , 17(5):1571, 2004.[7] J. Guckenheimer and P. Holmes.
Nonlinear Oscillations, Dynamical Sys-tems, and Bifurcations of Vector Fields . Springer-Verlag, 1983.[8] A. F. Haapala and K. C. Howell. Representations of higher-dimensionalpoincar´e maps with applications to spacecraft trajectory design.
ActaAstronautica , 96:23–41, 2014.[9] G. Haller. Distinguished material surfaces and coherent structures inthree-dimensional flows.
Physica D , 149:248–277, 2001.[10] C. Heine, H. Leitte, M. Hlawitschka, F. Iuricich, L. De Floriani,G. Scheuermann, H. Hagen, and C. Garth. A survey of topology-basedmethods in visualization.
Computer Graphics Forum , 35(3):643–667,2016.[11] K. Howell. Families of orbits in the vicinity of the collinear librationpoints.
J. of the Astronautical Sciences , 49(1):107–125, 2001.[12] R. S. Laramee, H. Hauser, L. Zhao, and F. H. Post. Topology-based flowvisualization, the state of the art. In H. Hauser, H. Hagen, and H. Theisel,editors,
Topology-Based Methods in Visualization 2005 . Springer, 2005.[13] A. J. Lichtenberg and M. A. Lieberman.
Regular and Chaotic Dynamics .Springer-Verlag, New York, New York, 2nd edition, 1992.[14] H. L¨offelmann, H. Doleisch, and E. Gr¨oller. Visualizing dynamical sys-tems near critical points. In , pages 175–184, 1998.[15] H. L¨offelmann and E. Gr¨oller. Enhancing the visualization of character-istic structures in dynamical systems. In
Visualization in Scientific Com-puting ’98 , pages 59–68. Springer, 1998.[16] H. L¨offelmann, T. Kuˇcera, and E. Gr¨oller. Visualizing poincar´e mapstogether with the underlying flow. In H.-C. Hege and K. Polthier, editors,
Mathematical Visualization: Proceedings of the International Workshopon Visualization and Mathematics ’97 , pages 315–328. Springer, 1998.[17] H. L¨offelmann, L. Mroz, E. Gr¨oller, and W. Purgathofer. Stream arrows:enhancing the use of stream surfaces for the visualization of dynamicalsystems.
The Visual Computer , 13(8):359 – 369, 1997.[18] R. Peikert and F. Sadlo. Visualization Methods for Vortex Rings and Vor-tex Breakdown Bubbles. In A. Y. K. Museth, T. M¨oller, editor,
Proceed-ings of the 9th Eurographics/IEEE VGTC Symposium on Visualization(EuroVis’07) , pages 211–218, May 2007.[19] R. Peikert and F. Sadlo. Flow Topology Beyond Skeletons: Visualiza-tion of Features in Recirculating Flow. In H.-C. Hege, K. Polthier, andG. Scheuermann, editors,
Topology-Based Methods in Visualization II ,pages 145–160. Springer, 2008.[20] A. Pobitzer, R. Peikert, R. Fuchs, B. Schindler, A. Kuhn, H. Theisel,K. Matkovi´c, and H. Hauser. The state of the art in topology-based visu-alization of unsteady flow.
Computer Graphics Forum , 30(6):1789–1811,2011.[21] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery.
Nu-merical Recipes: The Art of Scientific Computing . Cambridge UniversityPress, New York, 3rd edition, 2007.[22] A. Sagrista, S. Jordan, A. Just, F. Dias, L. G. Nonato, and F. Sadlo. Topo-logical analysis of inertial dynamics.
IEEE transactions on visualizationand computer graphics , 23(1):950–959, 2016.[23] A. R. Sanderson, G. Chen, X. Tricoche, D. Pugmire, S. Kruger,and J. Breslau. Analysis of recurrent patterns in toroidal magneticfields.
IEEE Transactions on Visualization and Computer Graphics ,16(6):1431–1440, 2010.[24] W. Schlei.
Interactive Spacecraft Trajectory Design Strategies FeaturingPoincar´e map Topology . PhD thesis, Purdue University, 2017. [25] W. Schlei, K. C. Howell, X. Tricoche, and C. Garth. Enhanced visualiza-tion and autonomous extraction of poincar´e map topology.
The Journalof the Astronautical Sciences , 61(2):170–197, 2014.[26] C. Short, K. Howell, and X. Tricoche. Lagrangian coherent structuresin the restricted three-body problem. In
Proceedings of 21st AAS/AIAASpace Flight Mechanics Meeting, New Orleans, Louisiana, Paper No.AAS , pages 11–250, 2011.[27] C. R. Short and K. C. Howell. Lagrangian coherent structures in variousmap representations for application to multi-body gravitational regimes.
Acta Astronautica , 94(2):592–607, 2014.[28] S. H. Strogatz.
Nonlinear Dynamics and Chaos . Westview Press, PerseusBooks Publishing, Cambridge, Massachusetts, 1994.[29] X. Tricoche, C. Garth, and A. Sanderson. Visualization of topologicalstructures in area-preserving maps.
IEEE Transactions on Visualizationand Computer Graphics , 17(12):1765–1774, 2011.[30] X. Tricoche, C. Garth, A. Sanderson, and K. Joy. Visualizing invari-ant manifolds in area-preserving maps. In
Topological Methods in DataAnalysis and Visualization II , pages 109–124. Springer, 2012.[31] J. R. Yim, J. L. Crassidis, and J. L. Junkins. Autonomous orbit naviga-tion of interplanetary spacecraft. In