# Designing the low-energy lunar transfers trajectories which pass in the vicinity of the libration points of the Earth - Moon system. Part 2. Algorithm and numerical analysis

**Authors:**Konstantinov M.S.^{1}, Thant A.M.^{1}-
**Affiliations:**- Moscow Aviation Institute (National Research University)

**Issue:**Vol 24, No 2 (2023)**Pages:**111-120**Section:**Articles**URL:**https://journals.rudn.ru/engineering-researches/article/view/35136**DOI:**https://doi.org/10.22363/2312-8143-2023-24-2-111-120**EDN:**https://elibrary.ru/CQFLUD

Cite item

## Full Text

## Abstract

An algorithm for designing a low-energy lunar flight trajectory is presented. It is based on the assumption that low-energy flight trajectories pass through the vicinity of one of the collinear libration points of the Earth - Moon system (L1 or L2). It is also assumed that at the moment of spacecraft flight in the vicinity of the libration point, the elements of the osculating geocentric orbit of the spacecraft are close to the elements of the osculating geocentric orbit of the libration point itself. The results of a numerical analysis of the obtained low-energy lunar flight trajectory are presented. It is shown that the use of such a trajectory makes it possible to reduce the deceleration impulse of the velocity during the transition to a low lunar orbit to a value of 638 m/s (in the traditional flight scheme, this impulse turns out to be more than 800 m/s). The influence of solar gravitational disturbances on the flight trajectory is analyzed. It is demonstrated that these perturbations ensure the approach of the spacecraft to the Moon with a negative selenocentric energy constant and contribute to the temporary capture of the spacecraft by the Moon. The influence of the terrestrial gravitational perturbation on the circumlunar part of the trajectory is studied. It is displayed that on the trajectory found this perturbation effectively reduces the selenocentric velocity of the spacecraft. The conditions for spacecraft flight in the vicinity of the libration point are considered.

## Full Text

Introduction[2] The problem of analyzing low-energy trajectories to the Moon is considered in many publications [1-11]. The first part of the article [12] describes a new method for designing low-energy trajectories for a flight to the Moon with the spacecraft (SC) insertion into a low lunar orbit. The method is based on the assumption that a low-energy trajectory can be obtained by using as an initial approximation a trajectory passing through the vicinity of one of the collinear libration points of the Earth-Moon system (L1 or L2). Additionally, it is assumed that on this trajectory, at the moment when the spacecraft approaches the libration point, the elements of the geocentric osculating orbit of the spacecraft are close to the corresponding elements of the geocentric osculating orbit of the libration point. Thus, the method assumes a narrowing of the space of analyzed trajectories. The authors do not claim that there are no low-energy lunar flights that cannot be carried out obtained using a smooth continuation of the trajectories considered in the paper as an initial approximation. But we believe that many (if not most) low-energy lunar trajectories can be obtained using a trajectory passing through the vicinity of the considered libration points as an initial approximation. The problem statement assumes that the impulse trajectory of the flight from low Earth’s orbit (LEO) to low Moon orbit (LMO) is being analyzed. It is assumed that the altitude and inclination of LEO and the altitude of LMO are known. The selectable (optimized) characteristics of the flight trajectory are: the date of start (Tst) at the analyzed given epoch, the longitude of the ascending node of the LEO (Ω), the latitude argument of the starting point (uo), the magnitude of the accelerating velocity impulse at the start (∆V1), the flight time to the target orbit of the artificial satellite of the moon (tp), the magnitude and direction of the braking velocity impulse at the end point of the flight trajectory to the moon (∆Vbr). The listed characteristics should be chosen so as to: 1) ensure the solution of the transport problem (insertion into the target LMO) and 2) the cost of solving the transport problem is minimal. The considered optimization criterion is either the summary velocity impulse or the magnitude of the velocity impulse, which ensures the insertion of the spacecraft into LMO when SC approaching this orbit (∆Vbr). To analyze the passage of a spacecraft in the vicinity of the libration point, the sum of three non-negative quantities (three distances) is considered: where the first term ΔrL is the SC distance from the libration point (it is found as the difference between the geocentric vectors of the SC and the libration point); and - radius of perigee and apogee of the osculating geocentric orbit of the SC; and - perigee and apogee radius of the osculating geocentric orbit of the libration point. The J-function depend on the four parameters of the flight pattern Tst, Ω, uo, ra, which determine the conditions for the motion of the SC after its launch from LEO, and the current time of motion of the SC t: J(Tst, Ω, uo, ra, t). On each flight trajectory, there is a time t1 when J is minimal. Let us define this minimum value as I and call it the total miss of the libration point: (1) The paper describes the developed algorithm and analyzes the numerical results of the obtained low-energy lunar trajectory. 1. Algorithm for designing a low-energy flight trajectory to the circumlunar orbit The developed method for designing low-energy flight trajectories assumes the following sequence of operations. At the first stage, such parameters of the flight pattern Tst, ra, are found, which ensure that the SC enters the vicinity of the libration points and the minimum total miss of this point. More specifically, the task of the stage is formulated as follows: in the space of two of the listed parameters Tst, ra, find such a set of them that minimizes the total miss of the libration point I. In this case, the values of the parameters Ω and uo do not vary and are taken equal to the values described in Section 3.1. The enumeration of, ra is used. The parameters changed with a fairly small step. The launch date increment is one hour; parameter ra increment is 5 thousand km. The launch date change range is one year for the considered launch epoch. Range of parameter ra is 1-1.5 million km. For each pair of values of these parameters, a system of differential equations is integrated that describes the geocentric motion of the SC in the restricted 4-body problem. There is such moment of time t1 when the sum J is minimal. An analysis of the dependence of I as a function of Tst and ra makes it possible to choose a relatively small number of start dates for the epoch under consideration, using which it is possible to fly over the vicinity of the libration point and enter the vicinity of the Moon. These dates and the values of the ra parameter for them are considered as an initial approximation when searching for low-energy flight trajectories. At the next 2nd stage, I is considered as a function of four variables Tst, ra, Ω, uo. The unconditional minimum I is found as a function of these variables (the method of local search is used) and the time when the SC hits the vicinity of the libration point t1 corresponding to this minimum. Further, the SC trajectory is considered to consist of geocentric and selenocentric sections. At the point of transition from the geocentric to the selenocentric section, an intermediate velocity impulse is introduced into consideration. In this case, the space of the chosen parameters of the flight pattern increases by 5 units (the magnitude of the velocity impulse ΔVc, its declination α and right ascension δ, the time (date) of this velocity impulse tgeo, and the time of movement in the selenocentric segment t2). Thus, the condition for the fulfillment of the transport problem [12] is transformed into the condition: (2) At the next 3rd stage, such ΔVc, α, δ and t2 are determined, at which the SC flies up to the Moon at a distance equal to the height of the circumlunar orbit. The remaining parameters of the flight pattern (the function arguments in (2)) do not vary, while tgeo is assumed to be equal to t1. The next three stages of the developed algorithm are based on the gradient projection method, which ensures the constant achievement of a given approach altitude to the Moon. The optimization criterion is the sum of the values of the intermediate velocity impulse and the velocity impulse during the transition to a circumlunar orbit. The stages differ in the number of optimized parameters of the flight pattern. From stage to stage, this number increases from 4 (ΔVc, α, δ and t2) to 5 (tgeo, ΔVc, α, δ and t2) and finally to 9 (Tst, ra, Ω, uo, tgeo, ΔVc, α, δ and t2). At the last stage of the algorithm, the value of the intermediate velocity impulse is excluded from the number of the selected parameters of the flight scheme. This value iteratively decreases. The fulfilment of the condition (1) is ensured by the choice of the remaining arguments of this function. In this case, a situation is possible, in which it is not possible to bring the intermediate velocity impulse to a zero value. Such a solution can be quite good from the point of view of the summary velocity impulse. That is, if the sum of the intermediate impulse and the impulse that transfers the SC to a LMO is sufficiently small, then the flight trajectory can be considered low-energy. 2. Numerical analysis of the flight trajectory to the Moon As an example, the problem of a flight to a LMO with a height of 100 km from a LEO with a altitude of 200 km and an inclination of 51.6о is considered. The start date is 2024. The position of the plane of the LMO is not fixed. The situation in which the libration point L2 is passed is analyzed. Figure 1 shows the level lines of the total miss of libration point I as a function of the start day (X-axis, the first 2920 hours of 2024, from January 1st to May 1st are considered) and ra (Y-axis). The level lines are shown, on which I is less than 300 thousand km. An analysis of the figure shows that for the considered range of launch dates in 2024, there are several areas of launch dates, using which it is possible: 1) to ensure that the SC enters the vicinity of the libration point; 2) to ensure that the shape and size of the osculating geocentric orbit of the SC are close to the shape and size of the osculating orbit of the libration point. It is proposed to explore each of these areas for the possibility of implementing a low-energy lunar flight. In particular, it is possible to provide a small total miss of the libration point in the range of start dates from March 23 to April 15 of the year under consideration. The minimum value of the total miss turned out to be 38.6 thousand km. Such a miss is obtained if we choose March 29, 2024 as the launch date (2133 hours of this year) and ra equal to 11.5 million km. Fig 1.tiff Figure 1. The level lines of the total miss of the libration point I on the plane: launch date (X-axis, hours of 2024) - ra (Y-axis, million km) The table shows the main characteristics of one of the obtained low-energy flight trajectories. The velocity impulse during the transition to the LMO turned out to be 638.1 m/s. The velocity impulse during the transition to the circumlunar orbit has a small radial component (3.15 m/s). The transverse component of the velocity impulse is negative (-638.1 m/s). Main characteristics of the low-energy transfer trajectory Characteristic Value Launch date 5 April 2024 Julian launch date 2460405.8865 Intermediate geocentric orbit apogee radius (ra), thousand km 1269.6106 The value of the first velocity impulse, m/s 3197.702 The ascending node longitude of LEO Ω -0.075797o Perigee argument of intermediate geocentric orbit uo -12.403506o Total flight time, days 87.52329 Maximum distance of the SC from the Earth, thousand km 1735.6 Minimum SC distance from the libration point L2 on the flight trajectory, thousand km 18.620 Selenocentric orbital energy constant at the moment of maximum approach of the SC to the libration point, km2/s2 -0.138 The magnitude of the velocity impulse during the SC transition to LMO ΔVbr, m/s 638.09556 Radial component of velocity impulse ΔVbr, m/s 3.14854 Transverse component of the velocity impulse ΔVbr, m/s -638.08779 Thus, the use of such a trajectory makes it possible to reduce the deceleration impulse of the velocity during the transition of the spacecraft to a low circumlunar orbit to the value of 638 m/s. Let us pay attention to the fact that in the traditional scheme of flight to the Moon with access to a low circumlunar orbit this impulse turns out to be more than 800 m/s. That is, the gain in the velocity impulse turns out to be very large (at least it is 140 m/s). 2.1. Characteristics of a low-energy trajectory Figure 2 shows the projections of the geocentric trajectory of the flight to the Moon on the x-y plane of the ecliptic and on the z-y plane. The dotted line shows the projections of the geocentric Moon’s orbit. In the scale adopted in the figure, the trajectory begins virtually from a point with zero coordinates and ends at a point in the Moon’s orbit. The black diamond shows the position of the libration point L2 at the moment of the maximum approach of the SC to this point. The maximum distance of the SC from the Earth occurs on the 37.064th day of the flight and is equal to 1.408 million km. It should also be noted that the radius of the apogee of the intermediate orbit is 1.270 million km. That is, solar gravitational disturbances ensured an increase in the SC distance from the Earth to the region where these disturbances are large. The SC stays in this region for a long time, which contributes to a large deformation of the geocentric orbit by the solar gravitational disturbance. a b Figure 2. Flight trajectory in the geocentric ecliptic coordinate system, distance unit is 1 million km: a - projection onto the x-y plane of the ecliptic; b - projection onto the z-y plane The SC flight time from the point of greatest distance from the Earth to the vicinity of the libration point is 48.9 days. At this time interval, the declination of the geocentric radius vector of the SC with respect to the plane of the ecliptic is very small (it varies in a narrow range from +4о to -7о). Figure 3 shows the projections of perturbing gravitational accelerations in the study of the geocentric trajectory of the SC. The time interval from the launch of the SC from a low Earth orbit to the moment of approach of the SC to the libration point (85.93 days) is considered. The perturbing solar acceleration ΦSun is analyzed in the left figure, and the lunar perturbing acceleration ΦMoon is analyzed in the right figure. The thick lines show the projections of perturbing accelerations on the direction of the geocentric velocity of the SC (ΦSun_V and ΦMoon_V). Thin solid lines show transversal projections of perturbing accelerations (ΦSun_T and ΦMoon_T). Dashed and dash-dotted lines show the radial (ΦSun_R and ΦMoon_R) and normal (ΦSun_N and ΦMoon_N) components of the perturbing accelerations. Fig 3.tiff a b Figure 3. Projections of perturbing solar gravitational acceleration ΦSun, mm/s2 (a) and projections of perturbing lunar gravitational acceleration ΦMoon, mm/s2 (b) as the functions of flight time (day); the geocentric part of the trajectory It can be seen that significant part of the considered part of the trajectory (0-71.7 days) the projection of the solar gravitational acceleration in the velocity direction ΦSun_V is positive. The maximum value of this acceleration (0.078 mm/s2) is reached on the 35.15th day of the flight, when the SC is at a great distance from the Earth (1.405 million km). Solar gravitational acceleration actively increases the geocentric velocity of the SC, increasing the radii of the apsidal points of the geocentric osculating orbit of the SC. The transversal component of the perturbing acceleration ΦSun is also positive over a long-time interval. This contributes to an increase in the semilatus rectum of the osculating geocentric orbit of the SC. The radial component of the perturbing acceleration ΦSun is somewhat smaller than the transversal one, but is positive over a longer time interval. The normal component of the perturbing solar acceleration with respect to the other components is quite small. The projections of the perturbing lunar acceleration (right figure) on a significant part of the considered trajectory are sinusoidal. Because of this, they do not create significant perturbations of the elements of the SC's geocentric orbit. As the SC approaches the vicinity of the libration point, perturbing lunar acceleration becomes very large. The projection of this acceleration onto the velocity direction reaches a value equal to -0.3 mm/s2. At this moment, the transversal, radial, and even normal components of the perturbing acceleration are also large. Figure 4 shows the projections of perturbing gravitational accelerations in the study of the selenocentric trajectory of the SC. The time interval from 76 days of flight to the moment the SC enters the LMO is considered. The left figure analyzes the disturbing terrestrial acceleration ΦEarth, and the right figure analyzes the solar disturbing acceleration ΦSun. The same notation is used as in the previous figure (Figure 3). Note the following properties of the given characteristics. Solar perturbing gravitational accelerations are less than perturbing Earth's accelerations by two orders of magnitude. They have very little effect on the trajectory of the SC. The projection of the perturbing terrestrial acceleration on the direction of the selenocentric velocity (thick line) is negative over the entire trajectory under consideration. On almost the entire trajectory (except for its final section), the value of this projection is significant (of the order of 1 mm/s2). This ensures a decrease in the energy of the SC's selenocentric motion and a temporary capture of the SC by the Moon. The transversal component of the perturbing terrestrial acceleration (thin solid line) is also negative. This contributes to a decrease in the semilatus rectum of the selenocentric osculating orbit of the SC. The radial and normal components of the perturbing acceleration (dotted and dash-dotted lines) have less effect on the SC trajectory. Figure 5 shows the change in some osculating elements of the geocentric trajectory of the SC. The time interval from the launch of the SC from a low Earth orbit to the moment of approach of the SC to the libration point is considered. The left figure (a) shows the change in the osculating eccentricity. It can be seen how the solar gravitational acceleration reduces the eccentricity from the eccentricity of the intermediate orbit 0.989702 to 0.268581 at the moment of the maximum approach of the SC to the libration point. The central figure (b) shows the change in the semilatus rectum (dashed line) and perigee radius (solid line) of the osculating geocentric orbit. The effective increase in these elements due to the solar disturbance on a large initial part of the trajectory is replaced by some decrease as the SC approaches the libration point. The main reason for this is the lunar gravitational perturbation. The right figure (c) shows the change in the apogee radius of the osculating geocentric orbit. This element actively decreases as the SC approaches the libration point due to lunar disturbances. Figure 6 shows the change in the SC distance from the libration point D over the entire flight trajectory (a) and on the last six days of the flight (b). The minimum SC distance from the libration point (18 620 km) is reached on the 85.93th flight day. Fig 4.tiff a b Figure 4. Projections of the perturbing Earth's gravitational acceleration ΦEarth, mm/s2 (a) and projections of the perturbing solar gravitational acceleration ΦSun, mm/s2 (b) as a function of the flight time (day); the selenocentric section of the trajectory Fig 5.tiff a b c Figure 5. The osculating elements of the geocentric orbit of the SC as the functions of time (day) on the trajectory until the SC approaches to the libration point: a - the eccentricity; b - the semilatus rectum p and perigee radius rπ, thousand km; c - the radius of the apogee, million km Fig 6.tiff a b Figure 6. SC distance from the libration point, thousand km, as a function of time (day) on the entire flight trajectory (a) and on the last seven days of flight (b) At the moment of closest approach of the SC to the libration point, the radius vector of the SC is 422.8 thousand km. The geocentric radius vector of the libration point at this moment of time is 432.4 thousand km. That is, the SC is located closer to the Earth than the libration point by almost 10 thousand km. This explains why the radii of the apsidal points of the osculating geocentric orbit of the SC at this moment of time are less than the radii of the apsidal points of the osculating geocentric orbit of the libration point. For example, the perigee radius of the SC orbit is less than the perigee radius of the osculating orbit of the libration point by 30 400 km. Figure 7 shows the change in the energy constant hSel of the osculating selenocentric orbit as a function of flight time. It can be seen that on 81.371 day of flight this constant becomes negative and continues to decrease. There is a “capture” of the SC by the Moon. At the moment when the SC is at the minimum distance to the libration point (85.930 days of flight), the energy constant of the SC's selenocentric motion is -0.138 km2/s2. This point is shown as a black diamond on the graph. It is noteworthy that the energy constant of the selenocentric orbit changes significantly even on the last day of the flight. The right Figure 7 shows the change in this element of the osculating orbit during the last three hours of the flight. During these 3 hours, the value of the selenocentric radius of the SC vector decreases from 11 405.6 km to 1838 km. And gravitational perturbations from the Earth non-monotonically change the energy constant of the selenocentric motion. At the end of the flight (before the implementation of the braking impulse of velocity), the energy constant of the osculating selenocentric orbit is -0.176 km2/s2. At this moment, the elements of the osculating selenocentric orbit turned out to be as follows: eccentricity 0.934028, periapsis radius 1837.996 km, apoapsis radius 53882.842 km; the true anomaly of the osculating orbit is 359.836o. The projections of the selenocentric trajectory when the spacecraft approaches the LMO and the projections of this orbit are shown in Figure 8. The approaching trajectory of the SC (a highly elongated elliptical orbit) practically touches the LMO (shown by a thin line). Fig 7.tiff a b c Figure 7. The change in the energy constant of osculating selenocentric orbit, km2/s2, as a function of time (day) on the flight trajectory: a - during of the entire trajectory of the flight; b - during of the last six days of the flight; c - during of the last three hours of the flight Fig 8.tiff a b Figure 8. Projections onto the x-y plane of the ecliptic (a) and onto the z-y plane (b) of the selenocentric trajectory when the SC approaches the LMO; distance unit is 1 thousand km 2.2. Conclusions on the numerical analysis The analyzed low-energy trajectory of the flight to the Moon, as we assumed, pass through the vicinity of the L2 libration point of the Earth - Moon system. The minimum distance to this libration point turned out to be 18.6 thousand km. At the time of the minimum approach of the SC to the libration point, the SC was closer to the Earth than the libration point (by 9.60 thousand km). Apparently, this explains the fact that the radii of the perigee and apogee of the osculating orbit of the SC at the moment of minimal approach to the libration point seemed to be less than the radii of the perigee and apogee of the osculating geocentric orbit of the libration point. This did not prevent us from using the trajectory with the minimum total miss of the libration point I (1) as the initial approximation to obtain the trajectory of the temporary capture of the SC by the Moon. At the moment of passage of the libration point, the selenocentric energy constant is negative (-0.138 km2/s2). An analysis of the level lines in Figure 1 makes it possible to assume that the launch window is wide enough for the considered type of flight trajectories in the considered time interval of 2024. The duration of the launch window is about 30 days (March 20 - April 20). Within the launch window, the intermediate orbit apogee radius (ra) is a nonmonotonic function of the launch date. The developed method is quite labor-intensive, and we expect future improvement. The solution of the following problems is considered. When using local search methods (at the second and subsequent stages of the developed algorithm), derivatives of the motion conditions at the end point of the flight trajectory are used. They are currently calculated using the central difference algorithm. Under conditions of high sensitivity of the considered trajectory, it is very difficult to achieve high accuracy of these derivatives. A possible solution to the problem of the accuracy of these derivatives is to use the apparatus of complex numbers or dual numbers [13; 14]. Calculation of derivatives with high accuracy can make it possible to use the necessary optimality conditions for the constrained optimization problem of the flight trajectory. Since the optimality conditions themselves contain derivatives of the characteristics of the trajectory with respect to the choosing parameters of flight pattern, the use of local optimization methods will be impossible without finding the second derivatives. A new difficult problem arises, i.e., finding the second derivatives of the characteristics of the trajectory at the end of the flight with respect to the parameters of the flight pattern. Its solution can lie in the use of dual complex numbers [15]. Conclusion The numerical analysis showed the operability of the proposed method for designing trajectories for a low-energy flight to the Moon with the SC insertion into the low lunar orbit. The main feature of the method is the assumption that the flight trajectory lies in the vicinity of the libration point of the Earth - Moon system and during this passage some restrictions are introduced on the magnitude and direction of the SC velocity. We do not claim that all low-energy flight trajectories satisfy these conditions. The authors argue that there are low-energy transfer trajectories passing through the neighborhood of libration points and propose a method for finding such trajectories. On the flight trajectory, obtained using the developed method, it was possible to reduce the decelerating impulse of the velocity when entering a circular circumlunar orbit with a height of 100 km (in relation to the traditional flight scheme) by more than 140 m/s.## About the authors

### Mikhail S. Konstantinov

Moscow Aviation Institute (National Research University)
**Author for correspondence.**

Email: mkonst@bk.ru

ORCID iD: 0000-0002-0138-6190

SPIN-code: 3030-7494

Scopus Author ID: 55396771600

Doctor of Sciences (Techn.), Professor of the Space Systems and Rocket Science Department, Aerospace Institute

4, Volokolamskoye Shosse, Moscow, 125993, Russian Federation### Aung Myo Thant

Moscow Aviation Institute (National Research University)
Email: aungmyothant4696@gmail.com

ORCID iD: 0009-0000-1159-3292

PhD student, Space Systems and Rocket Science Department, Aerospace Institute

4, Volokolamskoye Shosse, Moscow, 125993, Russian Federation## References

- Parker JS, Anderson RL. Low-energy lunar trajectory design. Hoboken, New Jersey: John Wiley & Sons, Inc.; 2014. https://doi.org/10.1002/9781118855065
- McCarthy BP, Howell KC. Cislunar transfer design exploiting periodic and quasi-periodic orbital structures in the four-body problem. 71st International Astronautical Congress, The CyberSpace Edition, October 12-14, 2020. Paris; 2020.
- Scheuerle ST, McCarthy BP, Howell KC. Construction of ballistic lunar transfers leveraging dynamical systems techniques. AAS/AIAA Astrodynamics Specialist Conference, Lake Tahoe, California (Virtual), August 9-12, 2020. South Lake Tahoe, California; 2021.
- McCarthy BP, Howell KC. Trajectory design using quasi-periodic orbits in the multi-body problem. Proceedings of the 29th AAS/AIAA Space Flight Mechanics Meeting, 2019. Maui; 2019.
- Ivashkin VV. On the Earth-to-Moon trajectories with temporary capture of a particle by the moon. 54th International Astronautical Congress, Bremen, Germany, September 29 - October 3, 2003. Paper IAC-03-A.P.01. https://doi.org/10.2514/6.IAC-03-A.P.01
- Ivashkin VV. Low energy trajectories for the Moon-to-Earth space flight. Journal of Earth System Science. 2005;114:613-618. https://doi.org/10.1007/BF02715945
- Belbruno EA, Carrico JP. Calculation of weak stability boundary ballistic lunar transfer trajectories. Proceedings of the AIAAJ'AAS Astrodynamics Specialist Conference, August 14-17, 2000, Denver, Colorado. Paper AIAA 2000-4142. https://doi.org/10.2514/6.2000-4142
- Belbruno EA, Miller JK. Sun-perturbed Earthto-Moon transfers with ballistic capture. Journal of Guidance, Control, and Dynamics. 1993;16(4):770-774. https://doi.org/10.2514/3.21079
- Koon WS, Lo MW, Marsden JE, Ross SD. Low energy transfers to the Moon. Celestial Mechanics and Dynamical Astronomy. 2001;81(1):63-73. https://doi.org/10.1023/A:1013359120468
- Miller JK. Lunar transfer trajectory design and four body problem. 13th AAS/AIAA Space Flight Mechanics Meeting at Ponce, Puerto Rico, 2003. American Astronomical Society, American Institute of Aeronautics and Astronautics; 2003.
- Miller JK, Hintz GR. Weak stability boundary and trajectory design. Spaceflight Mechanics. Conference paper AAS 15-297. Williamsburgh, VA; 2015.
- Konstantinov MS, Thant AM. Designing the low-energy lunar transfers trajectories which pass in the vicinity of the libration points of the Earth - Moon system. Part 1. Theory and method. RUDN Journal of Engineering Re-search. 2023;24(1):7-16. http://doi.org/10.22363/2312-8143-2023-24-1-7-16
- Martins JRRA, Sturdza P, Alonso JJ. The complex-step derivative approximation. ACM Transactions on Mathematical Software. 2003;29(3):245-262. https://doi.org/10.1145/838250.838251
- Konstantinov MS, Nikolichev IA, Thein M. Optimization of low thrust multi-revolution orbital transfer using the method of dual numbers. Proceedings of the 6th International Conference on Astrodynamics. Tools and Technics (ICATT-2016). Available from: https://indico.esa.int/indico/event/111/session/21/contribution/99/material/paper/0.pdf (accessed: 12.11.2022).
- Petukhov VG, Yoon SW. Optimization of perturbed spacecraft trajectories using complex dual numbers. Part 1. Theory and method. Cosmic Research. 2021;59(5):401-413. https://doi.org/10.1134/S0010952521050099