Reachable set analysis of practical trans-lunar orbit via a retrograde semi-analytic model

Different from the Apollo flight mode, a safer trans-lunar flight mode for the crews is preferred. The previous-arrived lunar module rendezvouses with the crew exploration vehicle at a low lunar destination orbit, and then the crews ride the lunar module to descend the lunar surface sampling destination. The lunar module, which includes the descent and ascent stages, flies from a low Earth orbit to the low lunar destination orbit with two tangential impulses. The low lunar destination orbital reachable set of the practical trans-lunar orbit limits the feasible lunar surface sampling region. Therefore, this paper addresses the low lunar destination orbital reachable set of the practical trans-lunar orbit. A retrograde semi-analytic model is proposed for rapidly computing the practical two-impulse trans-lunar orbit firstly, which refers the ephemeris table twice for more precision perilune orbital elements. The reachable set is generated using the multiple-level traversal searching approach with this retrograde semi-analytic mode. Its envelope is re-checked by the continuation theory with the high-precision orbital model. Besides, some factors that affect the reachable set are also measured. The results show that neither the Earth–Moon distance nor the trans-lunar duration affects the reachable set. However, if the trans-lunar injection inclination is smaller than the inclination of the moon’s path, the reachable set becomes smaller or even reduces into an empty set. In brief, the proposed retrograde semi-analytic model for computing the reachable set provides a helpful and fast tool for selecting an applicable lunar surface sampling site for the manned lunar mission overall design.


Introduction
The moon as the nearest celestial body is also the compulsory training camp for human's interstellar travels. After the Apollo project, especially when the evidences of water in lunar south pole are discovered by many lunar polar orbiters (Colaprete et al. 2010), new plans for the manned lunar mission by USA and other countries have become more diverse. The National Aeronautics and Space Administration (NASA) have published three important reports on the Constellation Program after 2005 (Stanley et al. 2005;Condon 2007; Garn et al. 2008). Although the Constellation Program was replaced later by the Lunar Orbital Platform-Gateway Program (Davis 2018), the aims of lunar Global Access and Anytime Return remained unchanged. Different from the lunar low latitude area detection of the Apollo project, the lunar polar and high-latitude districts would be the next manned lunar mission destination for water and possible lunar life. In the Constellation Program, NASA proposed a new idea called crew and cargo separation that the lunar module (i.e., descend and ascend stages are merged) transferred to a low lunar orbit (LLO) by a practical trans-lunar duration about 4-5 day, to wait for a crew exploration vehicle for rendezvous and assembling. Compared with the Apollo flight mode, this new flight mode is safer for crews (Stanley et al. 2005). When it comes to the mode of a 3-day trans-lunar orbit used in the Apollo Project, the new mode is also more fuel-saving for translunar injection (TLI) and low lunar orbit insertion (LOI) (He and Shen 2020). To save fuel, TLI and LOI are hoped to be tangential to the periapsis velocity vector. With the two constraints (i.e., the trans-lunar duration and tangential velocity increments), the orbital elements reachable set of LLO for the lunar module rendezvous with the crew exploration vehicle becomes a significant factor to select a feasible lunar surface sampling site. Therefore, this paper studies the reachable set of the practical twoimpulse trans-lunar orbits.
In February 1959, the Luna-1 lunar impact probe launched by Soviet Union did not reach the goal of hitting the lunar surface because of its inaccurate translunar orbital dynamics model and designed method. The Luna-2, launched in September 1959, successfully hit the lunar surface and became the first human lunar probe. After that, dozens of lunar probes including the unmanned and manned cases have had successfully missions to the Moon. However, the well-known translunar orbital dynamics models could be divided into only four items, namely, the circular restricted threebody problem (CR3BP), the double two-body model, the pseudo-state model, and the high-precision model. Issac Newton described the CR3BP in his Philosophiae Naturalis Principia Mathematica, which shows the mathematical formulas of the orbits in a three-body system. CR3BP still need numerical integration (Lee et al. 2014). The double two-body model was proposed based on the principle of computing the celestial body's gravitational sphere in the Pierre Simon Laplace's monumental work Celestial Mechanics (Egorov 1958). Wilson (1970) and Byrnes and Hooper (1970) proposed the pseudo-state model around 1970s. It is also named the multi-conic method. It calculates a trans-lunar orbit using both the Earth's and the Moon's gravitational forces acting on the spacecraft gradually. Such calculation consumes about 1% of the numerical integration time, and its model errors are no more than 5% compared with the double two-body model. The perilune attitude error for a trans-lunar orbit is about 20 km (Zhou et al. 2019). Apart from the above trans-lunar orbital dynamical models, the high-precision model considers all the known perturbation forces acting on the spacecraft. Thus, it is the most accurate dynamic model that reflects the real scenario with a huge potential in the near future (Yan and Gong 2011). To design a trans-lunar orbit with a high-precision model, a simplified orbital model is applied usually to provide an initial value of the design variable (Mohammad and Roberto 2019). The double two-body model is selected usually to play the role of the simplified orbital model because of its semi-analytical feature. However, the selections of the parameters are different. Peng et al. (2011) and Gao et al. (2018) selected the entry point parameters on the lunar influence sphere to play the orbital patched conic parameters. Li et al. (2015) selected the exit point parameters to do this. The above orbital patched conic techniques that select the entry or exit point are difficult to compute the precision orbital position and velocity at the epochs of the perilune and perigee. Because they referred the Jet Propulsion Laboratory (JPL) ephemeris table for the Moon's position and velocity at the epochs of the entry or exit point, but ignore the Moon's position and velocity deviations when the spacecraft flies inside the Moon's gravitational influence sphere. And the lambert iteration is needed to compute the orbital elements in the Earth-centric segment, so its calculation efficiency is low. A set of the perilune parameters selected in our previous works ) is an effective measure for avoiding the orbital perilunestate deviation.
The reachable set of non-linear differential equations refers to the dynamic system as Eq. (1). If the system is consecutive within t ∈ [t 0 , t f ] , and there is an initial state set x(t 0 ) ∈ � n ⊆ R n at the moment t 0 . There exists a control set u(t) ∈ U m ⊆ R n that causes the final state x(t f ) ∈ � n ⊆ R n at the moment t f , then the n can be called the reachable set of the initial set n (Grantham. 1981).
Here, x(t) denotes the state variables of the consecutive system; ẋ(t) is its differential variable. f [] denotes their dynamic model function. y(t) denotes the observation variables or the variables which are concerned. c[] denotes their function. In a physical dynamic system, all of the elements in the final reachable set are sometimes (1) not intuitive or concerned. However, its partial elements are concerned. Therefore, a relevant elements' equation y(t) = c[t, x(t)] is built, which can be used for calculating the final concerned parameters' set Y k ⊆ R n . If there are some inequality or equality constraints in the control process as Eq. (2), then Y k will become small, dimensionless, even an empty set Y k → .
Here, g[] and ζ [] denote the inequality or equality constraints. To review the approaches of computing the orbital reachable set, the previous literatures mainly focused on the following aspects, such as the spacecraft formation flight (Lee and Hwang 2018), collision probability calculation (Bai et al. 2012), geostationary satellites collocation (Li 2014), engagement analysis of exo-atmospheric interceptor (Chai et al. 2014), Earth re-entry (Lu and Xue 2010), and Moon and Mars entry landing footprints (Arslantas 2016;Benito and Mease 2012). Methods vary in different cases, but all of them could be categorized essentially into two types.
Type I: When a spacecraft is dominated by one force and the orbital state transition matrix can be derived by the linearization hypothesis in some form as ẋ = Ax + Bu , and if there is no control force or uncertain process model error in the duration t ∈ [t 0 , t f ] , a state transition matrix Φ (t 0 ,t f ) = e A(t f −t 0 ) can be utilized to compute the final state x t f = Φ (t 0 ,t f ) x t 0 . Therefore, the final state reachable set, affected by an initial state uncertainty or a small orbital maneuver Wen et al.2018;Yang et al. 2019), can be approximately calculated by the state transition matrix (i.e., δx t f = Φ (t 0 ,t f ) δx t 0 ). If there are cases such as some perturbation forces action on the spacecraft (Li et al.2011), large uncertainty range of the initial state, long initial-final duration (He et al. 2013), or a complex flight process (Feng et al.2019), the calculation of the reachable set by the utilization of linearization hypothesis leads to an inaccuracy, or even a wrong result.
Type II: When a spacecraft has some orbital control powers or there are some uncertain processes, such as the general entry or re-entry reachable set problems with uncertain atmospheric model parameters, the final state becomes more complex as Because of the unknowable controlling or uncertain perturbation forces t f t 0 Bu · dt during the process, the reachable set cannot be computed by the state transition matrix approach. Komendera et al. (2012) had suggested an intelligent feedback-adjust idea to search the initial parameters of a reachable set's envelope, but it also needs more computations and is not practical for the complex engineering problems. The practical approach is the two-step method which avoids a huge amount of calculations. Take the Earth re-entry for example (Lu and Xue 2010). The first step computes the maximum range, and the second step computes the maximum cross-ranges with some fixed-ranges which are between the minimum and maximum ranges. The high-precision dynamical model is adopted instead of a linearization assumption. By this way, the reachable set calculated by the high-precision method is undoubtedly more accurate than that calculated by the linearization hypothesis. This method is a calculation strategy operated at two levels. The outer layer gradually gives the virtual objective point, while the inner layer is a general single-objective non-linear constrained optimization problem. According to the specific conditions, the inner optimization problem can be solved by an indirect method or a direct optimization algorithm. The basic model is shown in Eq. (3).
Here, n, m, l, p denote the number of the optimization variables, optimization objectives, inequality constraints, and equality constraints, respectively. f , g j , h k denote the objective function, inequality constraint function, and equality constraint function, respectively. The s.t. is the abbreviations for the subjected to some constraints. The optimization results of the inner layer calculation constitute the envelope parameters of the reachable set. As for the shortcoming of this methodology, it is impossible to prove whether the reachable set is continuous or not. The structure of the remainder of the article is as follows. After reviewing the above-mentioned trans-lunar orbital dynamical models and different approaches of the orbital reachable set generation, details of problem formulation including the generation strategy proposed in this paper are presented in "Problem statement and strategy" section. The retrograde semianalytical model and the multiple-layer traversal searching approach are described in "Reachable set rapid generation" section. In "Precision envelopes generation" section, the retrograde semi-analytical model and the envelope of the reachable set are re-checked by the high-precision dynamics model. The property of the reachable set changes due to some influence factors are exhibited in "reachable set property measurement" section. Finally, conclusions are drawn in "Conclusion" section.

Problem statement
The practical two-impulse trans-lunar orbit has duration about 4 ~ 5 days, projected orbital altitudes at the epochs of TLI and LOI, and tangential incremental velocities of v TLI and v LOI as shown in Fig. 1. Tangential impulses are reasonable for computing a trans-lunar orbit. It is uneconomical that the trans-lunar burn from a low earth orbit (LEO) with a non-tangential thruster vector. It is same uneconomical for lunar orbit insertion with a non-tangential thruster vector. The TLI burn duration of Apollo-11 is 320 s, it is an impulse compared to the trans-lunar duration about 4-5 days (Berry 1970). To describe the orbital reachable set of the practical two-impulse trans-lunar orbit via the mathematical formula in Eq. (1), t 0 and t f represent the epochs of TLI and LOI, respectively. The orbital initial state set is described as Eq. (4).
Here, κ represents the radius of periapsis; subscript 'EJ2' represents the J2000.0 earth-centric coordinate system; h LEO represents the altitude of LEO; R E represents the radius of Earth; e , i , and f represent the eccentricity, inclination, and the true anomaly of the trans-lunar orbit, respectively. The superscript of 'lb' and 'ub' represent the low and upper boundaries. The expression of e EJ2 < 1 (4) . Fig. 1 Illustration of the practical two-impulse trans-lunar orbit implies the trans-lunar orbit is an elliptical orbit, which orbital energy is lower than the parabolic and hyperbolic orbits. The interval of i lb EJ2 , i ub EJ2 depends on the launching angle A 0 and the latitude of launch of launch site of B 0 , i.e., cos The nominal orbit does not need any mid-course corrections. The reachable set of the orbital elements at t f can be described as Eq. (5). Here, the subscript of 'MJ2' represents the J2000.0 mooncentric coordinate system; h LLO represents the altitude of LLO; and R M represents the radius of Moon. The expression of e MJ2 > 1 implies that the trans-lunar orbit is a hyperbolic orbit when it arrives at the perilune point. The equation of f MJ2 = 0 is equivalent to that v LOI acts on the perilune point. The interval of t lb , t ub implies that the acceptable trans-lunar duration is practical and useful in an engineering task.
Obviously, the reachable set as described in Eq. (5) is a multi-dimensional element set, which is difficult to present clearly and simply understand. However, only the inclination and the longitude of ascending node (LAN) in the lunar centric-fixed coordinate system (i.e., the subscript of 'LCF') in Eq. (6) are the elements concerned, which affect the ground track of satellite (GTS) after lunar orbit insertion as shown in Fig. 2

Two-step strategy
It implies that there are countless trans-lunar orbits with constraints as Eq. (7) to be calculated for generating a valid Y k t f . The difficulty lies in two aspects, one is the orbital dynamics model and another is the optimization algorithm. If the high-precision orbital mode is applied, its computing-time for a large number of the orbits is unaffordable. And if the evolutionary algorithm is applied, it leads to the same problem of the cost of computing.
Reference Type II. mentioned in the bottom of the Introduction, a strategy is suggested, which contains two steps as shown in Fig. 3. The simple explanation is as follows: Step 1: A retrograde semi-analytical model is established based on the double two-body concept, wherein the orbital elements at the perilune point are selected to play the role of the traversal searching variables. After traversal searching, the orbital elements, which satisfy the constraints, are recorded. Then, the topological structure and the influence relations of the elements in the reachable set are analyzed.
Step 2: The orbital design variable, which affects the reachable set obviously, is selected to play the role of the numerical continuation element. Then, the precision envelops of the reachable set are re-checked with the high-precision orbital dynamical model under the continuation frame. Every point of the envelope implies a constrained optimal problem. Its solution is solved via the initial values from Step 1 via SQP optimization algorithm (Gill et al. 2005).
Step 1 not only provides the topological structure information and the influence relation of the design variables in the reachable set but also provides the initial value of the orbit design variables for Step 2.
Step 2 checks the (7) reachable set solved by Step 1 using the proven high-precision orbital dynamical model.

A retrograde semi-analytical model
The latitude and longitude on the Moon's influence sphere, and flight duration from TLI to the moment when the lunar probe passes into the Moon's influence sphere are selected as the trans-lunar orbital design Here, r prl is equals to κ MJ2 ; M y ,M z , and unused M x are the basic rotation matrix of the coordinate systems. If t f (i.e., the epoch of LOI) and r prl (i.e., r prl = R M + h LLO ) are set to constants, the instantaneous lunar-centered LVLH coordinate system can be treated as an inertial system, the position, and velocity vectors of the trans-lunar orbit at this epoch before LOI in J2000.0 moon-centered coordinate system can be computed by Eq. (9).
Here, M LVLH2MJ2 is the rotation matrix from the lunarcentered LVLH coordinate system to the J2000.0 lunarcentered coordinate system which is expressed of The expressions of (u M , i M , � M ) are the Moon's argument of latitude, inclination, and right ascension of ascending node (RAAN), respectively, in the J2000.0 lunar-centered coordinate system. Therefore, the trans-lunar orbit can be only computed by , ϕ, i prl , v prl . When r MJ2 prl , v MJ2 prl is translated into the modified classical orbital elements , the radius of the Laplace influence sphere ρ is 66200 km. Hence, the true anomaly f MJ2 in at the moment of t in can be computed by Eq. (10).
Here, the subscript of 'in' represents the epoch when the lunar module enters the Laplace influence sphere, and p prl = κ prl 1 + e prl . Then, the position and velocity vectors of the trans-lunar orbit r MJ2 in , v MJ2 in at this epoch in J2000.0 lunar-centered coordinate system can be computed by the principle of the two-body orbital state transfer matrix described as Eq. (11).
Here, µ M represents the lunar gravitational constant, and According to Gudermann Christoph's transformation principle, the lunar-centered hyperbolic anomaly H at this epoch can be computed by Eq. (13).
The flight duration from t in to t f can be computed by Eq. (14).
Earth-centered coordinate system can be computed by Eq. (15).
When r EJ2 in , v EJ2 in is translated into the modified classical orbital elements κ EJ2 , e EJ2 , i EJ2 , � EJ2 , ω EJ2 , f in EJ2 , the eccentricity of e EJ2 is smaller than 1, and f in EJ2 is the true anomaly of the trans-lunar orbit at the epoch of t in . According to the Kepler's transformation principle, the flight duration from the moment of TLI to the moment t in can be computed by Eq. (16).

Multiple-layer traversal search
Based on the retrograde semi-analytical model, a multilayer traversal search method is utilized to generate the reachable set of the trans-lunar orbits rapidly. The flow chart is shown in Fig. 5. , ϕ, i prl , v prl are selected as the traversal searching elements in a four-layer frame respectively. Then, κ EJ2 , e EJ2 , i EJ2 , and, f in EJ2 are tested by the constraints. If any constraint is not satisfied, the next loop will be re-started. Considering the error of the double two-body model, an error of the radius of perigee �κ EJ2 is allowable as Eq. (17).

Fig. 5 Flow chart of the multiple-layer traversal search method
Until the searching process of the four layers is accomplished, the set of the orbital elements at the epoch of perilune constitutes the reachable set.

The reachable set based on the retrograde semi-analytical model
To verify the effectiveness of the proposed retrograde semi-analytical model and the multi-layer traversal search frame, an example is tested. The useful orbital elements and constraints are optioned as follows: A) The moment of perilune t f is set as 1 Jan 2025 0:00:00.000 UTCG (UTC, in Gregorian Calendar format). B) Refer to Apollo-11 mission (Berry 1970 After a pure analytical searching process, a mass of the constrained trans-lunar orbits is obtained. The distributions of the design variables are shown in Fig. 6a, b,  The orbital elements concerned as shown in Fig. 6d are the inclination and the LAN in the lunar-centric LVLH and the lunar fixed coordinate systems. The inclination in the lunar-centric LVLH coordinate system is scattered in  deg. However, when the inclination is less than about 167 deg, the LAN only distributes between about [80-140] deg and  deg. The case in the lunar fixed coordinate system is similar, except for minor differences caused by the Moon's rotation and libration.

Optimize a precise trans-lunar orbit
The position and velocity vectors of the trans-lunar orbit at the moment of perilune before LOI in J2000.0 lunar-centered coordinate system r MJ2 prl , v MJ2 prl can be computed by Eq. (8)  The trans-lunar orbit is computed via a 6-day retrograde-time numerical integration with the initial states of r EJ2 prl , v EJ2 prl using the Runge Kutta 7-8 integrator. The epoch of the perigee is found via the interpolation. Then, the position and velocity are translated into the modified classical orbital elements κ EJ2 , e EJ2 , i EJ2 , � EJ2 , ω EJ2 , f TLI EJ2 . Wherein, sin f TLI EJ2 < ε , which is the determined by the interpolation accuracy. κ EJ2 and i EJ2 are set as the equality and the inequality constraints due to the rocket. Therefore, the inclination and the RAAN at the epoch of LOI are the same with the target values. The model of searching a precise trans-lunar orbit via an optimal algorithm is shown in Eq. (19).
Here, the superscript 'tar' represents the target values of the orbital elements. To verify the model above, the same constraints and the interval of the design variables are set as that in section (the reachable set based on the retrograde semi-analytical model). The difference is the consideration of the precision perturbation orbital model. It contains the Sun's and the Moon's perturbation, A set of proximate design variables and the lunarcentered modified classical orbital elements in different lunar-centered coordinate systems are listed in Table 1 and Table 2. The Earth-centered modified classical orbital elements at the moment of TLI computed by the retrograde semi-analytical model are listed in Table 3. Table 1 are selected to play the role of the initial values of the design variables in Eq. (19). The optimal values of the design variables, the lunar-centered modified classical orbital elements in different lunar-centered coordinate systems, and the Earth-centered modified classical orbital elements computed via Eq. (19) with the high-precision numerical integration are listed in Tables 4, 5, and 6, respectively.

The values in
The iterative process of the optimal objective function and the constraints' fitness are shown in Fig. 7. The iteration converges quickly within 10 steps. It has two major    reasons: One is the utilization of the SQP_snopt optimization kit (Gill et al. 2005), which has good ability to solve non-linear programming problems (Song et al. 2020).
Another is the design variables obtained from the retrograde semi-analytical model, which provides an effective initial values.

Re-check the envelope of the reachable set
The reachable set is generated rapidly and clearly based on the retrograde semi-analytical model, and the accuracy of its envelope still has errors. According to the methodology in Eq.
(3) and Fig. 6(c), i prl is selected to play the role of the continuation element. In the continuation frame, there are three design variables left. i min EJ2 , i max EJ2 is still set as an inequality constraint, and the trans-lunar orbit is numerical integrated via a timeretrograde manner with a fixed t . The error absolute value between the perigee radius and the radius of LEO is set as the optimal objective function. The frame is shown in Eq. (20).
The initial values x 0 of the design variables are found from the reachable set in section (the reachable set based on the retrograde semi-analytical model). x opt n−1 plays the role of the initial values for the iteration optimization that n ≥ 2 as shown in Eq. (21). The orbital (20)  . 7 The iterative process of the optimal objective function and the constraints fitness elements are solved quickly and easily using this continuation theory, the envelope of the reachable set are constituted with them.
Corresponding to the time-retrograde integrate fixed time t ∈ [3, 6] day, respectively, the envelope of the reachable set in the lunar-centered LVLH and the lunar fixed coordinate systems is shown in Fig. 8.
It shows that the model error of the retrograde semianalytical model for the practical trans-lunar orbit with 6-day duration is larger than that with 3-day flight duration. When the flight duration is 3 days, the (21) if n = 1 : envelope solution of the reachable set is more consistent with the solution computed with the high-precision model. The inclinations in the lunar-centered LVLH and the lunar-centered fixed coordinate systems both approach to 180 deg. When the inclination is close to 180 deg, the LAN has a large error than other cases due to its numerical singularity.

Reachable set property measurement
The property of the reachable set changes due to some influence factors. It is useful to understand the relation of them for designing a manned lunar task in the overall design phase. In consideration of the effectiveness of the retrograde semi-analytical model, the property of the reachable set is measured. The distance of the Earth-Moon, the transfer duration, and the declination of the Moon are tested as follows, respectively.

The distance of the Earth-Moon
The eccentricity of the Moon's path changes from 1/23 to 1/5. It leads a result that the distance of the Earth-Moon changes from about 3.6 × 10 8 to 4.1 × 10 8 m in a lunation. After 1 Jan 2025 0:000.000 UTCG, 8 Jan 2025 0:000.000 UTCG, and 21 Jan 2025 5:000.000 UTCG in this lunation are selected as the epochs for testing, which correspond the nearest and the farthest Earth-Moon distances, respectively. The results are shown in Fig. 9a, b. It shows outwardly that the Earth-Moon distance leads non-obvious influence to the reachable set except for the orbital inclination is approaching to 180 deg in the lunar-centered LVLH frame. In practice, however, when the orbital inclination is approaching to 180 deg in the lunar-centered LVLH frame, the singularities of LAN emerge. Arbitrary value of LAN does not have significance. In a word, the Earth-Moon distance leads non-obvious influence to the reachable set.

The transfer duration
The reachable sets with the trans-lunar duration of 3 days, 4 days, 5 days, and 6 days are shown in Fig. 10, respectively. All of them arrive at the epoch of 1 Jan 2025 0:00:00.000 UTCG as perilune.
It shows that the LANs in the lunar-centered LVLH and the lunar-centered fixed coordinate systems have distinguishing features. Both of them has two manners (i.e., ascent orbit and descent orbit) arriving their perilune. When the trans-lunar duration is longer than 4 days, the Fig. 9 a The inclination and LAN in LVLH frame and LCF frame with the nearest value of the Earth-Moon distance. b The inclination and LAN in LVLH frame and LCF frame with the farthest value of the Earth-Moon distance perilune point is on the front of the Moon (i.e., the side is directly visible from the Earth). This property is significant for directly tracking and controlling from the Earth's surface tracking stations.

The declination of the Moon
The Moon's position in the J2000 Earth-centric coordinate system is described via the right ascension and declination. The declination of the Moon describes the Moon's position along the perpendicular direction of the J2000.0 mean equatorial plane. After 1 Jan 2025 0:000.000 UTCG, 7 Jan 2025 00.00.00.000 UTCG and 12 Jan 2025 00.00.00.000 UTCG in this lunation are selected as the epochs, which correspond the smallest and the biggest cases of the moon's declinations. The values of the Moon's declinations are 1.359 deg and 28.447 deg, respectively. The reachable sets are shown, respectively, in Fig. 11a, b. Refer to the discussion of the singularities of LAN when the orbital inclination is approaching to 180 deg, when we compare the difference between Fig. 9a, b. The reachable sets have non-obvious difference if just looking at these two pictures. However, the distributions of the LEO inclination as shown in Fig. 12 have a significant difference. It implies that if the inclination of LEO is less than the declination of the Moon, no trans-lunar orbit exists.
In other words, the reachable set is an empty set. The long period of the declination of the Moon is the Metonic cycle (i.e., 18.6 years), which periodically changes from 18.3 to 28.6 deg. The maximum declination of the Moon will be 28.6 deg in 2025 and the minimum is 18.3 deg in 2034. As long as the inclination of LEO is greater than 18.3 deg in 2034, the reachable set of the trans-lunar orbit will not be an empty set. Moreover, the reachable set of the trans-lunar orbit does not become a non-empty set at any time in 2025, and the inclination of the LEO must be greater than 28.6 deg.

Conclusion
A retrograde semi-analytical model is suggested to generate the reachable set of the practical two-impulse translunar orbit with two tangential maneuvers. Compared with the traditional two-body conic method in the literatures, this model does not have internal iteration, and its perilune point altitude can be an accurate constant. The multi-layer traversal searching frame and the precision envelops re-check have exhibited its rapidity and effectiveness. The classical influence factor test exhibits the changeable property of the reachable set. Some conclusions can be drawn as follows: