CR3BP module: comments and suggestions for expansion

[MiguelAvillez]

Hi all,

Just leaving some comments on the current state of the CR3BP module (this branch) and some suggestions for future developments.

Initial guess:

• The current code generates the initial guess for Lagrange point orbits (LPOs) using a Lindstedt-Poincaré method, which requires the user to provide the initial guess for the amplitudes of the orbits. At the moment, the amplitudes are hard-coded for the Earth-Moon system without the possibility for the user to select other amplitudes, so it would be nice if the user was able to select them. I'm not sure how much these amplitudes vary between systems, but perhaps it is possible to find some paper/book with recommended amplitudes for different systems.
• I don't know how sensitive the Lindstedt-Poincaré method is to poor initial guesses, but a more robust (in the sense that it would always work) way to generate the LPOs would be to compute the families as they emanate from the Lagrange points. For example, the horizontal and vertical Lyapunov families always bifurcate from the colinear Lagrange points, so it is quite easy to compute them (see Bardem and Howell, 1998; Meyer, 2017, section 9.2). The halo family then bifurcates from either the vertical or horizontal family; however, the family from which it bifurcates depends on the dynamical model (Scheeres et al., 2019), so it might be difficult to setup a code to generate the halo from with this method without any kind of user input.
• Currently, the user also needs to be able to manually generate two periodic orbits; this can be easily reduced to just one by modifying the continuation algorithm; see below.
• Finally, it would be very nice if it was possible to generate other families of orbits, e.g. QSOs, DROs (QSOs and DROs correspond to the same family; I'm calling QSOs to the orbits close to the secondary, and DROs to the ones further apart from the secondary), DPOs... and whatever else it out there :). Initial guesses for QSOs can be easily generated according to an ellipse (Henon, 1969). DROs can be computed by starting from a low-altitude QSO and continuing it to high altitude; alternatively, Dominik Stahl in his thesis focused on modeling DROs using semi-analytical methods, so might be work checking it to see whether he offers some more direct option. DPOs can also be generated according to Henon, 1969; Mingotti, 2010 generates them using that initial guess.

Differential correction:

• Currently done with single shooting and taking advantage of CR3BP symmetries.
• Would be nice to also have the option to compute orbits in a fully generic formulation (i.e. without symmetries), especially if computing orbits with irregular gravity fields (e.g. SH or polyhedron).
• Would also be cool to have the option of using multiple shooting, especially for computing multi-periodic orbits.

Continuation algorithm

• Regarding the current code, Koen Langenmeijer says it uses pseudo-arclength continuation, though I would consider that a very generous use of the term.
• The current code approximates the tangent to the family as the difference between two previously-known orbits. It doesn't really make a lot of sense to that, since the actual tangent to the family can be very easily calculated, as the variational equations are being propagated. Also, calculating the true tangent would mean only one orbit is necessary for the initial guess instead of two.
• The current code applies the same continuation step size to the orbit initial state and orbit period, which is rather questionable. It would make more sense to integrate the orbit state and period into a single vector and apply the step size to that vector.
• The current code doesn't really fix the direction of continuation of the differential corrector along the family, meaning that in principle, the differential corrector might jump forwards and backwards while continuing the family. This direction of continuation would be generally fixed using the natural parameter for natural-parameter continuation or a pseudo-arclength constraint for (you guessed it) pseudo-arclength continuation.
• Rant over. For a description of pseudo-arclength continuation see Doedel, 2003; Doedel, 1991.

Continuation step-size control

• Currently a fixed continuation step-size is used, which needs to be manually selected. It would be very pleasant if it were possible to use a variable step-size. For example, Allgower and Geog (2003) (side note: this is a really good book regarding continuation-related topics) present one step control method based on an error model for the differential corrector and one which is very similar to the step control used in the numerical integrators. Both seem good options; I used the former in my thesis because its implementation is very simple, it worked quite well.

Computation of manifolds

• I might be mistaken, but I think it isn't checked anywhere whether the orbits are unstable before trying to compute their manifolds (it should be, as otherwise there are no manifolds).
• Regarding the propagation of the manifolds, it should be possible to select which manifolds are to be propagated and to select the termination conditions to use in the propagation of the manifolds.

Further suggestions:

• Allow easily selecting the gravity model to use in the CR3BP. For example, secondary modeled by SH or polyhedron when studying the motion near binary asteroids or small moons.
• Include homotopy procedure. To deal with convergence issues when using models with irregular gravity fields.

Ok, I think I'm gonna stop here... already typed enough 😊
I'll add a link to my thesis in a few weeks, has some more useful references and information on continuation, bifurcation and branch switching