Revisiting the computation of the critical points of the Keplerian distance

We consider the Keplerian distance d in the case of two elliptic orbits, i.e., the distance between one point on the first ellipse and one point on the second one, assuming they have a common focus. The absolute minimum dmin\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_{\textrm{min}}$$\end{document} of this function, called MOID or orbit distance in the literature, is relevant to detect possible impacts between two objects following approximately these elliptic trajectories. We revisit and compare two different approaches to compute the critical points of d2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d^2$$\end{document}, where we squared the distance d to include crossing points among the critical ones. One approach uses trigonometric polynomials, and the other uses ordinary polynomials. A new way to test the reliability of the computation of dmin\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_{\textrm{min}}$$\end{document} is introduced, based on optimal estimates that can be found in the literature. The planar case is also discussed: in this case, we present an estimate for the maximal number of critical points of d2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d^2$$\end{document}, together with a conjecture supported by numerical tests.


Introduction
The distance d between two points on two Keplerian orbits with a common focus, that we call Keplerian distance, appears in a natural way in Celestial Mechanics.The absolute minimum of d is called MOID (minimum orbital intersection distance), or simply orbit distance in the literature, and we denote it by d min .It is important to be able to track d min , and actually all the local minimum points of d, to detect possible impacts between two celestial bodies following approximately these trajectories, e.g. an asteroid with the Earth [16,17], or two Earth satellites [19].Moreover, the information given by d min is useful to understand observational biases in the distribution of the known population of NEAs, see [13].Because of the growing number of Earth satellites (e.g. the mega constellations of satellites that are going to be launched [2]) and discovered asteroids, fast and reliable methods to compute the minimum values of d are required.
The computation of the minimum points of d can be performed by searching for all the critical points of d 2 , where considering the squared distance allows us to include trajectory-crossing points in the results.
There are several papers in the literature concerning the computation of the critical points of d 2 , e.g.[20,7,15,9,10,3].Some authors also propose methods for a fast computation of d min only, e.g.[21,14].
We will focus on an algebraic approach for the case of two elliptic trajectories, as in [15], [9].In [15] the critical points of d 2 are found by computing the roots of a trigonometric polynomial g(u) of degree 8, where u is the eccentric anomaly parametrizing one of the trajectories.The polynomial g(u) is obtained by the computation of a Groebner basis, implying that generically we can not solve this problem by a polynomial with a smaller degree.In [9], resultant theory is applied to a system of two bivariate ordinary polynomials, together with the discrete Fourier transform, to obtain (generically) a univariate polynomial of degree 20 in a variable t, with a factor (1 + t 2 ) 2 leading to 4 pure imaginary roots that are discarded, so that we may have at most 16 real roots.Note that the trigonometric polynomial g(u) of degree 8 corresponds to an ordinary polynomial of degree 16 in the variable t through the transformation t = tan(u/2).These methods were extended to the case of unbounded conics with a common focus in [3], [10].
In this paper we revisit the computation of the critical points of d 2 for two elliptic trajectories by applying resultant theory to polynomial systems written in terms of the eccentric or the true anomalies.We obtain different methods using either ordinary or trigonometric polynomials.Moreover, we are able to compute via resultant theory the 8-th degree trigonometric polynomial g(u) found by [15], and its analogue using the true anomalies (see Sections 4,5).Some numerical tests comparing these methods are presented.We also test the reliability of the methods by taking advantage of the estimates for the values of d min introduced in [13] when one trajectory is circular.For the case of two ellipses, since we do not have such estimates for d min , we use the optimal bounds for the nodal distance δ nod derived in [11].
After introducing some notation in Section 2, we deal with the problem using eccentric anomalies and ordinary polynomials in Section 3. In Sections 4 and 5 we describe other procedures employing trigonometric polynomials and, respectively, eccentric or true anomalies.Some numerical tests and the reliability of our computations are discussed in Section 6.Finally, we present results for the maximum number of critical points in the planar problem in Section 7, and draw some conclusions in Section 8.Additional details of the computations can be found in the Appendix.

Preliminaries
Let E 1 and E 2 be two confocal elliptic trajectories, with E i defined by the five Keplerian orbital elements a i , e i , i i , Ω i , ω i .We introduce the Keplerian distance where X 1 , X 2 ∈ R 3 are the Cartesian coordinates of a point on E 1 and a point on E 2 , corresponding to the vector V = (v 1 , v 2 ), where v i is a parameter along the trajectory E i .
In this paper we will parametrize the orbits either with the eccentric anomalies u i or with the true anomalies f i .Let (x 1 , y 1 ) and (x 2 , y 2 ) be Cartesian coordinates of two points on the two trajectories, each in its respective plane.The origin for both coordinate systems is the common focus of the two ellipses.We can write p = (p x , p y , p z ) , q = (q x , q y , q z ) , where If we use the eccentric anomalies u i we have for i = 1, 2, while with the true anomalies f i we have and set K = ⟨P, p⟩, L = ⟨Q, p⟩, M = ⟨P, q⟩, N = ⟨Q, q⟩.

Eccentric anomalies and ordinary polynomials
We look for the critical points of the squared distance d 2 as a function of the eccentric anomalies u 1 , u 2 , that is we consider the system where ∇d 2 = ∂d 2 ∂u 1 , ∂d 2 ∂u 2 .We can write where System (3) can be written as with , Following [9], we can transform (4) into a system of two bivariate ordinary polynomials in the variables t, s through Then, (4) becomes where and be the Sylvester matrix related to (5).From resultant theory [6] we know that the complex roots of det(S 0 (t)) correspond to all the t-components of the solutions (t, s) ∈ C 2 of (5).
The determinant det(S 0 (t)) is in general a polynomial of degree 20 in t.We notice that it can be factorized as det where Therefore, to find the t-components corresponding to the critical points, we can look for the solutions of det( Ŝ) = 0, which in general is a polynomial equation of degree 16.We can follow the same steps explained in [10,Sect. 4.3] to obtain the coefficients of the polynomial det( Ŝ) by an evaluation/interpolation procedure based on the discrete Fourier transform.Then, the method described in [4] is applied to compute its roots.We substitute each of the real roots t of det( Ŝ) in (5) and use the first equation p(t, s) = 0 to compute the two possible values of the s variable.Finally, we evaluate q at these points (t, s) and choose the value of s that gives the evaluation with the smallest absolute value.

Angular shifts
Define a shifted angle v 1 = u 1 − s 1 and let Then, system (4) becomes The coefficients α, β, γ, Ã, B, D are written in Appendix A. If T 0 (z) is the Sylvester matrix related to (10) we get det(T 0 (z We find the values of z by solving the polynomial equation det( T ) = 0, which again has generically degree 16.We compute the values of v 1 from (9) and shift back to obtain the u 1 components of the critical points.Substituting in (4) and applying the angular shift u 2 = v 2 + s 2 , we consider the system where the first equation corresponds to the first equation in ( 4), and For each value of u 1 we compute two solutions for cos v 2 , sin v 2 and the corresponding values of cos u 2 , sin u 2 .We choose between them by substituting in the second equation in (4).

Eccentric anomalies and trigonometric polynomials
To work with trigonometric polynomials, we write system (3) as where Inserting relation into cos 2 u 1 + sin 2 u 1 − 1 = 0 and into the first equation in (12), we obtain We call p 1 , p 2 the two trigonometric polynomials appearing on the left-hand side of ( 14).The Sylvester matrix of p 1 and p 2 is which corresponds to the resultant of p 1 , p 2 with respect to cos u 1 and is a trigonometric polynomial in u 2 only.The u 2 component of each critical point satisfies G (u 2 ) = 0.
Proposition 1.We can extract a factor β 2 from det S .
Proof.Using simple properties of determinants, we can write det S as a sum of different terms.The terms independent from β in this sum are given by , and both determinants are 0. The linear terms in β are given by and this sum is 0, because the two determinants are opposite.Therefore, G (u 2 ) is made by terms of order higher than 1 in β.It results where Their explicit expressions read The trigonometric polynomial has total degree 8 in the variables cos u 2 , sin u 2 , and corresponds to the polynomial g introduced in [15] with Groebner bases theory.For this reason, generically, there is no polynomial of smaller degree giving all the u 2 components of the critical points of d 2 .
We now explain the procedure to reduce the problem to the computation of the roots of a univariate polynomial.We set We find that for some polynomial coefficients g j such that Then, we consider the polynomial system g(x, y) = 0, Using relations obtained from the second equation in (15), we can substitute g in system (15) with where We can also write with Note that a and b have degree 7 and 8, respectively.We eliminate y from system g(x, y) = 0, by computing the resultant u(x) of the two polynomials with respect to y, and obtain which is a univariate polynomial of degree 16.Each of the real roots x of u, with |x| ≤ 1, is substituted into the equation g(x, y) = 0 to get the value of y.Finally, we evaluate α, β, γ, µ, ν at the computed pairs (x, y) and solve system (12) by computing the values of cos u 1 and sin u 1 from ( 14) and (13), respectively.

Finding the roots of u with Chebychev's polynomials
To compute the roots of the polynomial u(x) in a numerically stable way, we need to express u in a basis ensuring that the roots are well-conditioned functions of its coefficients.This can be achieved using Chebyshev's polynomials [18] in place of the standard monomial basis.
In the monomial basis we have for some coefficients p j .The same polynomial can be written as where T j are Chebyshev's polynomials, recursively defined by which are a basis for the vector space of polynomials of degree at most n.The coefficients c j are obtained from the p j as follows.Setting with where N n = 0, that is N is a nilpotent matrix of order n.Relation Let us introduce the vectors made by the coefficients of the polynomials in ( 17), (18), and the diagonal matrix From ( 20) and ( 21) we can write Therefore, the relation between the coefficients c j and p j is given by C = (D Ã−t )P.
Searching for the roots of u(x) corresponds to computing the eigenvalues of an n × n matrix C , called colleague matrix [8].We use the form of the colleague matrix described in [5]: The computation of the roots of a polynomial using the colleague matrix and a backward stable eigenvalue algorithm, such as the QR algorithm, is backward stable, provided that the 2-norm of the polynomial is moderate (see [18]).

True anomalies and trigonometric polynomials
The same steps described in Section 4 can be applied to look for the critical points of the squared distance function expressed in terms of the true anomalies f 1 , f 2 .Note that using true anomalies allows to deal with both bounded and unbounded trajectories [20,10].We write the system ∇d 2 (f 1 , f 2 ) = 0 where We also set Inserting relation into cos 2 f 1 + sin 2 f 1 − 1 = 0 and into the second equation of (23), we obtain As in Section 4, we consider the Sylvester matrix of the two polynomials in (25) and define that we are able to factorize.In particular, we can write where with We can show that h(f 2 ) has degree 8 in (cos We find that h(x, y) = Then, we consider the system h(x, y) = 0, Proceeding as in Section 4 we can substitute h(x, y) with with where We apply resultant theory to eliminate the dependence on y as in Section 4 and obtain a univariate polynomial v of degree 16.The real roots of v with absolute value ≤ 1 correspond to the values of cos f 2 we are searching for.We compute sin f 2 from (28) and substitute cos f 2 and sin f 2 in (25).Finally, cos f 1 and sin f 1 are found by solving (25) and using (24).

Angular shifts
Also for the method presented in this section we consider the application of an angular shift.If we define the new shifted angle v 2 by for some s 2 ∈ [0, 2π), the coefficients of the polynomial (26) written in terms of v 2 are derived following the computations of Appendix C.Then, following a procedure analogous to Section 5, we find the values of v 2 and shift back to get f 2 .Finally, we can apply an angular shift also to the angle f 1 when solving system (25).Defining the shifted angle as for s 1 ∈ [0, 2π), system (25) becomes where with α, β, γ, κ, λ, µ, ν defined at the beginning of this section.

Numerical tests
We have developed Fortran codes for each of the methods presented in this paper.We denote these methods with (OE, OES, TE, TEC, TT, TTS), see Table 1.Moreover, we denote by (OT) the method presented in [10].Numerical tests have been carried out for pairs of bounded trajectories to compare the different methods.
Taking the NEA catalogue available at https://newton.spacedys.com/neodys/,we applied these methods to compute the critical points of the squared distance between each NEA and the Earth, and between all possible pairs of NEAs.We applied a few simple checks to detect errors in the results: • Weierstrass check (W): for each pair of trajectories we have to find at least one maximum and one minimum point; • Morse check (M): for each pair of trajectories, let N be the total number of critical points, and M and m be the number of maximum and minimum points, respectively.Then (assuming d 2 is a Morse function) we must have N = 2(M + m); • Minimum distance check (d min ): we sample the two trajectories with k uniformly distributed points each (we used k = 10), and compute the distance between each pair of points.We check that the minimum value of d computed through this sampling is greater than the value of d min obtained with our methods.
For each method a small percentage of cases fail due to some of the errors above.However, in our tests, for each pair of orbits, at least one method passes all three checks.The angular shifts (see Sections 3.1, 5.1) allow us to solve the majority of detected errors for the methods of Sections 3 and 5 without shift (OE, TT).Applying a shift could also be a way to solve most of the errors detected by the method of Section 4 (TE, TEC).Some data on the detected errors for each method is reported in Table 1.Here we show the percentages of cases failing some of the three checks described above.We note that the NEA catalogue contains 31,563 NEAs with bounded orbits (to the date of February 25, 2023).Therefore, for the NEA-Earth test we are considering 31,563 pairs of orbits, while for the NEA-NEA test the number of total pairs is 498,095,703.From Table 1 we see that the method TE is improved with Chebychev's polynomials (TEC).In the same way, the methods OE, TT are improved by applying angular shifts in case of detected errors (OES, TTS).Indeed, the method TTS turns out to be the most reliable for the computation of d min .
Two additional ways to check in particular the computation of d min are discussed below.

Reliability test for d min
Although all the presented methods allow us to find all the critical points of d 2 , we are particularly interested in the correct computation of the minimum distance d min .For this reason, we introduce two different tests to check whether the computed values of d min are reliable.
The first test is based on the results of [13], where the authors found optimal upper bounds for d min when one orbit is circular.Let us denote with A 1 and A 2 the two trajectories.Assume that A 2 is circular with orbital radius r 2 , and call q 1 , e 1 , i 1 , ω 1 the pericenter distance, eccentricity, inclination and argument of pericenter of where we used q max = 1.3, which is the maximum perihelion distance of near-Earth objects.
We compare this optimal bound with the maximum values of d min computed with the method OE, for a grid of values in the (q 1 , ω 1 ) plane.The results are reported in Figure 1.
Here we see that the maximum values of d min obtained through our computation appear to lie on the grey surface corresponding to the graph of max C d min (q 1 , ω 1 ) defined in (31).
This test confirms the reliability of our computations.Similar checks were successful with all the methods of Table 1.
To test our results also in case of two elliptic orbits, we consider the following bound introduced in [11] for the nodal distance δ nod defined below.Let where q 2 , e 2 are the pericenter distance and eccentricity of A 2 , and ω M 1 , ω M 2 are the mutual arguments of pericenter (see [11]).
We introduce the ascending and descending nodal distances The (minimal) nodal distance δ nod is defined as Figure 1: Comparison between the graph of max (e 1 ,i 1 ) d min (q, ω) defined in (31) and the maximum values of d min computed with the algorithm of Section 3 for a grid of values of (q 1 , ω 1 ).
For each choice of (q 1 , ω 1 ) ∈ D, defined as in (30), we have max where, denoting by Q 2 the apocenter distance of A 2 and by p 2 = q 2 (1 + e 2 ) its conic parameter, .
We compare the computed values of d min with the bound (34) on the maximum nodal distance.The results are displayed in Figure 2 where, for four different values of e 2 , the Figure 2: Comparison of the maximum MOID obtained with the method of Section 3 and the bound on the nodal distance found in [11].These plots were drawn using values of e 2 equal to 0.2 (top-left), 0.3 (top-right), 0.4 (bottom-left) and 0.5 (bottom-right).grey surface represents the bound of [11], while the black dots correspond to the maximum value of d min for a grid in the (q 1 , ω 1 ) plane computed with the method OE.Since the value of δ nod is always greater than or equal to the value of d min , for the test to be satisfied, we need all the black dots to fall below or lie on the grey surface.From Figure 2 we see that this is indeed what happens.
Similar checks done with the methods OT, OE, OES, TT, TTS were successful.

The planar case
Let us consider the case of two coplanar conics parametrized by the true anomalies to the first and second conic at X 1 (f 1 ) and X 2 (f 2 ), respectively, are parallel.If one trajectory, say the second one, is circular, then the tangent vector τ 2 is orthogonal to the position vector X 2 for any value of f 2 .Therefore, to find critical points that do not correspond to trajectory intersections, it is enough to look for values of f 1 such that r 1 • τ 1 = 0.By symmetry, we can write that is we can assume ω 1 = 0. Thus, up to a multiplicative factor, we have Therefore, in general, we have the four critical points ( f1 , f2 ) = (0, 0), (0, π), (π, 0), (π, π).
We may have at most two additional critical points that correspond to trajectory intersections, see [12,Sect. 7.1].In conclusion, the maximum number of critical points with a circular and an elliptic trajectory in the planar case is 6.
We consider now the case of two ellipses.The position vectors can be written as Up to a multiplicative factor, we have The critical points that do not correspond to trajectory intersections are given by the values of f 1 , f 2 such that τ 1 , τ 2 are parallel and both orthogonal to X 2 − X 1 .These two conditions lead to the system Multiplying (35) by r 2 and subtracting (36) we get Equations ( 35) and (37) can be written as    From (38) we obtain which is replaced into relation cos 2 f 2 + sin 2 f 2 − 1 = 0 to give Moreover, after replacing in (40) which follows from (39), we obtain Since α 2 + β 2 = 1 + e 2 1 + 2e 1 cos f 1 , the trigonometric polynomial in (41) has, in general, degree 5 in cos f 1 , sin f 1 .Therefore, we can not have more than 10 critical points which do not correspond to trajectory intersections.Then, the maximum number of critical points of d 2 (including intersections) for two elliptic orbits in the planar case is at most 12.However, we remark that this bound has never been reached in our numerical tests, where we got at most 10 critical points that we think is the maximum number.This conjecture adds a new question to Problem 8 in [1].
In Table 2 we write a set of orbital elements giving 10 critical points.We draw the level curves of d 2 in Figure 3, as function of the eccentric anomalies u 1 , u 2 , where the position of each critical point is highlighted: we use an asterisk for saddle points, and crosses for local extrema.Finally, the critical points, the corresponding values of d and their type (minimum, maximum, saddle) are displayed in Table 3.

Conclusions
In this work we investigate different approaches for the computation of the critical points of the squared distance function d 2 , with particular care for the minimum values.We focus on the case of bounded trajectories.Two algebraic approaches are used: the first employs ordinary polynomials, the second trigonometric polynomials.In both cases we detail all the steps to reduce the problem to the computation of the roots of a univariate polynomial of minimal degree (that is 16, in the general case).The different methods are compared through numerical tests using the orbits of all the known near-Earth asteroids.We also perform some reliability tests of the results, which make use of known optimal bounds on the orbit distance.Finally, we improve the theoretical bound on the number of critical points in the planar case, and refine the related conjecture.

A Coefficients of shifted polynomials for method with ordinary polynomials and eccentric anomalies
The coefficients of system (10) are

B Factorization of H
and define H (f 2 ) = det T .
Proof.We first prove that we can extract the factor β 2 .Noting that we consider We can write det T as a sum of terms where the only one that is independent on β is which is equal to 0, as previously proved at the beginning of Proposition 1.The terms that are linearly dependent on β are given by and their sum is equal to 0, because the two determinants are opposite.Therefore, H (f 2 ) is made by terms of degree higher than 1 in β.Thus, we can write where Remark 1. H (f 2 )/β 2 is a trigonometric polynomial of degree 10 in (cos f 2 , sin f 2 ).

Table 2 :
Cometary orbital elements of a pair of coplanar elliptic orbits giving 10 critical points.Angles are in degrees.

Figure 3 :
Figure 3: Level curves of the squared distance for the case reported in Table2.The position of the critical points is highlighted: saddle points are represented by black asterisks, while the red and blue crosses correspond to maximum and minimum points, respectively.

Table 1 :
Percentages of detected errors with each method applied to the computation of all critical points between all NEAs and the Earth and between all pairs of NEAs.

Table 2 .
The position of the critical points is highlighted: saddle points are represented by black asterisks, while the red and blue crosses correspond to maximum and minimum points, respectively.

Table 3 :
Critical points and critical values for the pair of orbits displayed in Table2.