Abstract

The computation of the Minimum Orbital Intersection Distance (MOID) is an old, but increasingly relevant problem. Fast and precise methods for MOID computation are needed to select potentially hazardous asteroids from a large catalogue. The same applies to debris with respect to spacecraft. An iterative method that strictly meets these two premises is presented.

keywords:

Celestial mechanics — methods: numerical — astronomical data bases

††pubyear: 2018††pagerange: On the minimum orbital intersection distance computation: a new effective method–C

1 Introduction

The MOID is the distance between the closest points of the osculating orbits of two bodies. It can be computed as the minimum distance between two sets of points embedded in the Euclidean 3-space (R3) – see Fig. 1. One of the orbits is considered to be the reference or target (primary). Any other orbit facing the primary is called secondary, although the roles can be exchanged as needed.

MOID computation is a classical but ever relevant problem. The risk of collision of a Potentially Hazardous Object (PHO) such as Apophis with the Earth (Wlodarczyk, 2013) and of space debris with an operational spacecraft (Casanova
et al., 2014) are two important applications.

Since the sets of secondaries are usually large, MOID and its time derivative can be used as a prefilter – for example, the second level filter proposed by Hoots (Hoots
et al., 1984) – to discard those that will not give problems in the immediate future. Those selected can be studied with more sophisticated methods – see, for example, Bonnano (2000).

In most MOID calculations, orbits can be considered as Keplerian ellipses. Some bodies originate outside of the Solar System. Recently the interstellar object 1I/2017 U1, later named Oumuamua, reached its perihelion (0.2483ua) in a hyperbolic trajectory. In this case, MOID calculation with Earth involves a hyperbola. Thus, the main task is to compute the minimum distance between two confocal Keplerian orbits, the primary – Earth or spacecraft – and one or more secondaries – near-Earth objects or space debris.

Current methods for MOID computation are sufficient for identifying NEOs as Potentially Hazardous Asteroids, or space debris coming close to a spacecraft orbit. Faster algorithms for MOID computation could help to identify NEOs coming close to perturbing planets (Sitarski, 1968), which in time could change their orbits significantly; study evolution of perturbed NEO’s MOID; manage growing catalogues of space debris; and compute covariance of NEO’s and debris’ MOID.

There is an abundant literature on MOID computation from the mid-20th century to the present day. Some authors use an algebraic approach, obtaining all critical points of the distance function. Others use numerical methods, finding the global minimum by iterations. Others use a hybrid approach, such as Derevyanka (2014). Some key results are summarized in Table 1.

Evaluating the collision risk of asteroids or comets with the Solar system planets

Application of geometric algebra techniques to search for the real roots of a univariate polynomial of degree 16 that gives the critical points of the distance between two orbits using Fast Fourier Transform.

Study of parabolic comet approaches to major planets and their possible capture.

Numerical resolution by bisection method of the two gradient cancellation equations of the three-dimensional distance for seeking stationary points, followed by the type checking of the critical point and the selection of the global minimum.

Classification of asteroids as potentially hazardous (PHA) by estimating the MOID

A numerical method based on approximating the actual distance between orbits by the distance measured exclusively in each ecliptic meridian plane and using parallel programming to accelerate the computation time of the MOID.

Table 1: Timetable of some MOID computation methods.

The Space Dynamics Group (SDG) computation method of the MOID (from now on SDG method) can be classified as numerical iterative and it is based on the following two algorithms:

Algorithm 1:

For a given point P0, the minimum distance d to the primary ellipse e1 is computed (see Fig. 2).

Algorithm 2:

As we move P0 along the secondary ellipse e2, by using Algorithm 1 we obtain a set S of minimum distances to e1. The absolute minimum of set S is then selected, which is the MOID (see Fig. 1).

Both algorithms have strong theoretical foundations, as shown below.

A preliminary version of this work was presented at 2016 Stardust Final Conference on Asteroids and Space Debris (Hedo
et al., 2016).

The remainder of this article is organized as follows: section 2 gives a more detailed description of both algorithms of the SDG method, starting with its mathematical foundation. Section 3 deals with the calculations made during the development of the software, which justify the choice of its individual components. Section 4 includes numerical tests for execution time and accuracy. All results are compared with those of a well accredited method (Gronchi, 2005). Indirect comparisons with other methods are also carried out. Finally, section 5 presents the conclusions.

2 Mathematical foundations of the SDG method

2.1 Distance between two confocal ellipses

Let e1,e2⊂R3 (Euclidean 3-space) be two ellipses with a common focus F (Fig. 1). The MOID is the minimum Euclidean distance between these ellipses considered as sets of points

d(e1,e2)=minE1∈e1,E2∈e1d(E1,E2)

(1)

Distance computation is simpler if we take eccentric anomalies u1 and u2 of each ellipse as variables:

d(e1,e2)=minu1,u2∈[0,2π]d(E1(u1),E2(u2))

(2)

This requires transferring vectors between four reference frames:

∙Fx0y0z0, the common inertial reference frame, with origin in focus F and orthogonal directions given by right-handed vectrix [i0j0k0]. Classical element sets aj,ej,Ωj,ij,ωj of both orbits refer to this frame (not shown in figures).

∙Fpjqjhj, perifocal reference of orbit j, with directions given by the right-handed vectrix [upjuqjuhj]: upj towards pericentre and uhj along angular momentum vector. That of orbit 1 is represented in Fig. 1.

∙Ox1y1z1, central frame of primary ellipse e1 (see Figs. 1 and 2), with origin O in the geometric centre of the ellipse, and directions parallel to the perifocal frame of this orbit.

Rotation matrix {Q}0j from reference 0 to reference j is

[upjuqjuhj]

=[i0j0k0]{Q}0j=

=[i0j0k0]{R}z(Ωj){R}x(ij){R}z(ωj)

(3)

where {R}z(α) is the rotation matrix of angle α about the axis Oz, in the classical Euler angle sequence.
The relative rotation matrix from the perifocal reference of e2 to that of e1 is

[up2uq2uh2]=[up1uq1uh1]{Q}12{Q}⊺01%Q02

(4)

{Q}12={R}⊺3(ω1){R}⊺1(i1){R}2(Ω2−Ω1){R}1(i2){R}2(ω2)

(5)

In the central frame of e1, the position vector of point E2∈e2 is OE2=FE2+OF, and its components are

(6)

The Euclidean distance between points E1 and E2 is

d(u1,u2)

=|OE1(u1)−OE2(u2)|=

=
⎷3∑i=1(xE1i(u1)−xE2i(u2))2

(7)

Therefore, MOID computation is an unconstrained minimization of the function of two variables given by (7). The geometric meaning of the eccentric anomaly makes the function d(u1,u2)2π-periodic in each of its two arguments, so the final domain will be restricted to (u1,u2)∈[0,2π]×[0,2π] (a topological torus).

\pstThreeDPut(0,0,0)O\pstThreeDLine[linewidth=.5pt,linecolor=black,arrows=->](-2.5,0,0)(5,0,0)\pstThreeDPut(5,0,0)x1\pstThreeDLine[linewidth=.5pt,linecolor=black,arrows=->](0,-1.5,0)(0,5,0)\pstThreeDPut(0,5,0)y1\pstThreeDLine[linewidth=.5pt,linecolor=black,arrows=->](0,0,0)(0,0,3)\pstThreeDPut(0,0,3)z1\pstThreeDEllipse[linewidth=0.5pt,linecolor=blue](0,0,0)(4.25,0,0)(0,4,0)\pstThreeDPut(1.436,0,0)\pstThreeDDot[drawCoor=true,linecolor=red](-2,-0.9,2.35)\pstThreeDPut(-3,-1.8,2.4)P0\pstThreeDPut(-3,-1.8,0)P\pstThreeDPut(-3.7,-1.95,0)E∗\pstThreeDLine[linecolor=brown](-2,-0.9,0.0)(-3.7,-1.95,0)\pstThreeDPut(-3.03,-0.83,0)d∥\pstThreeDLine[linecolor=darkgreen](-2,-0.9,0.0)(-2,-0.9,2.35)\pstThreeDPut(-2,-0.55,1.25)d⊥\pstThreeDLine[linecolor=black](-3.7,-1.95,0)(-2,-0.9,2.35)\pstThreeDPut(-2.8,-1.2,1.25)d\pstThreeDPut(-3,2,0)e1Figure 2: Distance from a point to an ellipse in R3.

2.2 Distance between a point and an ellipse

Let P0∈R3 be a point, e⊂R3 an ellipse, and u∈[0,2π] its eccentric anomaly. Distance is defined as

d(P0,e)=minu∈[0,2π]d(P0,E(u))=minu∈[0,2π]d(u)

(8)

Let P be the orthogonal projection of P0 onto the plane of e (see Fig. 2). P0E∗ can be divided into components normal to the orbit plane and parallel to it

d(P0,e)2=d2⊥+d2∥=|P0P|2+|PE∗|2

(9)

where E∗ is the point of e which gives minimum distance. Note that the plane P0PE∗ is normal to e at E∗; this property was exploited by Milisavljevic (2010) to compute MOID.

Keeping P0 fixed, the out-of-plane distance d⊥ is constant. Therefore only the in-plane distance has to be minimized. This is the core of Algorithm 1 which provides the point E∗ where in-plane distance reaches its minimum value d∥.

Since Ox1y1z1 are symmetry axes of the ellipse, distance to point P0 is invariant to symmetries respect to coordinate planes. We only consider the first octant of Ox1y1z1. The same applies to the projection P: we only need to consider the first quadrant.

2.3 Distance between an ellipse and a coplanar point: in-plane distance

Let E be a point of ellipse e, and P a point contained in its plane. The minimum distance is

d(P,e)=minE∈ed(P,E)

(10)

and it can be computed in Cartesian coordinates or with parametric equations. The first approach gives useful insight on the number and distribution of distance minima, while the second is better suited for numerical minimization.

2.3.1 First approach: Cartesian coordinates

Let Oxy be the symmetry axes of ellipse e (see Fig. 3), with semiaxes a and b, parallel to Ox and Oy, respectively. Let (α,β) be the Cartesian coordinates of point P in Oxy. The minimization problem to solve is

d2(P,e)=min(x,y)∈e[(x−α)2+(y−β)2]

(11)

subject to x2a2+y2b2=1

(12)

We use the method of Lagrangian multipliers. Eliminating the multiplier from the objective function derivatives, we obtain

xy(a2−b2)−αa2y+βb2x=0

(13)

which is the implicit Cartesian representation of an equilateral hyperbola h of asymptotes parallel to the axes, passing through P and O. This property is known since the third century BC (Apollonius of
Perga, 1896).

eOxyP(α,β)Ed(P,E)abc=aeFaFigure 3: In-plane distances between an ellipse and a coplanar point.

Each solution of the system (12, 13) is an intersection point of e and h, where the distance to P has a stationary value.
The intersections of two second order curves consist of 0, 2 or 4 points (counted with their multiplicities). Fig. 4 shows an example with the relative disposition of both curves. Notice that

e has its centre in the origin O(0,0)

h is an equilateral hyperbola. The asymptote parallel to Oy has the equation x=αe2>α≥0; the asymptote parallel to Ox has the equation y=−β1−e2e2≤0.

where e is the eccentricity of the ellipse.

The upper hyperbolic branch comes down from infinity through P(α,β), through the origin O, and again to infinity towards the left. Thus h and e cut each other at least in two points, and system (12, 13) always has at least two real roots: the problem of extreme distances always has solution (Armellin
et al., 2010). Besides, conics are simple curves, without inflection points, so this branch can only intersect the ellipse at two points.

The upper hyperbolic branch is located to the left of the vertical asymptote and above the horizontal one; the same happen with points P and O that belong to this branch. As a consequence,

in the limit x→(αe2)− the upper branch of h approaches the vertical asymptote. By continuity, it cuts the ellipse e in a single point E∗ of the first quadrant. This is the global minimum of the distance.

in the limit x→−∞ the upper branch of h approaches the horizontal asymptote. By continuity, it cuts the ellipse e in a single point E∗∗ of the third quadrant. This is the global maximum of the distance.

The lower hyperbolic branch is located to the right of the vertical asymptote, and below the horizontal one. It can only cut the ellipse in the fourth quadrant: two cuts, one (tangent), or none. Fig. 4 shows a case with two cuts, E− and E+. The first one gives a relative minimum and the other a relative maximum, unless both join into a single tangent point. Figs. 8–8 show the three possibilities. In any case, they satisfy

d(P,E∗)

<min(d(P,V+),d(P,W−))<d(P,E−)≤

≤d(P,E+)<d(P,E∗∗)

(17)

The intersection points of e and h can be obtained by eliminating y between Eqs. (13) and (12):

x4e4−2αe2x3−(a2e4+β2e2−[α2+β2])x2++2αe2a2x−a2α2=0

(18)

Eq. (18) is a quartic with real coefficients, which in general can have 0, 2, or 4 real roots (counting the multiplicity index). As shown by the intersections with the hyperbola, in this case there are at least 2 real roots. The exact number can be determined by geometric considerations.

It is well known that, for critical distances, PE is normal to the ellipse e and tangent to its plane evolute p – see Milisavljevic (2010). Since all normals are tangent to the evolute, there are as many real roots of (18) as tangents to p from P. The number depends on the position of P relative to p.

Substituting the coordinates (α,β) of P into Π shows its relative position.

If P is inside p (Fig. 8), Π(α,β)<0 and there are four different tangents to p, unless P is upon one of the axes. In this case there are only three, but one is the axis itself, which is tangent at two points. Therefore, there are four roots of Eq. (18).

If P is on p (Fig. 8), Π(α,β)=0 and there are only three different tangents from P. So the number of different real roots of Eq. (18) is three but one, the tangent at P itself, is double (two tangents to the same branch join into one).

If P is outside p (Fig. 8), Π(α,β)>0 and only two tangents to p can be traced from P; then the number of different real roots of Eq. (18) is two.

Two consequences can be drawn from the previous results:
(a) the problem always has at least one solution with minimum distance and one solution with maximum distance; and (b) there may be two different real roots that give minimum distance because of the symmetries.

Figure 7: Tangent lines to the plane evolute of an ellipse from a point on the evolute itself, Π(α,β)=0.

Figure 6: Tangent lines to the plane evolute of an ellipse from an interior point, Π(α,β)<0.Figure 7: Tangent lines to the plane evolute of an ellipse from a point on the evolute itself, Π(α,β)=0.Figure 8: Tangent lines to the plane evolute of an ellipse from an exterior point, Π(α,β)>0.Figure 6: Tangent lines to the plane evolute of an ellipse from an interior point, Π(α,β)<0.

2.3.2 Second approach: parametric representation

Consider the parametric equations of the ellipse e in terms of the eccentric anomaly u∈R, as shown in Fig. 9

where u is the independent variable and a,b,α,β are the four parameters of the problem. For any parameter values, the function d(u) is:

2π-periodic; as a consequence u∈[0,2π];

always non-negative, zero if and only if P∈e.

Using the squared distance simplifies the computation of stationary points u∗, since the square root is avoided, and the stationary points and their character are preserved (see appendix A). Stationary points are related to zeros of the first derivative.

We introduce the univariate function

f(u;a,b,α,β)

=d[d2(u;a,b,α,β)]du=

=αasinu−βbcosu−c2sinucosu

(22)

whose real roots provide the eccentric anomalies of the points of ellipse e with stationary distance to P(α,β). These are the same points obtained from Eq. (18), but here the variable is the parameter u.

2.3.3 Study of minimum distance cases

Let u∗ be any real root of Eq. (22), that is, an extremum of distance. Let u∗ be a real root with global minimum distance (there could be more than one due to symmetries). From now on we will assume that P is in the first quadrant of Oxy, that is, α≥0, β≥0. We can consider two types of solutions:

Immediate closed-form solutions, such as a circular orbit, or P upon one of the principal axes. These are studied in appendix B, and summarized in Table 2.

Those requiring the application of classic quartic polynomial solution, or numerical methods.

The coefficients in m3 and m of Eq. (28) can become very large when β is very small. In that case, a change of variable improves numerical accuracy. If the relations

sinu=1−q21+q2,cosu=2q1+q2,q=tan(π4−u2)

(30)

are introduced into (22) and squared, the following equation is obtained

q4+(~Λ−~Σ)q3+(~Λ+~Σ)q−1=0

(31)

where the positive parameters (~Λ,~Σ) are given by

~Λ=2bβaα=4Λ,~Σ=2ae2α=12Σ~Λ

(32)

The classic solution can be applied to any of the above quartic polynomials.

In general, the quartic equations (18), (25), (26), (28), and (31) have four complex solutions that can be obtained by explicit formulas, but the useful solutions are real.

The analytical solution of a quartic equation has been coded and tested, but it is slower than the iterative approach, and less accurate (see Table 3 in subsection 3.1 below). Therefore this approach was discarded.

2.3.5 Numerical iterative method

When P(α,β) is in the first quadrant and outside of the coordinate axes, Eq. (22) gives

f(0)=−βb<0;f(π2)=αa>0

Since f(u) is a continuous function in the closed interval [0,π2], Bolzano’s theorem assures that there is, at least, one root u∗∈[0,π2].

As shown in section 2.3.1, there is only one root u∗ in the open interval (0,π2). It is a simple root because f′(u∗)≠0 and it defines the ellipse point with minimum distance to P, so f′(u∗)>0.

Therefore, the problem of computing the minimum distance between an ellipse e and a coplanar point P(α>0,β>0) is reduced to finding the root of the univariate function f(u) in the open interval (0,π2), because the other cases (α=0 and/or β=0) have closed-form trivial solutions.

is non-linear and smooth (of class C∞ in its definition domain). Its derivatives can be easily computed, so high-order iterative methods can be used.

Being transcendental, an iterative root-finding method must be used. Convergence is assured: Brouwer’s fixed-point theorem states that, for any continuous function f mapping a compact convex set to itself, there is a point x0 such that f(x0)=x0 – see Agarwal
et al. (2009). To apply the fixed-point theorem, we rewrite it as

u

=arctan(β√1−e2+ae2sinuα)=Φ(u)

(34)

Φ(u) satisfies the two conditions of Brouwer’s theorem:

Φ(u) is obviously continuous: Φ(u)∈C∞

Φ([0,π2])⊆[0,π2].

To prove the second statement, we see first the trivial cases (aeα=0), where Φ is constant and the fixed point uo is:

If α=0(P∈Oy):Φ(u)=π2(const.)⟹u0=π2

If e=0(e circle):Φ(u)=arctan(βα)(const.)⟹m⟹u0=arctan(βα)

If a=0(e≡O):Φ(u)=arctan(βα)(const.)⟹m⟹u0=arctan(βα)

For non trivial cases ∀α>0,β≥0,a>0,e>0, we must study Φ and its derivative:

Instead of simply iterating u=Φ(u), a high convergence method is more efficient. Some have been considered: Newton–Raphson, Halley’s (Weisstein, 2018), and Improved Newton–Raphson (Danby, 1992). The last two are always faster. However, the specific machine determines which is the best, because of a balance between speed of convergence and number of function evaluations. A quantitative comparison will be shown in subsection 3.2. As an example, appendix C collects Halley’s method characteristics, with its advantages and the pitfalls to avoid.

2.4 SDG method

The approach used in this paper transforms the complex problem of minimization in two variables with ten parameters into a set of simpler minimization problems, each in one variable, with easier mathematical treatment. The objective is speeding up computation.

The main idea of the method is to separate the two components of the distance between a point and an ellipse in Euclidean 3-space: normal to the ellipse plane and contained in it. The procedure is as follows:

Algorithm 1:

finds minimum distance of an arbitrary point P0 to the primary orbit e1.

P0 is projected onto the orbital plane of e1 as P

The minimum distance from P to e1 and the corresponding eccentric anomaly are found. It uses Halley’s iterative method to find the root, which corresponds to the minimum value.

Finally, in-plane and normal components are combined to find the distance.

Algorithm 2:

MOID computation is reduced from a two-variable problem to a one-variable function minimization. Steps are:

Discretise secondary orbit e2, call Algorithm 1 to compute minimum distance to e1 for each anomaly, and bracket local minima.

Refine by a closed interval search algorithm.

Select global minimum between the former values.

The detailed steps are: define a uniform grid of N+1 values {u02,u12,…,uN2} in the interval u2∈[0,2π] where u2 is the eccentric anomaly of e2. The constant interval size is Δu2=2πN. A corresponding grid of N image points {E02,E12,…,EN−12} is defined in e2. The last value, uN2, is ignored because E02≡EN2. The choice of N will be discussed in subsection 3.4.

For each u2, we compute the minimum distance to e1 through Algorithm 1 – as described in subsection 2.2 – to obtain a reticular approximation of the function

Φ(u2)=minu1∈[0,2π]d(u1,u2)

(37)

which is used to bracket the abscissas of relative minima of d(u1,u2) within an interval.

The abscissa of each local minimum is determined by decreasing the size of the uncertainty interval through golden section search. The smallest local minimum is the global one, that is, the MOID.

3 Software development

The NEODyS database corresponding to epoch March 23rd, 2018 = MJD 58200, from the Near Earth Objects - Dynamic Site, has been used to compute the Earth MOID.

Development includes the selection of the best algorithms for each component, and optimisation of free parameters.

3.1 Algorithm 1: algebraic vs iterative

First, a study was made on how to calculate the distance in the plane. The approaches described in subsections 2.3.4 and 2.3.5 have been compared. The complete database has been computed 20 times to minimize the effect of background machine load, and find the average unit time. Test results show:

The algebraic method is 66 per cent slower than the iterative approach, as shown in Table 3. This may be due to the arithmetic of complex numbers, frequent calls to slow functions, and many conditional bifurcations on the code.

In a few cases, the characterization of a root as minimum or maximum fails at double precision. The accumulation of rounding and truncation errors leads to wrong negative radicands and other problems. Quadruple precision has been necessary to accurately cover all cases, making the algebraic method even slower.

3.2 Iterative procedure: Halley vs Danby

We still have to choose the iterative method. Low order ones such as Newton–Raphson have been outright discarded. As mentioned above, Halley’s (Weisstein, 2018) and Danby’s (Danby, 1992) methods have been tested for incorporation into the SDG code. As in the previous section, the whole database of NEOs has been computed 20 times, to find the unit average time.

Comparisons have been carried out in an Intel I7-4771 PC @ 3.5GHz with 16GiB of RAM. Table 4 shows mean times for each MOID in the database, and Halley’s method is slightly faster. Later on, when running SDG method against Gronchi’s in a different machine, Danby’s was ahead (see Table 10 in subsection 4.1 below).

Results show that the most appropriate iteration method depends substantially on the technical characteristics of the computer used. Danby’s method uses fewer iterations, but more function calls for each. Conversely, Halley’s needs more iterations, but with fewer function calls. Thus, a higher frequency favours Halley, but differences in architecture may favour Danby. It may be necessary to test both on the target machine.

Table 5 shows the iterations needed to compute the MOID of the first 30 database objects, with seed arctan(αβ), in the same machine as Table 4.

3244 Halley iterative processes are needed, with an average of 1.73 iterations per process. In most cases (72.66 per cent) two iterations are needed while in a very small part (0.31 per cent) up to three iterations are needed.

3242 Danby iterative processes are needed, with an average of 1.005 iterations per process. In most cases (99.51 per cent) one iteration is needed while in a very small part (0.49 per cent) up to two iterations are needed.

Comparing the columns of Table 5 the advantage of Danby’s method over Halley’s in convergence speed is clear.

Average unit time

µs

Samples

Halley

Danby

#1

33,73

35,53

#2

33,76

35,77

#3

33,71

35,83

#4

33,22

35,65

#5

33,10

35,53

#6

33,22

35,58

#7

33,09

35,56

#8

34,30

35,66

#9

33,72

35,56

#10

33,26

35,62

mean

33,511

35,629

st.dev.

0,394

0,102

Rel. diff to Halley

0,0632

Sample size

20×17648

20×17648

Table 4: Results of unit execution times of the two iterative procedures tested for Algorithm 1 in the SDG method: Halley vs Danby.

Halley

Danby

No.

Abs.

Rel.

Abs.

Rel.

Iterations

Freq.

Freq.

Freq.

Freq.

1

877

0.2703

3226

0.9951

2

2357

0.7266

16

0.0049

3

10

0.0031

0

0.0000

TOTAL :

3244

1.0000

3242

1.0000

AVERAGE :

1.733 iter/proc

1.005 iter/proc

Table 5: Number of Halley and Danby iterative processes and number of iterations for the 30 first objects of NEODyS database.

3.3 Iteration tolerance of Algorithm 1

One or more tolerance values must be selected to end the iterations. Two examples of stop conditions are

|xi−xi−1|/|xi|<xTOL;|fi−fi−1|/|fi|<yTOL;

(38)

The values of these tolerances determine the error admitted and must obviously be greater than the machine epsilon (ϵ): xTOL>ϵ,yTOL>ϵ. In C++ floating point arithmetic with double precision ϵ=2−52≈$2.22×10−16$.

These values influence the average computation times. In SDG method, xTOL=yTOL has been chosen for simplicity, and some studies have been carried out to measure the influence of this value on computation times. Table 6 collects the most significant results. Since the number of iterations carried out by SDG method is small, the influence of these tolerances is very low as can be seen in Table 6.

Tolerance

Unit Time

Non-dim.

µs

MEAN

STD. DEV.

1.00×10−14

30.6761

0.198927374

7.50×10−15

31.1439

0.373606284

5.00×10−15

31.5604

0.391580615

2.50×10−15

31.9417

0.182431753

1.00×10−15

32.9995

0.375676057

5.00×10−16

33.5142

0.251684459

Table 6: Influence of error tolerance on average computation times.

Finally, xTOL=2ϵ≈4.44×10−16has been adopted, more restrictive than the value of 1×10−14chosen by other authors.

3.4 Grid size of Algorithm 2

It is well known – see Gronchi (2002), Milisavljevic (2010), or Wisniowski &
Rickman (2013) – that there can be up to 4 minimum values of the distance between the points of two confocal ellipses. For the method to be fast and accurate, grid size N has to be the smallest number that allows to detect all minima. The number of Earth-distance minima in the full NEODyS database has been computed for increasing values of N. For N≥48, all minima – shown in Table 7 – are detected. A more conservative N=50 has finally been chosen.

No. Minima

No. NEOs

1

8186

2

9444

3

18

4

0

Total:

17648

Average:

1.537 min/object

Table 7: Number of NEOs with a certain number of minimum values of the distance between its orbit and the Earth orbit.

The code is able to redo this parametric analysis if the database is updated.

3.5 Algorithm 2: golden section search iterations

Algorithm 2 needs to find the value and anomaly of bracketed local minima. Golden section search has been selected to refine the values. No comparisons have been made, since we consider it clearly superior to other methods, such as bisection and ternary search.

The choice of grid divisions N affects the number of iterations needed to locate the minimum to required precision. A finer grid needs fewer iterations. But N has other requirements: the lowest value that allows to locate all minima.

Golden Section Search

No.

Abs.

Rel.

Iterations

Freq.

Freq.

28

1

0.02128

29

1

0.02128

30

2

0.04255

31

9

0.19148

32

17

0.36170

33

5

0.10638

34

5

0.10638

35

2

0.04225

36

4

0.10638

37

1

0.08511

TOTAL :

47

1.0000

AVERAGE :

32.467 iter/proc

Table 8: Golden Section algorithm calls and number of iterations for the 30 first objects of NEODyS database

Table 8 shows the number of iterations needed to refine all bracketed local minima. This has been done for the first 30 database objects, on the Intel I7 computer aforementioned.

Golden section search has been called 47 times – as many as minima in the set of 30 NEOs. The number of iterations per call ranges from 28 to 37, with an approximate average of 32.5.

4 Comparison of methods

4.1 Gronchi’s MOID method

Giovanni Gronchi’s FORTRAN software for MOID computation of 1995 (Gronchi, 2005) has been provided by the author’s kindness. It is highly regarded for its speed and precision.

Two computer systems have been used to compare SDG and Gronchi’s methods:

No parallelization has been implemented, only one core is used on each machine. Earth MOID has been computed for the 17648 asteroids in the set. Computations have been repeated at different times to discount machine load.

Tables 9 and 10 show average execution times per MOID for the two computers, each with its set of compilers: C++ for SDG and FORTRAN for Gronchi’s.

Average unit time

µs

Run

Gronchi

SDG

#1

168.43

134.16

#2

170.42

136.09

#3

167.69

135.83

#4

167.55

135.14

#5

167.82

134.20

#6

167.58

134.29

#7

168.25

134.49

#8

167.27

133.93

#9

166.16

135.33

#10

168.96

136.73

mean

168.01

135.02

st.dev.

1.13

0.96

Diff. rel. to Gronchi

0.1964

Sample size

15000

20×16748

Table 9: Average MOID computation times in the Fujitsu PC. Run ten times at different moments.

Average unit time

µs

Run

Gronchi

SDG

Halley

Danby

#1

82.78

64.72

61.18

#2

77.47

64.01

68.79

#3

77.96

68.53

63.57

#4

81.19

61.62

60.47

#5

81.23

68.53

67.73

#6

81.59

63.13

61.00

#7

79.24

63.48

61.54

#8

77.46

61.89

60.83

#9

78.97

63.84

62.07

#10

77.51

65.69

61.09

Mean

79.54

64.54

62.83

St. dev.

2.00

2.42

3.00

Diff. rel. to Gronchi

0.1885

0.2101

Sample size

20×16748

20×16748

20×16748

Table 10: Average MOID computation times in the Supermicro server. Run ten times at different moments.

As for accuracy, the results of SDG and Gronchi’s method are essentially the same. They have been compared asteroid by asteroid, but the spreadsheet would be too large to reproduce. Table 11 shows the eleven NEOs with the greatest differences. Even among those, the difference is less than 1.1×10−15 except for the first entry (2016VB1), which seems to be anomalous. This needs clarification.

Throughout the development of the SDG method, a reference tool was needed to check the results. The classic two-variable minimisation shown in subsection 2.1 has been coded in Maple™ 2017.3. Minimisation has been performed over the [0,2π]×[0,2π] anomaly domain, with 30-digit precision. Maple routines were used. This method is very accurate, but slow.

To solve the discrepancy, global minimum for 2016VB1 has been computed with Maple. The result is the same as that of SDG’s method to 16 digits, the latter’s accuracy. The complete values, primary and secondary anomalies and distance, are shown in Table 12.

Additionally, a contour map of the distance as a function of anomalies – as in Armellin
et al. (2010) – has been generated with Maple (see Fig. 10). This allows to see all extrema and their character. The minima found by SDG and Gronchi’s methods are marked on the plot. These are two very close minima, in a narrow valley, one smaller than the other. It is easy to slip from one to the other.

NEO’S NAME

MOID\earth

ΔMOID\earth

ua

ua

GRONCHI

SDG

|SDG-GRONCHI|

2016VB1

0.0019763068271733

0.0012415773444360

0.0007347294827373

2014SE

0.0342838294450500

0.0342838294450511

0.0000000000000011

159560

0.4733258977758740

0.4733258977758730

0.0000000000000011

2005QO11

0.3352818313812930

0.3352818313812920

0.0000000000000011

2005TS45

0.2892358593469260

0.2892358593469250

0.0000000000000011

2006AV2

0.2668343096980130

0.2668343096980140

0.0000000000000011

2010BH2

0.3557077619773020

0.3557077619773030

0.0000000000000011

2012BD27

0.4626880416478720

0.4626880416478730

0.0000000000000011

2012TA79

0.3278131305485040

0.3278131305485050

0.0000000000000011

2013GY73

0.3041646227553150

0.3041646227553160

0.0000000000000011

2017FO2

0.2802988699612530

0.2802988699612520

0.0000000000000011

Table 11: List of NEOs with the greatest MOID differences between Gronchi’s and SDG methods.Figure 10: Maple contour plot of distances between the Earth and the asteroid 2016VB1. Red is larger and green smaller. Minimum point symbols: SDG (+), Gronchi (∘).

Table 12: Two different minima of the distance between the Earth and the asteroid 2016VB1 orbits with the corresponding anomalies.

Remarkably, only one out of 17648 cases shows discrepancy between both methods. For all others, differences are equal to or smaller than 1.1×10−15ua. Note that both methods use double precision floats, and SDG’s tolerance is xtol=4.44×10−16.

4.2 Other MOID methods

Since we have no access to codes for other methods, direct comparisons are not possible. Two authors include tests against Gronchi’s method, both in speed and accuracy, which allow some indirect conclusions.

Wisniowski &
Rickman (2013) report average unit times of 400µs per MOID for Gronchi’s methods and from 90µs to 150µs for their method. But no indication is given about the equipment with which calculations have been performed.

Derevyanka (2014) obtains average unit times of 120µs without activating parallelization, versus 400µs for Gronchi’s. They use a PC with Intel Core i7-4800MQ CPU @ 2.7GHz, 8.0GiB of DDR3 RAM, and a 1TB HDD at 7200min−1.

The results match Gronchi’s up to at least 9 exact figures, so it could be inferred that a tolerance of about 5×10−10 is being used. Derevyanka’s code is not available. If this tolerance is set in SDG method and a computer similar to Derevyanka’s is used – that mentioned in subsection 3.2 – average unit times are less than 23.5µs.

We cannot directly compare accuracy, but both authors acknowledge that their speed involves some loss of precision. Since our results match Gronchi’s to 15 exact figures, we can conclude that SDG is more precise than the two mentioned above.

5 Conclusions

The mathematical justification given in section 2 makes unnecessary the exhaustive calculation of all critical points of the distance between Keplerian orbits to obtain MOID. The method described in this paper computes MOID faster than the other methods considered, with a smaller risk of omitting or confusing minima.

Depending on the computer used, this method is between 19 and 20 per cent faster than that in Gronchi (2005), as Tables 9 and 10 show. It is not a dramatic improvement, but it also gives an equal or superior precision in locating global minima, which is a common weakness of some methods for multivariate optimization.

Acknowledgements

This work is part of the project entitled Dynamical Analysis of Complex Interplanetary Missions with reference ESP2017-87271-P, supported by Spanish Agencia Estatal de Investigación (AEI) of Ministerio de Economía, Industria y Competitividad (MINECO) and by European Found of Regional Development (FEDER).

We are grateful to Prof. Roberto Armellin for his detailed and helpful review of the manuscript.

The roots of f(u) are u∗∈{0,π2,π,−π2}. Both roots u∗=±π2 give the same minimum distance d(±π2)=b, because f′(±π2)>0 (relative minima) and f′(0)<0,f′(π)<0 (relative maxima). It is easy to see that d(0)=d(π)=a>b=d(±π2).

Subcase α=0,β>0 (P in the positive side of Oy)

The hyperbola degenerates into its asymptotes. Equations (13), (12), and (11) give