Review of Lambert's problem

REVIEW OF LAMBERT’S PROBLEM David de la Torre Sangr` a(1) and Elena Fantino(2) (1)Polytechnic University of Catalonia (UPC), E.T.S.E.I.A.T. calle Colom 11, 08222 Terrassa (Spain), david.de.la.torre.sangr[email protected] (2)Space Studies Institute of Catalonia (IEEC), Polytechnic University of Catalonia (UPC), E.T.S.E.I.A.T., Colom 11, 08222 Terrassa (Spain), [email protected] Abstract: Lambert’s problem is the orbital boundary-value problem constrained by two points and elapsed time. It is one of the most extensively studied problems in celestial mechanics and astrodynamics, and, as such, it has always attracted the interest of mathematicians and engineers. Its solution lies at the base of algorithms for, e.g., orbit determination, orbit design (mission planning), space rendezvous and interception, space debris correlation, missile and spacecraft targeting. There is abundance of literature discussing various approaches developed over the years to solve Lambert’s problem. We have collected more than 70 papers and, of course, the issue is treated in most astrodynamics and celestial mechanics textbooks. From our analysis of the documents, we have been able to identify five or six main solution methods, each associated to a number of revisions and variations, and many, so to say, secondary research lines with little or no posterior development. We have ascertained plenty of literature with proposed solutions, in many cases supplemented by performance comparisons with other methods. We have reviewed and organized the existing bibliography on Lambert’s problem and we have performed a quantitative comparison among the existing methods for its solution. The analysis is based on the following issues: choice of the free parameter, number of iterations,generality of the mathematical formulation, limits of applicability (degeneracies, domain of the parameter, special cases and peculiarities), accuracy, and suitability to automatic execution. Eventually we have tested the performance of each code. The solvers that incorporate the best qualities are Bate’s algorithm via universal variables with Newton-Raphson and Izzo’s Householder algorithm. The former is the fastest, the latter exhibits the best ratio between speed, robustness and accuracy. Keywords: Lambert, Orbits, Two-Body Problem, Transfer Time Equation, Root Finding Algorithms. 1. Introduction This contribution deals with a prominent piece of Celestial Mechanics and Astrodynamics, the TwoBody Orbital Boundary-Value Problem (TBOBVP), also known as Gauss’ or Lambert’s Problem. It is the problem of determining the Keplerian orbit connecting two positions in a given time. Figure 1 illustrates its definition: F is the primary (the Sun in the case of interplanetary orbits, the Earth in the context of geocentric motion), origin of the reference frame, while r1 and r2 are the positions occupied by the secondary at times t1 and t2 , respectively (with t2>t1 ), thus being ∆t=t2−t1 the flight time. The angle θ between r1 and r2 indicates the direction of motion. Finding the conic section that connects the two positions in the given time is equivalent to determining the velocity v1 to be imparted to the secondary at r1 to execute the transfer. The case shown as an example in Fig. 1 is an arc of an ellipse with one focus at F . The TBOBVP problem arose towards the end of the 18 th century in connection with the determination of the orbits of celestial bodies from observation. Carl Friedrich Gauss (1777-1855) used three observation times instead of two (and 1 Figure 1. Definition of Lambert’s problem. with this method he was able to correctly determine the orbit of the newly-discovered Ceres). Gauss understood the potential of the TBOBVP and his developments on this subject were remarkable, with the result that the problem now bears his name. The work of Gauss merges with the previous findings of Johann Heinrich Lambert (1728-1777) who deduced the homonymous theorem: Given the gravitational parameter µ=GM , the time ∆t required to accomplish a given transfer is a function of the semimajor axis a of the orbit, the sum r1+r2 of the distances from the primary at the beginning and at the end of the transfer and the length c of the chord that connects such positions, i.e.: √µ∆t=f(a,r1+r2,c).(1) When the geometry of the radius vectors is fixed, there is only one free parameter left that wholly defines the transfer time between r1 and r2 . In the original formulation, such parameter is the unknown semimajor axis a. In general, once a suitable parametrization of the orbit passing trough the given points is obtained, it is possible to rephrase Lambert’s problem in terms of the evaluation of the parameter’s value, such that the corresponding orbit is characterized by a transfer time that matches exactly the prescribed one. There exist several ways to express f , but no analytical closed-form solution of the transfer time equation (Eq. 1) is possible. Especially valuable and useful are the unified forms, i.e., providing one formulation valid for all three types of conic sections. In modern Astrodynamics, Lambert’s problem has direct application in the solution of intercept and rendezvous, ballistic missiles targeting and interplanetary trajectory design. The relevance of this issue is confirmed by the vast related literature. Several authors since the time of Gauss have studied alternative formulations of the transfer time equation, adopting a variety of independent variables (free parameter) and solution methods. Due to the transcendental character of the transfer time equation, all the available approaches are based on a numerical procedure in which the value of the free parameter is searched iteratively. Then, the orbital elements or the state vector at departure and arrival are obtained by means of orbital mechanics relations. In this work, we present a review of the algorithms for the solution of Lambert’s problem, based on quantitative comparisons. We start by discussing the relationships among the many methods, and we indicate the main representatives of the several research lines, as we identified them from our study of the literature (Sect. 2.). Then, in Sect. 3. we illustrate the tests that we carried out in order to compare the selected methods in terms of performance. We discuss the results and we draw 2 conclusions in Sect. 4. 2. Solution methods: comparison and discussion The first formulations of Lambert’s problem date back to Lagrange [ 1 ] and Gauss [ 2 ]. However, studies on the subject turned intensive in the mid 1960s. As of today, more than 60 authors have proposed formulations and solutions to the problem. All these methods can be grouped into a number of major lines of research on the basis of the free parameter adopted. Here we highlight the most productive ones, and for each we indicate the progenitor algorithm and the successive improvements upon it: •Universal variables: –Lancaster & Blanchard [3], Gooding [4, 5], Izzo [6]. –Bate [7], Vallado [8], Luo [9], Thomson [10], Arora [11]. –Battin-Vaughan [12], Loechler [13], Shen [14], MacLellan [15]. •Semi-major axis: –Lagrange, Thorne [16, 17], Prussing [18], Chen [19], Wailliez [20]. •Semi-latus rectum (p-iteration): –Herrick-Liu [21], Boltz [22]. –Bate [7]. •Eccentricity vector: –Avanzini [23], He [24], Zhang [25, 26], Wen [27]. •Kustaanheimo-Stiefel (K-S) regularized coordinates: –Sim´ o [28]. –Kriz [29], Jezewsky [30]. The work of Sun [ 31 ] is also worth mentioning, being the first to extensively investigate and analyse the multi-revolution Lambert’s problem and derive an expression for the minimum transfer time in this case. For brevity, in what follows we shall indicate [ 3 ] with LB69, [ 5 ] with G90, [ 6 ] with I15, [ 7 ] with B71, [8] with V97, [12] with BV84 and [28] with S73. We have selected representatives of each line of work and we have studied them. The methods due to Lagrange and Gauss have been retained for historical reasons. From the Universal variables line, the selection includes B71, BV84, and the work by G90 and I15, both superseding the original algorithm of LB69. Finally, the work of S73 has been selected as representative of the K-S branch. In the remainder of this section, we will describe these algorithms. For reasons of space, the reader is referred to the original publications for details and more quantitative aspects. 2.1. Generic Lambert procedure In general, the procedure to be adopted to solve Lambert’s problem consists in: 1. computing the geometric parameters of the transfer; 2. obtaining an initial guess for the free parameter; 3. iterating on the transfer time equation until convergence; 3 4. computing the velocity vectors. 2.2. Transfer time equation The algorithm by Lagrange expresses the transfer time as a function of the semi-major axis a . This choice is inconvenient for two reasons: the solution contains a pair of conjugate orbits (i.e., the solution is not unique), and the derivative of the transfer time is singular when a corresponds to the minimum-energy orbit ( am ). To overcome this difficulty, [ 32 ] expresses the transfer time equation as a function of the universal variable x ( x2=1−am/a ) which is valid for all types of conic sections, thus obtaining a single-valued and monotonic function for the transfer time. In his Theoria Motus, Gauss develops a method to solve Lambert’s problem with a system of two equations and using the sector-to-triangle area ratio y . Then, B71 expands the original Gauss definition to any type of conic section using a power series expansion developed by [33]. In his book, B71 starts from the Gauss f and g functions and expresses the transfer time in terms of the universal variable z using Stumpff functions. z is defined as ∆E2 for elliptical orbits and as −∆F2 for hyperbolic orbits, being ∆E and ∆F the eccentric anomaly in the ellipse and in the hyperbola, respectively. Note that Bate’s time equation is singular at z= (2nπ)2 with n=1,2,3... , and [2nπ]2 and [2(n+1)π]2 are the limits of the domain of n -revolution elliptical orbits. Note also that there is a lower limit for z when the transfer angle is less than π . Such limit corresponds to a hyperbolic orbit with transfer time ∆t=0 . Lower values of z imply negative values for the transfer time. S73 derives the transfer time in regularized space by means of the Levi-Civita coordinate transformation. His free parameter z has a universal definition, valid for the three types of conic section. The transfer time is singular at z=nπ2 with n=1,2,3,... , and [2nπ]2 and [2(n+1)π]2 are the limits of the domain of n -revolution elliptical orbits. Also in this case, there exists a lower limit for the value of z for transfers with θ<π . The work also provides some insight on how to avoid and correct certain scenarios where numerical precision gets degraded. BV84 improve the original Gauss formulation in their Elegant Algorithm, extending the convergence region and treating θ=π transfers. The new system of equations is also expressed in terms of the sector area-to-triangle area ratio y and makes use of a continued fraction. The method, however, is singular for 2πtransfers. G90 follows the work of LB69 in terms of the universal variable x ( x2=1−am/a ), implementing a Halley cubic iterator, initial guess value formulas to ensure fast convergence, a power-series expansion in the near-parabolic range to avoid round-off errors, and expanding its multi-revolution capabilities. In his paper, G90 provides a Fortran code with three routines: one computes the time of flight (and several of its derivatives) as a function of x and using the formulation by LB69. The second routine finds x , while implementing all the initial guesses for all orbits and multi-revolution cases. The main code computes the geometric parameters of the transfer, calls the second routine to compute xfor each multi-revolution case, and computes the velocity vectors. 4 I15 also starts from the work by LB69, implements a new initial guess based on empirical results and a House-Holder iteration scheme that converges in about two or three iterations. The algorithm has been implemented in a python code available as part of the ESA PyKEP toolbox. For values of x close to one (i.e., near-parabolic range), I15 makes use of the Lagrange and Battin [ 32 ] equations in terms of the universal variable z . As part of the same toolbox, I15 also implements another version of the algorithm using only the Lagrange equation, a set of transformed variables and a Regula-Falsi iteration scheme. 2.3. Solution algorithm Several methods can be used to numerically solve the equations resulting from the several treatments described in Sect. 2.2.. Each method has its benefits and shortcomings. For example, the NewtonRaphson method needs a single starter and has a high convergence rate but works well only with single-valued monotonic functions. Singularities or steep changes in the derivative slow down the convergence, make the method converge to the wrong solution, or even diverge. The bisection method, requires two starters bracketing the solution and it is usually slower than Newton-Raphson, but it will always converge to the the correct solution when properly bracketed. The iteration methods implemented for the selected Lambert algorithms are presented below. The Newton-Raphson scheme is a Householder root-finding method of 1 st order. The method requires the derivative of the transfer time with respect to the free parameter, which we denote by w . At the ith iteration, wis updated through: wi+1=wi−t−TOF dt dw .(2) TOF is the target value for the time of flight, t is current value of the transfer time and dt/dw is its derivative with respect to w. The Newton-Raphson method is used by Lagrange and B71. The bisection method does not require any derivatives, since it evaluates the transfer time equation at the arithmetic mean point (e.g., ¯w ) of a solution-bracketing interval (namely [wlow,wup] ). Then, the resulting value for the time of flight t is compared against the target TOF , and the interval bounds are updated as follows: wlow =¯w t ≤TOF,(3) wup =¯w t >TOF. The bisection method is implemented in the B71 and S73 algorithms. The Gauss algorithm uses a successive substitution technique called the Gauss Method. The scheme recomputes all parameters with the current value of the free variable and updates it using the recomputed value of the parameters. BV84 also develops a successive substitution solving method, exploiting the property that the largest positive real root of the cubic equation is always the correct solution, and using the resulting equation to directly compute the solution for the cubic. G90 uses the Halley method, which is a Householder root-finding method of 2nd order. The first and second derivatives of the transfer time equation with respect to the free parameter are then required 5 and used as follows: wn+1=wn+Tdt dt2+Td2t 2 (4) where T=TOF −t , and dt , d2t are the first and second order derivatives of the time equation, respectively. I15 uses a 3rd-degree Householder method to iterate on the transfer time equation: wn+1=wnTdt2−Td2t/2 dt(dt2−Td2t)+d3td2t/6(5) where T=TOF −t , dt and higher-order derivatives are computed using LB69’s version of the transfer time equation. In the alternative version of the algorithm at the PyKEP toolbox, I15 uses a Regula-Falsi scheme as iteration method. The Regula-Falsi method requires a solution-bracketing interval as a starter and updates the iteration variable as follows: wn+1=w1T2−T1w2 T2−T1 (6) where w1 , w2 are interval limits, while T1 and T2 are the values for the difference tn−TOF at each interval limit. Once the new value of the iterator ( wn+1 ) is computed, the corresponding value of the time of flight ( tn+1 ) is obtained via the transfer time equation. Then, the interval is updated as follows: w1=w2,w2=wn+1,T1=T2and T2=Tn+1. The process is repeated until convergence. 2.4. Choice of the initial guess The initial guess, or starter, is introduced as input to the solution algorithm. Note that each solution algorithm may require a specific number of initial guesses (e.g., Householder-type methods require a single starter, bisection requires two starters, etc.). The starters proposed by the authors of the selected algorithms are the following: Lagrange starts from x=0, corresponding to a minimum-energy orbit. Gauss starts from y=1, which corresponds to a 1:1 triangle-to-sector area ratio. B71 starts from z=0 for the Newton-Raphson scheme, corresponding to a parabolic orbit. For the bisection method, [ 8 ] suggests an upper bound of zup =4π2 , corresponding to the t=2π asymptote, and a lower limit of zlow =−4π , which he claims is valid for most orbits except for highly-eccentric ones. In such cases, the lower limit should be extended. S73 proves that his time equation has a single solution in the interval z∈[zf,π2] , where π2 is the t=2π asymptote and zf is the lower limit for the z parameter. The value of zf corresponds to a specific value for θ<π transfers and zf=−∞ for θ>π transfers. S73 recommends to 6 use z= (θ/2)2 as a starter for the elliptical orbits domain ( z∈[0,π2] ) and z=0 as a starter for hyperbolic orbits (z∈[zf,0]). BV84 provides a specific initial guess x for his successive substitution method if the normalized time of flight corresponds to an elliptical orbit, and x=0 otherwise. G90 derives a set of semi-empirical initial starters for both single-revolution and multi-revolution transfers. For the single-revolution case, he uses a bilinear approximation method to analytically obtain an approximate solution for the x<0 region of Lancaster’s transfer time equation, and a weighted combination of approximate solutions for the x>0 region. The multi-revolution starters are far more complex and, for reasons of space, the reader is referred to the corresponding publications. I15 performs a piece-wise linear approximation of the transfer time equation. Inverting the previous relation provides an approximate solution which is then used as a starter. For the multi-revolution cases, I15 approximates the lower and upper asymptotes of the time equation for the corresponding n -revolution scenario, then he inverts the approximated curves to obtain the appropriate starters. For his Regula-Falsi implementation, I15 uses an empirical set of constant starters for the singlerevolution case or multi-revolution case. 2.5. Computation of the velocity vectors Once the free parameter has been approximated, the velocity vectors at the boundaries of the transfer must be determined. The methods employed by the several authors are highlighted below. 2.5.1. Orbital elements This method consists in expressing the semi-latus rectum p , the semi-major axis a , the eccentricity e and the true anomaly ν of the selected point as functions of the free parameter. Then, the velocity vpwq at either r1and r2in the perifocal reference frame is determined as: vpqw =rµ p[−sinν,e+cosν,0].(7) Then, with the inclination i , the argument of the periapsis ω and the longitude of the ascending node Ω, the pqw vectors can be rotated into the xyz frame: vi jk =Rz(Ω)Rx(i)Rz(ω)vpqw,(8) where Rk ’s are the elementary rotation matrices following the right-hand sign convention. The orbital elements method is used by S73 and BV84. 2.5.2. Gauss fand gfunctions The values for the Gauss functions f , g and ˙g are expressed as functions of the free parameter. Then the velocity vectors v1 and v2 at the beginning and at the end of the transfer, respectively, are found 7 through: v1= (r2−fr1)/g(9) v2= ( ˙gr2−r1)/g(10) Note that there exists a singularity in the computation of the velocity when g=0 , which corresponds to θ=πtransfers. This method is used by Gauss and B71. 2.5.3. Radial and transversal components The radial vr and transversal vt components of the velocity are expressed as functions of the free parameter. Then, either velocity vector vi(i.e., v1or v2) is calculated as: vi=vri ri ri +vti h×ri |h×ri|.(11) Here, ri/ri is the unit vector in the radial direction at Pi , and h×ri/|h×ri| is the unit vector in the transversal direction. Lagrange, G90 and I15 use this method. 3. Performance tests The methods discussed in Sect. 2. have been implemented following the description of their authors. Only single-revolution cases have been considered for this contribution. The implementation has been done in Matlab. This programming framework allows for a fast prototyping and includes powerful code analysis and debugging tools, which favours the algorithm-oriented analysis that we aim at in this work. Although a preliminary performance-oriented analysis has also been carried in order to obtain a relative comparison between methods, the correct practice in this case should be to program the algorithms in a high-performance compiled code (e.g., FORTRAN or C). The latter approach will be addressed in a future development of this work. 3.1. Application region Several simulations have been carried out in order to test the algorithms in a wide range of scenarios. The goal is to stress each method and detect its range of applicability, slow-convergence zones, singularities, etc. The scenario sets a transfer between two radial distances, r1 and r2 , such that r2/r1=2 . Then the transfer angle θ is varied between 0 and 2π at constant steps, and for each value of the transfer angle, the transfer time is also varied between 0 and 2π . The normalization has been applied which consists in that µ and r1 are unitary and the circular orbit with radius r1 has unit mean motion. This setup ensures that the algorithms are tried also in the tough cases, including θ=0,π,2π , parabolic orbits, highly-eccentric hyperbolic orbits. The results of the simulations are discussed below. Battin’s implementation of Lagrange’s equation fails for highly-eccentric hyperbolic orbits with TOF <π/4 (Fig. 5), corresponding to values of the universal variable x1 . In this scenario, the Gauss hypergeometric function 2F1 does not converge because its convergence radius (unitary) is exceeded. As a result, the equation suffers a steep change in the derivative with which the 8 Newton-Raphson method can no longer cope (Fig. 2). The method performs 10-20 iterations in the elliptical orbits domain (best case near θ≈π ), and executes 40-60 iterations when approaching the non-convergence zone of 2F1. x -101234 t -2 0 2Equation t Equation dt/dx Target TOF Iteration t Iteration dt/dx Figure 2. Lagrange Newton-Raphson solver in the low-TOF scenario. The Gauss Method algorithm expanded by Bate only works for scenarios where y≈1 , which includes both elliptical orbits with a small transfer angle and hyperbolic orbits with a moderate transfer angle. In general, the convergence region is given by the approximate relation θ<π−TOF , as shown in Fig. 5. The method, however, fails completely in other cases. It is worth mentioning that for θ>π , the value of the triangle-to-sector area ratio y becomes negative, which corresponds to a physically impossible situation. Also, there is a small region for high values of θ where the algorithm converges but the results are obviously incorrect. The precision of the results in the convergence region is also relatively low compared with other algorithms (more than 10−2 relative error on the magnitude of the velocity vector). Due to the adoption of the f and g functions for the computation of the velocity vectors, the method also intrinsically fails for the θ=π transfer. The Gauss Method by Bate makes approximately 30 iterations in the elliptical orbits zone and 10-20 in the hyperbolic orbits domain, but the iteration number quickly increases to more than 100 when moving away from the convergence region. B71 method using universal variables and Newton-Raphson with a constant starter z=0 (parabolic orbit) has two conflictive regions. One corresponds to large values of TOF and high θ in the elliptical orbits range (Fig. 5), where the solution is a point close to the z=4π2 asymptote and the derivative is steep. In this situation, Newton-Raphson has a high chance to “jump” over the asymptote and fall into the multi-revolution region, where the error is propagated further (Fig. 3). The second conflictive zone is associated to θ<π (type I, short-way transfers) and TOF <π/8 in the hyperbolic orbits region, where the solution falls close to the minimum negative value for z . In this case, Newton-Raphson has a high chance of falling into the imaginary region of B71 equation (Fig. 4), thus yielding imaginary numbers. These two failure regions can be avoided by improving the initial guess, for example with a preliminary search to bracket the solution within bounds and then selecting the geometric mean between bounds as a starter. In this case, the convergence is greatly improved, except for a few cases with extremely low or extremely high values of TOF . The number of iterations in any case ranges from 5-10. B71 method with bisection and the starter proposed by [ 8 ] works for all orbits except those with θ>3/4π and TOF <π/2 in the hyperbolic orbits region, where the initial bounds no longer bracket the solution (see Fig. 5). In such cases, he suggests to decrease the value of the lower 9 25th International Symposium on Space Flight Dynamics ISSFD [1] Lagrange, Académie des Sciences, 1788. [2] Gauss, Brown and Co., 1857. [3] Herrick & Liu, Aeronutronic, 1959. [4] Godal, T. Astronautik, 1961. [5] Sconzo, Astron. J., 1962. [6] Escobal, TRW Space Technol., 1964. [7] Lancaster, J. Spacecr. Rockets, 1966. [8] Pitkin, J. Astronaut. Sci., 1968. [9] Price, Comput. J., 1968. [10] Lancaster, Greenbelt, 1969. [11] Lancaster, Greenbelt, 1970. [12] Bate, Dover Publ., 1971. [13] Monuki, 1973. [14] Simó, Collect. Math., 1973. [15] Jezewski, Celest. Mech., 1976. [16] Kriz, Celest. Mech., 1976. [17] Vinh, Univ. of Michigan, 1977. [18] Prussing, J. Guid. Control. Dyn., 1979. [19] Sun, Astr. Spec. Conf., 1979. [20] Battin & Vaughan, Astr. Spec. Conf., 1983. [21] Sun, Acta Astronaut., 1983. [22] Sun, R. Report AM-035, 1983. [23] Battin & Vaughan, J. Guid. Control. Dyn., 1984. [24] Boltz, J. Astronaut. Sci., 1984. [25] Moulton, Dover Publ., 1984. [26] Sun, R. Report AM-043, 1985. [27] Taff, Celest. Mech., 1985. [28] Sun, Astrodynamics, 1986. [29] Gooding, Farnborough, 1988. [30] Loechler, MIT, 1988. [31] Sarnecki, Acta Astronaut., 1988. [32] Gooding, Celest. Mech. Dyn. Astron., 1990. [33] Klumpp, Astrodyn. Conf., 1990. [34] Klumpp, JPL Tech. Rep., 1991. [35] Peterson, Adv. Astronaut. Sci., 1991. [36] Nelson, J. Guid. Control. Dyn., 1992. [37] Thorne, J. Astronaut. Sci., 1995. [38] Vallado, McGraw-Hill, 1997. [39] Reynolds, 1998. [40] Battin, AIAA, 1999. [41] Prussing, J. Astronaut. Sci., 2000. [42] Han, Chin. Sp. Sci. Tech., 2004. [43] Shen, Adv. Astronaut. Sci., 2004. [44] Maclellan, 2005. [45] Izzo, J. Guid. Control. Dyn., 2006. [46] Avanzini, J. Guid. Control. Dyn., 2008. [47] Arora, Sp. Fl. Mech. Meeting, 2010. [48] Avendaño, Celest. Mech. Dyn. Astron., 2010. [49] Bando, J. Guid. Control. Dyn., 2010. [50] He, J. Guid. Control. Dyn., 2010. [51] Leeghim, J. Guid. Control. Dyn., 2010. [52] Zhang, J. Guid. Control. Dyn., 2010. [53] Arlulkar, J. Guid. Control. Dyn., 2011. [54] Der, Adv. Maui Opt. and Sp. Surv. Tech. Conf., 2011. [55] Luo, Celest. Mech. Dyn. Astron., 2011. [56] Thompson, J. Guid. Control. Dyn., 2011. [57] Zhang, Celest. Mech. Dyn. Astron., 2011. [58] Aretskin, CalTech, 2012. [59] Parrish, CalTech, 2012. [60] Wagner, 2012. [61] Ahn, J. Guid. Control. Dyn.,2013. [62] Arora, 2013. [63] Arora, Astr. Sp. Conf., 2013. [64] Chen, J. Guid. Control. Dyn., 2013. [65] Usachov, Sol. Syst. Res., 2013. [66] Arlulkar, J. Guid. Control. Dyn., 2014. [67] Thorne, J. Guid. Control. Dyn., 2014. [68] Wailliez, Adv. Sp. Res., 2014. [69] Wen, Acta Astronaut., 2014. [70] Wen, J. Guid. Control. Dyn., 2014. [71] Ahn, J. Guid. Control. Dyn., 2015. [72] Izzo, Celest. Mech. Dyn. Astron., 2015. CPU time [log 10 µs] 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 Izzo 2015 Izzo 2010 Gooding 1990 Battin 1984 Simo 1973 Bate 1971 NR Bate 1971 BI Gauss 1971 Lagrange 1977 197.54 326.01 128.49 82.32 170.47 235.86 196.20 113.81 103.72 Lambert performance | Solver CPU time [µs] | Samples: 1E5 Microsoft Windows 10 Pro | Matlab R2015a | IntelCorei7 CPU 860 @ 2.80GHz | 1024 KB cache Conclusions: Fastest algorithm: Bate12 algorithm via universal variables with Newton-Raphson iterator and velocity computed via Gauss f,g functions. Best overall ratio between speed, robustness and accuracy: Izzo72 algorithm via universal variables with 3rd order Householder iterator and velocity computed via radial and transversal components. Code implementations are highly succeptible to performance losses due to overhead and non-optimised coding. Izzo 2015 (HH) | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Izzo 2015 (HH) | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 0 1 2 3 x -10123 log10 ∆t -1 0 1 2Izzo 2015 | Transfer time equation θ = 000 deg θ = 180 deg θ = 360 deg θ = 540 deg θ = 720 deg TOF = 0.1π TOF = 2.0π TOF = 4.0π TOF = 6.0π TOF = 8.0π ∆t=          η(x)3Q(x)+4λη(x)/2+nπ/ρ(x)3/2Battin a3/2[(α−sin α)−(β−sin β)+2nπ]/2 Lagrange, elliptic (−a)3/2[(sinh α−α)−(sinh β−β)] /2 Lagrange, hyperbolic (x−λz(x)−d(x)/y(x)) /E(x) Lancaster & Blanchard 1 Izzo72 2015 Gooding 1990 | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Gooding 1990 | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 0 1 2 3 x -10123 log10 ∆t -1 0 1 2 Gooding 1990 | Transfer time equation θ = 000 deg θ = 180 deg θ = 360 deg θ = 540 deg θ = 720 deg TOF = 0.1π TOF = 2.0π TOF = 4.0π TOF = 6.0π TOF = 8.0π ∆t=        4/3(1 −q3) Exact solution, parabolic σx(f(x)) −q3σx(q2f(x)) Series, near-parabolic 2 E(x)x−λz(x)− d(x) y(x)Lancaster & Blanchard 1 Gooding32 1990 Battin 1984 | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Battin 1984 | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 5 10 15 x 0 0.2 0.4 0.6 0.8 1 y 1 1.2 1.4 Battin 1984 | Sector to Triangle area ratio θ = 000 deg θ = 045 deg θ = 090 deg θ = 135 deg θ = 180 deg TOF = 0.1π TOF = 1.0π TOF = 2.0π TOF = 3.0π x= 1−l 2  2 + m y2 − 1+l y2 0=y 3 −y 2 −h 1 (x)y 2 −h 2 (x) 1 Battin23 1984 Simo 1973 | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Simo 1973 | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 30 35 40 45 z [π2] -4 -3 -2 -1 0 1 2 3 4 log10 ∆t 0 2 4 Simó 1973 | Transfer time equation θ = 000 deg θ = 090 deg θ = 180 deg θ = 270 deg θ = 360 deg TOF = 0.1π TOF = 1.0π TOF = 2.0π TOF = 3.0π ∆t=2Pc 3(4z)+Q[c1(z)c2(4z)−2c0(z)c3(4z)] c3 1(z)2P−Qc0(z) µ 1 Simó14 1973 Bate 1971 (NR) | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Bate 1971 (NR) | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 2 4 6 8 10 z [ π2 ] -8 -4 0 4 8 12 16 20 24 28 32 36 log10 ∆t 0 2 4 Bate 1971 | Transfer time equation θ = 000 deg θ = 090 deg θ = 180 deg θ = 270 deg θ = 360 deg TOF = 0.1π TOF = 1.0π TOF = 2.0π TOF = 3.0π √µ∆t=x(z)3S(z)+A  y(z) 1 Bate12 1971 (Newton-Raphson) Bate 1971 (BI) | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Bate 1971 (BI) | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 35 40 45 z [π2] -8 -4 0 4 8 12 16 20 24 28 32 36 log10 ∆t 0 2 4 Bate 1971 | Transfer time equation θ = 000 deg θ = 090 deg θ = 180 deg θ = 270 deg θ = 360 deg TOF = 0.1π TOF = 1.0π TOF = 2.0π TOF = 3.0π √µ∆t=x(z)3S(z)+A  y(z) 1 Bate12 1971 (Bisection) Gauss 1971 | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Gauss 1971 | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 0 50 100 x 0 0.2 0.4 0.6 0.8 1 y -20 -10 0 10 20 Gauss 1971 | Sector to Triangle area ratio θ = 000 deg θ = 060 deg θ = 120 deg θ = 180 deg θ = 240 deg θ = 300 deg θ = 360 deg x=w/y2 −s y=1+X(x)[s+x] 1 Gauss12 1971 Lagrange 1977 | Status TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 OK BAD IN NO SOL MAX IT π on F,G ∞ NaN CMPLX Lagrange 1977 | Iterations TOF [π] 0 0.5 1 1.5 2 θ [deg] 0 90 180 270 360 10 20 30 40 50 x 01234 log10 ∆t -1 0 1 2 Lagrange 1977 | Transfer time equation θ = 000 deg θ = 090 deg θ = 180 deg θ = 270 deg θ = 360 deg TOF = 0.1π TOF = 1.0π TOF = 2.0π TOF = 3.0π  µ a3 m ∆t=4 3  2F1(x)+λ3 2F1(y)  1 Lagrange40 1977 F P1 P2 c θ r1 r2 Lambert Theorem Given the gravitational parameter µ=GM, the time ∆t required to accomplish a given transfer (angle θ) is a function of the semimajor axis a of the orbit, the sum r1+r2 of the distances from the primary at the beginning and at the end of the transfer and the length c of the chord that connects such positions. Vallado38 Bate12 Luo55 Thomson56 Arora62 Battin23 Loechler30 Shen43 MacLellan44 Lancaster10 Gooding29 Izzo72 Boltz24 Liu3 Bate12 Lagrange1 Thorne37,67 Prussing41 Chen64 Wailliez68 Avanzini16 He50 Zhang52 Wen70 Kriz16 Jezewsky15 Simó14 Universal Variables Semi-latus Rectum (p) Semi-major Axis (a) Eccentricity Vector (e) K-S Regularized Space Lagrange1 Gauss2 Herrick3 Godal4 Sconzo5 Escobal6 Lancaster7 Pitkin8 Price9 Lancaster10 Lancaster11 Bate12 Monuki13 Simó14 Jezewski15 Kriz16 Vinh17 Prussing18 Sun19 Battin20 Sun21 Sun22 Battin23 Boltz24 Moulton25 Sun26 Taff27 Sun28 Gooding29 Loechler30 Sarnecki31 Gooding32 Klumpp33 Klumpp34 Peterson35 Nelson36 Thorne37 Vallado38 Reynolds39 Battin40 Prussing41 Han42 Shen43 Maclellan44 Izzo45 Avanzini46 Arora47 Avendaño48 Bando49 He50 Leeghim51 Zhang52 Arlulkar53 Der54 Luo55 Thompson56 Zhang57 Aretskin58 Parrish59 Wagner60 Ahn61 Arora62 Arora63 Chen64 Usachov65 Arlulkar66 Thorne67 Wailliez68 Wen69 Wen70 Ahn71 Izzo72 1788 1857 1959 1961 1962 1964 1966 1968 1969 1970 1971 1973 1976 1977 1979 1979 1983 1984 1985 1986 1988 1990 1991 1992 1995 1997 1998 1999 2000 2004 2005 2006 2008 2010 2011 2012 2013 2014 2015 David de la Torre Sangrà*, Elena Fantino§ * Space Studies Institute of Catalonia (IEEC), Polytechnic University of Catalonia (UPC), E.T.S.E.I.A.T., Colom 11, 08222 Terrassa (Spain), david.de.la.torre.sangr[email protected] § Space Studies Institute of Catalonia (IEEC), Polytechnic University of Catalonia (UPC), E.T.S.E.I.A.T., Colom 11, 08222 Terrassa (Spain), elena.fan[email protected] Review of Lambert’s problem