METHODS AND SYSTEMS FOR TRAJECTORIES, APPROACHES, FLYBYS, LANDINGS, AND ORBITS IN THREE-OR-MORE BODY SYSTEMS

Abstract

Systems and methods are described for computing a trajectory of an object in space to a secondary body (M2) in orbit around a primary body to land on, or capture into orbit, or flyby M2 in a Three-Or-More Body Problem. A special plotting of sampled vectors from M2 are integrated backward using a Poincaré Map to form a “Swiss Cheese plot” to find a nominal trajectory. A funnel-like set of trajectories can be constructed along the nominal trajectory for navigation purposes. A global resonant encounter map over a sphere around M2 can be constructed to provide trajectories to, for example, flyby any point near M2, capture into orbit over any point about M2, land on any point on M2. Besides space exploration, there are many applications to the development of Cislunar space commercialization and colonization including asteroid capture and mining.

Claims

1. A method to provide a nominal trajectory to land an object on a secondary body orbiting a primary body, the method comprising: selecting a Jacobi constant for a nominal trajectory, the nominal trajectory being a landing trajectory for the object to land on the secondary body at a nominal landing site; selecting initial conditions comprising a plurality of velocities tangent to the nominal landing site, each of the plurality of velocities having the Jacobi constant; selecting a surface of section for the primary body; propagating the initial conditions backwards in time for trajectories to intersect the surface of section; producing a Poincaré map from the intersections in the surface of section, keeping track of iterations of the trajectories passing through the surface of section; producing a Swiss Cheese plot of the k-th iterates by producing a Poincaré map of the k-th iterate in at least 2 steps: (i) plot all points using Delaunay variables in 2-dimension; (ii) plot points from the k-th iterate such that the points using Delaunay variables can be distinguished from the points from the k-th iterate; using the Swiss Cheese plot of the k-th iterate by locating a plurality of resonant trajectories of the k-th iterate by locating those points of the k-th iterate on the Poincaré map that are close to the center points between vertical holes in the Poincaré map which determine a resonance; selecting one of the plurality of resonant trajectories as the nominal trajectory.

2. The method of claim 1, wherein the object is a spacecraft.

3. The method of claim 1, wherein the primary body is a planet and the secondary body is a moon of the planet.

4. The method of claim 1, further comprising forming an invariant funnel of trajectories from the initial conditions, the invariant funnel describing a zone of flight for the object surrounding the nominal trajectory.

5. The method of claim 1, further comprising forming a resonant encounter map based on nominal trajectory, the resonant encounter map mapping possible resonances of landing trajectories on a surface of the secondary body.

6. The method of claim 1, wherein the Delaunay variables comprise L and g.

7. The method of claim 1, wherein the primary body, the secondary body, and the object are modeled using Circular Restricted Three-Body Problem.

8. The method of claim 1, wherein the Jacobi constant is determined from a state of a trajectory event at the nominal landing site.

9. A non-transient computer readable medium configured to perform the method of claim 1 when read by a computer.

10. A computer system configured to perform the method of claim 1.

11. A method to provide, for an object, a nominal trajectory to fly-by or orbit near a secondary body orbiting a primary body, the method comprising: selecting a Jacobi constant for a nominal trajectory, the nominal trajectory being a trajectory for the object to fly-by or orbit the secondary body at a nominal target site at a distance R from a center of the secondary body; selecting initial conditions comprising a plurality of velocities tangent to the nominal target site, each of the plurality of velocities having the Jacobi constant; selecting a surface of section for the primary body; propagating the initial conditions backwards in time for trajectories to intersect the surface of section; producing a Poincaré map from the intersections in the surface of section, keeping track of iterations of the trajectories passing through the surface of section; producing a Swiss Cheese plot of the k-th iterates by producing a Poincaré map of the k-th iterate in at least 2 steps: (i) plot all points using Delaunay variables in 2-dimension; (ii) plot points from the k-th iterate such that the points using Delaunay variables can be distinguished from the points from the k-th iterate; using the Swiss Cheese plot of the k-th iterate by locating a plurality of resonant trajectories of the k-th iterate by locating those points of the k-th iterate on the Poincaré map that are close to the center points between vertical holes in the Poincaré map which determine a resonance; selecting one of the plurality of resonant trajectories as the nominal trajectory.

12. The method of claim 11, wherein the object is a spacecraft.

13. The method of claim 11, wherein the primary body is a planet and the secondary body is a moon of the planet.

14. The method of claim 11, further comprising forming an invariant funnel of trajectories from the initial conditions, the invariant funnel describing a zone of flight for the object surrounding the nominal trajectory.

15. The method of claim 11, further comprising forming a resonant encounter map based on nominal trajectory, the resonant encounter map mapping possible resonances of landing trajectories on a surface of a sphere of radius R around the center of the secondary body.

16. The method of claim 11, wherein the Delaunay variables comprise L and g.

17. The method of claim 11, wherein the primary body, the secondary body, and the object are modeled using Circular Restricted Three-Body Problem.

18. The method of claim 11, wherein the Jacobi constant is determined from a state of a trajectory event at the nominal target site.

19. A non-transient computer readable medium configured to perform the method of claim 11 when read by a computer.

20. A computer system configured to perform the method of claim 11.

Description

BRIEF DESCRIPTION OF DRAWINGS

[0029] FIG. 1 shows an example of a CR3BP model.

[0030] FIG. 2 shows an example of a surface of section in a 3-body problem.

[0031] FIG. 3 shows an illustration of the Delaunay Variables.

[0032] FIG. 4 shows an example SOI-resonant orbit and exit/entry to the SOI.

[0033] FIG. 5 shows example landing trajectories of different resonances back-propagated (integrated back) to a surface of section.

[0034] FIG. 6 shows a close-up of the trajectories near the secondary body.

[0035] FIG. 7 shows an example resonant encounter map.

[0036] FIG. 8 shows a blow-up of a section (black rectangle) of FIG. 7.

[0037] FIG. 9 shows three example backward integrations of different resonances.

[0038] FIG. 10 shows an example Swiss Cheese plot.

[0039] FIG. 11 shows states on the Swiss Cheese plot of FIG. 10 for different number of passes.

[0040] FIG. 12 shows an example of a set of initial conditions centered around a location (e.g., landing spot).

[0041] FIG. 13 shows three examples of invariant funnels.

[0042] FIG. 14 shows an example of the initial conditions for a resonant encounter map.

[0043] FIG. 15 shows an example of a hemisphere of initial conditions shaded based on the flight path angle, t.

[0044] FIG. 16 shows a 5:6 resonant ring based on the initial conditions of FIG. 15.

[0045] FIG. 17 shows examples of integrating each resonant ring backwards in time.

[0046] FIG. 18 shows the same resonant encounter map as shown in FIG. 7 but generated with ε=5e-4 and C.sub.0=3:0020. The star shows the nominal landing site.

[0047] FIG. 19 shows a blow-up of a section (black rectangle) of FIG. 18.

[0048] FIG. 20 shows an example 5:6 resonant map

[0049] FIG. 21 shows a resonant ring of FIG. 20 integrated backward through the focus.

[0050] FIG. 22 shows an example of the evolution of the 5:6 resonant encounter map as it is integrated backward through the focus.

[0051] FIG. 23 shows an example of Poincaré sections taken along 9:10 resonant ring trajectories, showing a pinch-point in panel g.

[0052] FIG. 24 shows an example of Poincaré sections taken along 11:12 resonant ring trajectories, showing three gaps where the 10:11 and 11:12 rings lie before and after passing through the focus.

[0053] FIG. 25 shows an example resonant encounter map populated with only SOI-resonant orbits.

[0054] FIG. 26 shows the resonant orbits of FIG. 25.

[0055] FIG. 27 shows a blow-up of a region (thick rectangle) of FIG. 25.

[0056] FIG. 28 shows a blow-up of a region (thin rectangle) of FIG. 25.

[0057] FIG. 29 shows an example of a false SOI-resonant orbit.

[0058] FIG. 30 shows an example of two pseudo-resonant orbits tied with a heteroclinic connection

[0059] FIG. 31 shows the heteroclinic connection of FIG. 30 close to the secondary body.

[0060] FIG. 32 shows an example 9:10 resonant ring trajectories showing a phase offset.

[0061] FIG. 33 shows the trajectories of FIG. 32 close to the secondary body.

[0062] FIG. 34 shows examples of two groups of trajectories taken from FIG. 33.

[0063] FIG. 35 shows trajectories that seem to follow the stable manifolds of a halo orbit before exiting to join the two main groups of FIG. 34.

[0064] FIG. 36 shows an example trajectory that remains in the SOI for multiple revolutions before finally escaping without joining a group of FIG. 34.

[0065] FIG. 37 shows an example of sampling unit vectors for all azimuth angles but constrained in elevation.

[0066] FIG. 38 shows an example of sampling unit vectors constrained within a cone.

[0067] FIG. 39 shows an example of how the points passing through the surface of section become the Poincaré map.

DETAILED DESCRIPTION

[0068] As used herein, a “body” is an object with mass. Typically, herein it would be a planet, moon, other astronomical object, or a spacecraft or other artificial objects.

[0069] As used herein, a “spacecraft” is any vehicle or platform, either manned or unmanned, capable of travelling outside Earth's atmosphere. Examples include satellites, probes, landing modules, orbiters, rover spacecraft, penetrator spacecraft, and cargo transport spacecraft, service spacecraft, colonization spacecraft, spacefaring robots.

[0070] As used herein, an “orbit” is a trajectory of one body around or partially around another body. Examples include a spacecraft around a primary body, a spacecraft around a secondary body, and a secondary body around a primary body.

[0071] A “resonant orbit” is an orbit that has a period approximately commensurate with that of the secondary body, expressed as a ratio between two natural numbers p and q:

[00001]$\begin{array}{cc}\frac{p}{q}=\frac{T_2}{T_{\mathrm{SC}}}\pm ϵ& \left(1\right)\end{array}$

where T.sub.2 is the period of the secondary body, T.sub.SC is the period of the spacecraft, and E is some tolerance (normally ϵ=0.0000000001). For example if p:q=5:6 then the spacecraft completes 5 orbits about the primary in about the same time that the secondary completes 6 orbits. An exterior resonant orbit has p<q, while an interior resonant orbit has p>q.

[0072] A “periodic resonant orbit” is an orbit that returns to the same initial state after some time (the Poincaré section of a periodic orbit is a finite number of points). There are both stable and unstable periodic orbits.

[0073] A “pseudo-resonant orbit” is an orbit which shadows a periodic resonant orbit. Examples are the trajectories on the invariant manifolds of the periodic resonant orbits. A manifold is a high dimensional surface. Invariant manifolds are a tubular structure formed form trajectories that either approach a periodic orbit or depart a periodic orbit. This is a special feature of the Three-Body Problem which does not exist in the Two-Body Problem. Most orbits in this work are pseudo-resonant orbits. With perturbations and navigation errors and using the real ephemeris of the Solar System, there are no “true resonant orbits” to speak of. As time progresses, these orbits may evolve to different values of p and q through heteroclinic connections, especially if they make a close flyby of the secondary body.

[0074] An “SOI-resonant orbit” is a pseudo-resonant orbit that returns to the sphere of influence (SOI) of the secondary body retaining the p:q resonance condition. Flybys and landing orbits are included in this group. Most flyby orbits transition from one resonance to another while in the SOI.

[0075] A “flyby trajectory” is a trajectory that enters and exits the SOI without intersecting the surface of the secondary body.

[0076] A “landing trajectory” is a trajectory that intersects with the surface of the secondary body as it travels. Pseudo-resonant orbits can be landing trajectories given their use in connecting flybys of other small bodies to a landing site.

[0077] A “nominal trajectory” is a trajectory of interest for a particular design problem, e.g. this can be a landing orbit or a flyby orbit. The nominal trajectory is used for computing an invariant funnel or resonant encounter map. The nominal trajectory as a function of time is denoted by x(t). Typically, the landing condition is chosen to be at the periapsis of the orbit. The time at landing or periapsis is t=0, so the state at landing or periapsis is

[00002]$x\left(0\right)=x_0=\left[\begin{array}{c}{r}_0\\ v_0\end{array}\right].$

[0078] As used herein, a “resonant ring” is a region, normally an annulus or ring, on the resonant encounter map where all initial conditions share the same resonance.

[0079] In embodiments described herein, the Circular Restricted Three-Body Problem (CR3BP) dynamic model is used. Standard barycentric normalized coordinates can be used. The mass parameter p is defined as

[00003]$\begin{array}{cc}\mu =\frac{m_2}{m_1+m_2}& \left(2\right)\end{array}$

where m.sub.1 is the mass of the primary body and m.sub.2 is the mass of the secondary body, where m.sub.1 is greater than (or equal to) m.sub.2. This gives a range of p in the range of 0 to 0.5. Some example approximate values are μ(Jupiter/Europa)=2.5280176826×10.sup.−5 and μ(Earth/Moon)=1.2150577033×10.sup.−2.

[0080] For the purposes of this method, the primary and secondary bodies are defined to be in circular orbits about their barycenter with a period normalized to 2π. The distance between them is normalized to 1. The CR3BP can be plotted in the rotating frame, meaning the primary and secondary bodies are fixed on the x-axis at x.sub.1=−μ and x.sub.2=1−μ respectively.

[0081] FIG. 1 illustrates this model with the barycenter at the origin. [NON] refers to non-dimensional units. The primary and secondary bodies lie on the x-axis at −μ and 1−μ respectively. The equilibrium points (Lagrange points) are marked by black x's, with L1 and L2 lying immediately to the left and right of the smaller body. The grey region is the forbidden region, bounded by the zero-velocity curve. The shape of the forbidden region depends on the Jacobi constant, C. Setting μ=0.001 and C=3.035 produces this plot where the L2 gateway is clearly visible.

[0082] The six-dimensional state x of the third body (such as the spacecraft) is composed of the position r and velocity v, where:

[00004]$\begin{array}{cc}r={\left[xyz\right]}^{\top }& \left(3\right)\\ v={\left[\stackrel{.}{x}\stackrel{.}{y}\stackrel{.}{z}\right]}^{\top }& \left(4\right)\\ x=\left[\begin{array}{c}r\\ v\end{array}\right]& \left(5\right)\end{array}$

[0083] Since the mass of the third body, the spacecraft, is negligible compared to that of the primary and secondary, it is considered to have no gravitational effect on either of the large bodies, hence the restricted three-body problem. The equations of motion for a particle in the CR3BP rotating frame are

[00005]$\begin{array}{cc}\stackrel{\mathrm{.Math.}}{x}=-\left(\left(1-\mu \right)\frac{x+\mu }{{.Math.r_1.Math.}^3}+\mu \frac{x-1+\mu }{{.Math.r_2.Math.}^3}\right)+x+2\stackrel{.}{y}& \left(6\right)\\ \stackrel{\mathrm{.Math.}}{y}=-\left(\left(1-\mu \right)\frac{y}{{.Math.r_1.Math.}^3}+\mu \frac{y}{{.Math.r_2.Math.}^3}\right)+y+2\stackrel{.}{x}& \left(7\right)\\ \stackrel{\mathrm{.Math.}}{z}=-\left(\left(1-\mu \right)\frac{z}{{.Math.r_1.Math.}^3}+\mu \frac{z}{{.Math.r_2.Math.}^3}\right)& \left(8\right)\\ r_1={\left[\left(x+\mu \right),y,z\right]}^{\top }& \left(9\right)\\ r_1={\left[\left(x+\mu \right),y,z\right]}^{\top }& \left(10\right)\end{array}$

[0084] The augmented potential is

[00006]$\begin{array}{cc}U\left(r\right)=\frac{\left(x^2+y^2\right)}{2}+\frac{\left(1-\mu \right)}{.Math.r_1.Math.}+\frac{\mu }{.Math.r_2.Math.}& \left(11\right)\end{array}$

and the Jacobi constant is

C(r,v)=2U(r)−v.sup.Tv  (12)

[0085] The forbidden region is defined for a given Jacobi constant C as the area bounded by the zero-velocity surface v.sub.0={r∈.sup.3|2U(r)=C}. There are five equilibrium points in the CR3BP, L1-L5, known as libration points or Lagrange points. The Jacobi constant of a libration point is computed with v=0 (e.g., C.sub.L2=C(r.sub.L2; 0)). As the Jacobi constant decreases below C.sub.L2, the forbidden region shrinks and splits apart, first at L1, then at L2, as shown in FIG. 1. This effectively creates a gateway for spacecraft to approach the secondary body from the exterior or interior regions.

[0086] To reduce the dimensions of the system (for easier analysis), a Poincaré map can be utilized. When many orbits are present, it can be like spaghetti, very difficult to visualize. Poincaré thought of placing a plane cutting across the orbits and study the resulting discrete plot of points which is much easier to analyze. FIG. 39 illustrates the idea, with orbit trajectories (3910) being back integrated to a surface of section (3920) provide points (3930) that form the Poincaré map (3940). The plane is called a “surface of section” and it can be replaced by a surface when it is convenient. More generally, a Poincaré map is a discrete map between two surfaces of section in some high dimensional system. The spatial CR3BP has a 6-dimensional phase space with a 5-dimensional energy surface. The standard Poincaré map reduces the dimensions of this system by 1, hence, it is 4-dimensional. FIG. 2 shows an example of a generalized Poincaré map on which the concept of the Swiss Cheese plot is based. One of the novelties of the Swiss Cheese plot is the reduction in dimension of the standard Poincaré map from 4 dimensions to our 2-dimensional Swiss Cheese plot. The two dimensional “surface of section” S.sub.1 (210) is a half-plane associated with the primary body (220), the surface of section defined by

S.sub.1={(x,y,z)|y=0,x<0}  (13)

with a second surface of section, S.sub.2, being, for landing solutions, the surface of the secondary body (230) (e.g., moon) defined by

S.sub.2={(x,y,z)|(x−1+μ).sup.2+y.sup.2+z.sup.2=R.sub.2.sup.2}  (14)

where R.sub.2 is the radius of the secondary body. In the case where the destination is an orbit around the secondary body or a flyby trajectory around the secondary body, R.sub.2 can be the distance from the center of mass of the secondary body to the planned orbit/flyby nearest approach.

[0087] In embodiments of the invention, the system and method includes forming a Swiss Cheese plot.

[0088] A Swiss Cheese plot starts with a 2D Poincaré section taken at S.sub.1 with Delaunay variables L and g as axes. The plot is used to identify a nominal trajectory to a desired landing site or periapsis location. It can also be used to identify periodic resonant orbits. The 2D map is achieved by using a 2D surface of velocity directions centered at the landing site as the starting conditions for the surface of section. One of the key novelties here is the calculation of the 2D Poincaré section to form the Swiss Cheese plot in the 6D phase space. Whereas the typical Poincaré section in the 6D phase space has to be 4-dimensional, which is extremely difficult to work with due to its high dimension.

[0089] In order to form the Swiss Cheese plot, the values of the Delaunay variables L and g need to be determined for each event state. These can be calculated, for example, as shown below and are represented graphically in FIG. 3. Delaunay variables are an alternate system of orbital elements used to fully describe a state.

[00007]$\begin{array}{cc}r=\sqrt{x^2+y^2}& \left(15\right)\\ \theta ={\tan }^{-1}\left(\frac{y}{x}\right)& \left(16\right)\\ \stackrel{.}{r}=\stackrel{.}{x}\left(\frac{x}{r}\right)+\stackrel{.}{y}\left(\frac{y}{r}\right)& \left(17\right)\\ \stackrel{.}{\theta }=-\stackrel{.}{x}\left(\frac{y}{r^2}\right)+\stackrel{.}{y}\left(\frac{x}{r^2}\right)& \left(18\right)\\ G=r^2\left(1+\stackrel{.}{\theta }\right)& \left(19\right)\\ L={\left(-\left(\frac{G^2}{r^2}\right)+\left(\frac{2}{r}\right)-r^2\right)}^{-1/2}& \left(20\right)\\ \stackrel{\_}{g}=\theta -{\cos }^{-1}\left(\frac{\frac{G^2}{r}-1}{e}\right)& \left(21\right)\end{array}$

[0090] Note that in the normalized CR3BP, L=√{square root over (a)} and g=ω−t.

[0091] In embodiments of the system and method, forming a Swiss Cheese plot includes selecting initial conditions. A Swiss Cheese plot is formed by plotting back-integrated trajectories intersecting the surface of section against a plot of a standard Poincaré map formed from the Delaunay variables.

[0092] Define a set of initial conditions at the point of interest to find an orbit that reaches a given point on or near the secondary body. These 6-dimensional initial conditions share the same position and Jacobi constant. A landing site or periapsis location, r.sub.0, and a Jacobi constant, C.sub.0 are selected. The Jacobi constant must be less than that of the L2 libration point, C.sub.L2, to allow for both interior and exterior resonant orbits to appear. If the Jacobi constant is greater than C.sub.L2, then the L2 gateway is closed. With r.sub.0, and C.sub.0 chosen, v.sub.0 is constrained to have magnitude:

[00008]$\begin{array}{cc}.Math.v_0.Math.=\sqrt{x_0^2+y_0^2+\frac{2\left(1-\mu \right)}{.Math.r_1.Math.}+\frac{2\mu }{.Math.r_2.Math.}-C_0}& \left(22\right)\end{array}$

Therefore, there are two remaining degrees of freedom to define the velocity direction.

[0093] Sample N.sub.v unit vectors originating at the location r.sub.0 with different directions in those two degrees of freedom. The sampling is done with different methods depending on if this is targeting a landing location or a periapsis. For example, when targeting a landing location, constrain the flight path angle, ϕ, to be within some limits (e.g., between 0 to 5 degrees) and sample directions using Archimedes' theorem. If targeting a periapsis, constrain the velocity direction to a disc perpendicular to r.sub.2 and sample uniformly along the circumference. FIG. 37 shows an example of sampling unit vectors for all azimuth angles but constrained in elevation to −5° to 0°. FIG. 38 shows an example of sampling unit vectors to be within a cone with an opening angle α=5°.

[0094] FIG. 10 shows an example Swiss Cheese plot. Every point corresponds to an impact trajectory that lands on the surface of the secondary body at the south pole with flight path angle ϕ<5 degrees and passes through the surface of section. The gaps correspond to stable quasi-periodic resonant orbits surrounding the stable periodic orbit at the center. The p:q resonances at the top correspond to the L variable indicated by the vertical dashed lines. Stable orbits are usually avoided, as they make changing trajectory more difficult.

[0095] The initial conditions are integrated backwards in time and the state at each pass (a number (e.g., “k”) of iterations of intersections with the surface of section) is recorded to build the plot. The k-th iterate is plotted such that it can be distinguished from the standard plot (the plot points from the Delaunay variables). For example, this can be done graphically by plotting the k-th iterations in one color or tone and the Delaunay variable derived points in another color or tone (e.g., black points against grey points). A change of coordinates can be used to find exact events (see e.g., Henon [13]). Programs on parallel computers and GPUs can be used to speed up the integration, allowing for tens of thousands of trajectories to be integrated in parallel. As trajectories pass through 51 for the first time, smooth curves are shown. With each subsequent pass, the plot looks more chaotic. See e.g., FIG. 11.

[0096] The events can be converted to the Delaunay variables L and g. Planar Delaunay variables can be calculated in a rotating frame [7]. As stated above, L=√{square root over (a)}, and since a is a function of period T, each L will correspond to a specific resonance value as shown below:

[00009]$\begin{array}{cc}\frac{p}{q}=L^{-3}& \left(23\right)\end{array}$

[0097] Note that the unstable pseudo-resonant orbits where p is odd are at apoapsis in the Poincaré section, while those where p is even are at periapsis. This alternating effect can be seen in the gaps in the map, which alternate being centered around g=0 and g=π radians. The gaps correspond to stable quasi-periodic orbits with the stable periodic orbits at the center of the gap. The unstable periodic orbit for each resonance is located exactly π radians above.

[0098] From the Swiss Cheese plot, a nominal trajectory can be selected from a pseudo-resonant orbit of a desired resonance. In some embodiments, this is picked from a point that came during a first pass in the Swiss Cheese plot. The case for landing on the North Pole of Enceldus is used here as an example on how to use the Swiss Cheese Plot to find a nominal trajectory, in this case, a resonant trajectory landing at the North Pole. FIG. 37 illustrates how the initial conditions are determined. In this example the landing position is the North Pole, P, of Enceladus. Select a Jacobi constant, C, for this case. Then select n velocities V1, V2, . . . Vn tangent to P uniformly distributed in all directions. These states (P, V1), (P, V2) . . . (P, Vn) all have the same Jacobi constant C. Integrate these initial states backwards until they intersect the surface of section to form the black dots on the Poincaré Map of FIG. 11. The plot for Pass #1 shows a circle (1105) around the line of black dots where the 13:14 resonance can be determined. This is because the background of grey dots are the Poincaré Map of the resonances. Each resonant orbit is indicated by the vertical line of holes. The desired unstable resonant orbit is located exactly in the middle between the vertical holes. The circle in plot Pass #1 captures the location of the 13:14 resonance. Although this example describes this method for landing on the North Pole of Enceladus, this method is equally applicable for any point on the surface of Enceladus. It is also applicable for any point on the surface of a sphere centered at Enceladus. The case for the sphere is for the design of a close flyby of Enceladus where capture into orbit around Enceladus is also possible.

[0099] FIG. 5 shows examples of nominal landing trajectories back-propagated to S.sub.1. The orbits are classified p:q (e.g. 6:7) by their period with respect to the primary body. FIG. 6 shows a close-up of the secondary body (here, Europa) with the trajectories converging on the landing site.

[0100] In embodiments of the system and method, an invariant funnel is created based on the nominal trajectory.

[0101] An invariant funnel is a set of trajectories that converge to a nominal trajectory as t increases. In some embodiments, all trajectories of the funnel share the same Jacobi constant and have parallel velocities at the landing site or periapsis.

[0102] Beginning with a nominal pseudo-resonant trajectory (flyby or landing), perhaps found using the Swiss Cheese plot as shown above, the Jacobi constant C.sub.0 is determined from a state along the trajectory. The landing/periapsis state (r.sub.0, v.sub.0) can be used to avoid numerical errors.

[0103] A set of N.sub.r positions are sampled in a ring of radius R around the periapsis/landing site. The sampled points can be within the ring (circle) or just at the periphery. The area can be a circle, but other shapes (e.g., ellipse, square, etc.) can also be used. See e.g., FIG. 12. To each of these position vectors assign a velocity vector parallel to the original velocity vector. Constrain the magnitude of each velocity vector so that each state has the same Jacobi constant (see eq. 22).

[0104] As long as R is small enough, the initial conditions will trace out an invariant funnel about the original trajectory when integrated backward. The funneling effect is particularly strong near L1 and L2. Any initial condition inside the ring will stay within the funnel, unless R is too large. Then the funnel will diverge into chaos, with some trajectories exiting through the L1 gateway and others through L2. The radius where this begins to happen is R.sub.max which must be determined numerically. Examples of invariant funnels for various systems are shown in FIG. 13.

[0105] In embodiments of the system and method, a resonant encounter map is created based on the nominal trajectory.

[0106] A resonant encounter map is a mapping on the entire surface of the secondary body that shows all of the possible resonances of landing trajectories. The landing states all share the same Jacobi constant and velocity direction as some nominal trajectory. The landing trajectories are not necessarily SOI-resonant. An example resonant encounter map is shown in FIGS. 7 and 8 (FIG. 8 showing a blow-up of a portion of FIG. 7).

[0107] A “resonant ring” is a region, normally an annulus or ring, on the resonant encounter map where all initial conditions share the same resonance. For example, FIG. 8 shows a thin region of 5:6 pseudo-resonant orbits (710). The 11:12 region (720) near the center can also be called a resonant ring even though it is topologically not an annulus.

[0108] Similar to the invariant funnel, begin with a nominal landing trajectory with landing state (r.sub.0, v.sub.0) and corresponding Jacobi constant C.sub.0. Then generate a set of position vectors n for i=1 N covering the entire surface of the sphere. FIG. 14 shows an example of initial conditions for a resonant encounter map, with r.sub.0 being the nominal trajectory and n being the sampled position vector locations. Deserno's method [14] can be used to distribute the points nearly uniformly across the sphere. Other methods, such as the Archimedes method, can also be used—some producing a random distribution of points across the sphere. The goal is to use some method to produce a nearly uniform distribution of points on the sphere to discretize the sphere for numerical calculations. Next, assign a velocity to each point with the same direction as v.sub.0, but scaled as to give each state the same Jacobi constant (eq. 22). Almost half of the initial conditions are infeasible as their velocity points out of the surface of the secondary. Prune these by asserting that (r.sub.i−r.sub.2).sup.Tv.sub.i<0. As N goes to infinity, the initial conditions for the resonant encounter map can be written in set notation as

X={(r.sub.i,v.sub.i)∈.sup.3×.sup.3|∥r.sub.i−r.sub.2∥=R.sub.2,v.sub.i×v.sub.0=0,(r.sub.i−r.sub.2).sup.Tv.sub.i<0,C.sub.i(r.sub.i,v.sub.i)=C.sub.0}  (24)

[0109] The flight path angle, ϕ, for the initial conditions varies widely even across a single resonant ring (see e.g., FIGS. 15 and 16). In this example shown in FIG. 16, the nominal trajectory is 5:6 resonant and has a shallow flight path angle ϕ.sub.0≈−10°. The entire 5:6 resonant ring varies by almost 35° across the ring.

[0110] Just as with the invariant funnel, integrate each of these initial conditions backward in time (t<0) until they reach S1 (the XZ plane). At that point, compute the Delaunay variable L (see eq. 20) and the approximate resonance p:q (see eq. 25).

[00010]$\begin{array}{cc}\frac{p}{q}=L^{-3}+ϵ& \left(25\right)\end{array}$

[0111] FIG. 15 shows an example hemisphere of initial conditions on the surface of a secondary body shaded according to flight path angle. Since all initial velocities are parallel, the flight path angle varies from 0° to −90° across the hemisphere. FIG. 16 shows that, for FIG. 15, the 5:6 resonant ring does not have a uniform flight path angle, rather it varies from −10° to −44°.

[0112] FIG. 17 shows examples of these backward integrated trajectories near the secondary body, grouped by resonance. Each shade corresponds to a different resonant ring.

[0113] Plot the initial conditions on the surface of the sphere, differentiating them by resonance (e.g., color coding). Decreasing c from 0.01, a thinning of the ringed regions corresponding to each resonance is observed, as shown in FIGS. 18 and 19 (FIG. 19 being a blow-up of a portion of FIG. 18).

[0114] Every resonant ring passes through a region of convergence, or “focus”, close to the secondary body as it is integrated backward. See e.g., FIGS. 17 and 21 (FIG. 21 being the 5:6 trajectories of FIG. 20).

[0115] As used herein, a “focus” is a region of convergence that appears when integrating a resonant ring backward. FIG. 9 shows examples of foci for resonant rings integrated backward: 5:6 (left), 9:10 (center), and 12:13 (right).

[0116] Taking Poincaré sections along the trajectories, one can visualize how the topology of the resonant ring changes. FIG. 22 shows an example of a set of Poincaré sections for a 5:6 resonant ring. As the trajectories pass through the focus, they pass through each other in configuration space, effectively turning the ring inside out.

[0117] FIG. 23 shows an example of Poincaré sections taken along the 9:10 resonant ring trajectories. The ring can be seen twisting over itself, leaving a pinch point in the top left of panel g.

[0118] FIG. 24 shows an example of Poincaré sections taken along the 11:12 resonant ring trajectories. The three gaps where the 10:11 and 11:12 rings lie can be seen clearly before and after passing through the focus.

[0119] In some embodiments, SOI-resonant trajectories are used from the map. FIG. 4 shows an example of a SOI-resonant orbit shown in a rotating frame.

[0120] SOI-resonant trajectories are useful for missions that incorporate a flyby of the target body to pump down energy before landing. FIGS. 25 to 28 show examples of how the resonant encounter map becomes quite sparse when only accounting for SOI-resonant trajectories. Only the outer ring, where the flight path angle is most shallow, is SOI-resonant. This is convenient, since it is much easier to land with a shallow flight path angle.

[0121] When computing SOI-resonant orbits, check for a crossing of S.sub.1 between checking for exits and entries of the SOI. Otherwise, it is easy to include false SOI-resonant orbits like the example in FIG. 29, which exits and enters the SOI multiple times before finally leaving through the L1 gateway and never returning.

[0122] Most SOI-resonant orbits switch from one resonance to another during the flyby. This is illustrated in FIGS. 30 and 31, where there is an example of a 5:6 change into an 11:12 during a close encounter with Europa. This is called a heteroclinic connection.

[0123] In some embodiments, groups of trajectories are determined.

[0124] When integrating the trajectories of the resonant encounter map, each resonant ring tends to branch off into groups. For example, FIGS. 32 to 36 show how most of the 9:10 resonant ring trajectories can be divided into two mainstream groups. There are some orbits that seem to follow a halo-type trajectory before eventually making their way to one of the mainstream groups. Only a few trajectories don't follow the mainstream groups, and they tend to remain stuck around the secondary body for several orbits before finally escaping.

[0125] A number of embodiments of the disclosure have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the present disclosure. Accordingly, other embodiments are within the scope of the following claims.

[0126] The examples set forth above are provided to those of ordinary skill in the art as a complete disclosure and description of how to make and use the embodiments of the disclosure, and are not intended to limit the scope of what the inventor/inventors regard as their disclosure.

[0127] Modifications of the above-described modes for carrying out the methods and systems herein disclosed that are obvious to persons of skill in the art are intended to be within the scope of the following claims. All patents and publications mentioned in the specification are indicative of the levels of skill of those skilled in the art to which the disclosure pertains. All references cited in this disclosure are incorporated by reference to the same extent as if each reference had been incorporated by reference in its entirety individually.

[0128] It is to be understood that the disclosure is not limited to particular methods or systems, which can, of course, vary. It is also to be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting. As used in this specification and the appended claims, the singular forms “a,” “an,” and “the” include plural referents unless the content clearly dictates otherwise. The term “plurality” includes two or more referents unless the content clearly dictates otherwise. Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the disclosure pertains. All documents referenced in this specification are incorporated by reference in their entirety.

[0129] List of References and Related Matter, all of which are incorporated by reference in their entirety: [0130] [1] J. T. Blanchard, M. W. Lo, D. Landau, and B. D. Anderson, “Invariant Funnels for Resonant Landing Orbits,” AAS/AIAA Spaceflight Mechanics Meeting, 2021. [0131] [2] J. T. Blanchard, B. D. Anderson, M. W. Lo, and S. Close, “Low energy capture into high inclination orbits for ocean worlds missions,” AAS/AIAA Astrodynamics Specialist Conference, 2020. [0132] [3] V. Szebehely, “Theory of Orbits—The Restricted Three Body Problem,” 1967. [0133] [4] H. Pollard, Mathematical Introduction to Celestial Mechanics. Pollard, 1966. [0134] [5] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer, 2 ed., 2003, 10.1007/b97481. [0135] [6] T. S. Parker and L. O. Chua, “Practical Numerical Algorithms for Chaotic Systems,” Mathematics of Computation, Vol. 56, 4 2006, p. 377, 10.2307/2008550. [0136] [7] W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, Dynamical Systems, the Three-Body Problem and Space Mission Design. Koon, Lo, Marsden, Ross, 2006, 10.1142/9789812792617 0222. [0137] [8] T. M. Vaquero Escribano, Poincaré Sections and Resonant Orbits in the Restricted Three-Body Problem. PhD thesis, Purdue University, 2010. [0138] [9] E. Barrabés and G. Gómez, “A note on second species solutions generated from p-q resonant orbits,” Celestial Mechanics and Dynamical Astronomy, Vol. 88, No. 3, 2004, pp. 229-244, 10.1023/B:CELE.0000017163.37356.50. [0139] [10] E. Barrabés and G. Gómez, “Spatial p-q resonant orbits of the RTBP,” Celestial Mechanics and Dynamical Astronomy, Vol. 84, No. 4, 2002, pp. 387-407, 10.1023/A:1021137127909. [0140] [11] J. Font, A. Nunes, and C. Sim'o, “Consecutive quasi-collisions in the planar circular RTBP,” Nonlinearity, Vol. 15, 1 2002, pp. 115-142, 10.1088/0951-7715/15/1/306. [0141] [12] M.-z. Shao and N. I. Badler, “Spherical Sampling by Archimedes' Theorem,” Technical Reports (CIS), Univ. of Pennsylvania, No. January 1996, p. 11. [0142] [13] M. Henon, “On the Numerical Computation of Poincaré Maps,” Physica 5D, 1982, pp. 412-414. [0143] [14] M. Deserno, “How to generate equidistributed points on the surface of a sphere,” 2004, p. 55128. [0144] [15] R. L. Anderson and M. W. Lo, “Flyby design using heteroclinic and homoclinic connections of unstable resonant orbits,” AAS/AIAA Spaceflight Mechanics Meeting, Vol. 140, 2011. [0145] [16] D. L. Richardson, “Analytic construction of periodic orbits about the collinear points,” Celestial Mechanics, Vol. 22, No. 3, 1980, pp. 241-253, 10.1007/BF01229511.