# Orbit determination from one position vector and a very short arc of optical
observations

Erica Scantamburlo<EMAIL_ADDRESS>Department of Mechanical and
Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24,
10129, Italy  Giovanni F. Gronchi<EMAIL_ADDRESS>Dipartimento di Matematica, University of Pisa, Largo B. Pontecorvo 5, 56127,
Italy  Giulio Baù<EMAIL_ADDRESS>Dipartimento di Matematica, University
of Pisa, Largo B. Pontecorvo 5, 56127, Italy

###### Abstract

In this paper we address the problem of computing a preliminary orbit of a
celestial body from one topocentric position vector and a very short arc (VSA)
of optical observations. Using the conservation laws of the two-body dynamics,
we write the problem as a system of 8 polynomial equations in 6 unknowns. We
prove that this system is generically consistent, namely it admits solutions
at least in the complex field. From this system we derive a univariate
polynomial $\mathfrak{v}$ of degree 8 in the unknown topocentric distance at
the mean epoch of the VSA. Through Gröbner bases theory, we show that the
degree of $\mathfrak{v}$ is minimum among the degrees of all the univariate
polynomials solving this problem. The proposed method is relevant for
different purposes, e.g. the computation of a preliminary orbit of an Earth
satellite with radar and optical observations, the detection of maneuvres of
an Earth satellite, and the recovery of asteroids which are lost due to a
planetary close encounter. We also show some numerical tests in the case of
asteroids undergoing a close encounter with the Earth.

Keywords: Orbit determination – Keplerian integrals – Algebraic methods

##  1 Introduction

The problem of computing the orbit of celestial bodies has attracted the
interest of scientists since a long time, see Laplace (1780), Lagrange (1783),
Gauss (1809). The recent improvements in the observation technology have posed
new interesting mathematical problems in this field. The number of asteroids
observations performed by modern telescopes is very large. Usually they can be
grouped into very short arcs (VSAs) of optical observations, however it is not
easy to determine whether VSAs collected in different nights belong to the
same observed objects. In general, from a VSA we can not compute a reliable
preliminary orbit, but we can try to put together different VSAs to perform
this task, see e.g. Milani et al. (2005). Assuming that two VSAs belong to the
same asteroid, we can write polynomial equations to compute a preliminary
orbit using the conservation laws of the two-body dynamics. Recently,
different ways to combine these integrals of motion have been developed and
tested, see Taff and Hall (1977); Gronchi et al. (2010, 2011, 2015, 2017).

In this paper we address the problem of computing the orbit of a celestial
body (hereafter, the OD problem) from one topocentric position
$\mathcal{P}_{1}=(\alpha_{1},\delta_{1},\rho_{1})$ at epoch $t_{1}$ and a VSA
of optical observations, from which we can derive an attributable
$\mathcal{A}_{2}=(\alpha_{2},\delta_{2},\dot{\alpha}_{2},\dot{\delta}_{2})$ at
the mean epoch $\overline{t}_{2}$ of the arc. Here, $\alpha$, $\delta$, $\rho$
denote respectively right ascension, declination and topocentric distance of
the body, while $\dot{\alpha},\dot{\delta}$ stand for the angular rates.

Using the conservation of angular momentum, energy and Laplace-Lenz vector we
write the OD problem as a system of polynomial equations in the unknowns
$\dot{\rho}_{1}$, $\dot{\alpha}_{1}$, $\dot{\delta}_{1}$, $\rho_{2}$,
$\dot{\rho}_{2}$, $z_{2}$, where $z_{2}$ is an auxiliary variable, and the
other variables allow us to obtain an orbit in spherical coordinates at the
two epochs.

We prove that this polynomial system is consistent, namely it generically
admits solutions, at least in the complex field, and we obtain a univariate
polynomial $\mathfrak{v}$ of degree 8 in the unknown range $\rho_{2}$ to solve
the OD problem. Through a computer algebra software we are also able to show
that the degree of $\mathfrak{v}$ is minimum among the degrees of all the
univariate polynomials in $\rho_{2}$ solving this problem.

The proposed method is relevant for different purposes, such as the
computation of a preliminary orbit of an Earth satellite with radar and
optical observations, the detection of maneuvres of an Earth satellite, the
recovery of asteroids which are lost due to a planetary encounter, if the
latter can be modeled as an instantaneous change of direction of the velocity
vector like in Öpik theory (Öpik, 1976). Here we show the results of the
application of our algorithm to the orbits of some near-Earth asteroids, whose
positions have been changed to increase the effect of the close encounter with
the Earth.

The paper is organized as follows: after recalling the expressions of the
Keplerian integrals in Section 2, we introduce the OD problem in Section 3 as
an overdetermined polynomial system, whose consistency is shown in Section 4,
where the univariate polynomial $\mathfrak{v}$ is derived. The minimality of
the degree of $\mathfrak{v}$ is proved in Section 5. In Section 6 we discuss
the selection of the solutions. Finally, in Section 7 we present the numerical
tests.

##  2 Keplerian integrals

Let us consider a celestial body whose dynamics can be modeled by Kepler’s
problem

$\ddot{{\bm{r}}}=-\mu\displaystyle\frac{{\bm{r}}}{|{\bm{r}}|^{3}},$ (1)

where ${\bm{r}}$ and $\mu$ denote respectively its position and the
gravitational parameter. Equation (1) admits the first integrals

${\bm{c}}={\bm{r}}\times\dot{{\bm{r}}},\hskip 17.07182pt{\cal
E}=\frac{1}{2}|\dot{\bm{r}}|^{2}-\frac{\mu}{|{\bm{r}}|},\hskip
17.07182pt{\bm{L}}=\frac{1}{\mu}\dot{\bm{r}}\times{\bm{c}}-\frac{{\bm{r}}}{|{\bm{r}}|},$
(2)

corresponding to the angular momentum, the energy, and the Laplace-Lenz vector
of the body, respectively. We call Keplerian integrals these constants of
motion.

Note that we can write

$\mu{\bm{L}}=\Bigl{(}|\dot{\bm{r}}|^{2}-\frac{\mu}{|{\bm{r}}|}\Bigr{)}{\bm{r}}-(\dot{\bm{r}}\cdot{\bm{r}})\dot{\bm{r}}.$

##  3 The OD problem

Let us introduce a reference frame whose origin is the center of force of the
Keplerian dynamics. We consider the Keplerian motion of a celestial body
around the origin observed from a moving point of view. Assume the position
$\bm{q}$ and the velocity $\dot{\bm{q}}$ of the observer are known functions
of time. The position of the observed body is given by

${\bm{r}}=\bm{q}+\rho\hat{\bm{e}}^{\rho},$ (3)

where $\rho$ represents the topocentric distance, and

$\hat{\bm{e}}^{\rho}=(\cos\delta\cos\alpha,\cos\delta\sin\alpha,\sin\delta)$

is the line of sight.

Let us assume that we know the position vector

${\cal
P}_{1}=(\alpha_{1},\delta_{1},\rho_{1})\in[-\pi,\pi)\times(-\pi/2,\pi/2)\times\mathbb{R}^{+}$

at epoch $t_{1}$, and a very short arc of optical observations of the same
object, from which we compute the attributable

${\cal
A}_{2}=(\alpha_{2},\delta_{2},\dot{\alpha}_{2},\dot{\delta}_{2})\in[-\pi,\pi)\times(-\pi/2,\pi/2)\times\mathbb{R}^{2}$

at the mean epoch $\bar{t}_{2}$. We can combine these data to compute the
quantities $\dot{\rho}_{1}$, $\xi_{1}$, $\zeta_{1}$, $\rho_{2}$,
$\dot{\rho}_{2}$, where

$\xi_{1}=\rho_{1}\dot{\alpha}_{1}\cos{\delta_{1}},\qquad\zeta_{1}=\rho_{1}\dot{\delta}_{1},$

which are missing to have a complete set of orbital elements at both epochs.
For this purpose, we use the conservation of the integrals listed in (2).

Hereafter, we shall use subscripts $1,2$ for all the quantities relative to
epochs $t_{1}$, $\bar{t}_{2}$.

The dependence of position and velocity on the unknowns is given by

$\begin{split}\dot{\bm{r}}_{1}&=\dot{\bm{q}}_{1}+\dot{\rho}_{1}\hat{\bm{e}}^{\rho}_{1}+\xi_{1}\hat{\bm{e}}^{\alpha}_{1}+\zeta_{1}\hat{\bm{e}}^{\delta}_{1},\cr{\bm{r}}_{2}&=\bm{q}_{2}+\rho_{2}\hat{\bm{e}}^{\rho}_{2},\cr\dot{\bm{r}}_{2}&=\dot{\bm{q}}_{2}+\dot{\rho}_{2}\hat{\bm{e}}^{\rho}_{2}+\rho_{2}\hat{{\bm{e}}}^{\perp}_{2},\end{split}$
(4)

where

$\hat{{\bm{e}}}^{\perp}_{2}=\dot{\alpha}_{2}\cos\delta_{2}\hat{\bm{e}}^{\alpha}_{2}+\dot{\delta}_{2}\hat{\bm{e}}^{\delta}_{2}$

is a known vector. Note that
${\bm{r}}_{1}=\bm{q}_{1}+\rho_{1}\hat{\bm{e}}^{\rho}_{1}$ is known.

To write a polynomial system, we replace the term $\frac{\mu}{|{\bm{r}}_{2}|}$
with an auxiliary variable $z_{2}$, independent from $\rho_{2}$, and set

$\displaystyle\mu\widetilde{{\bm{L}}}_{2}$
$\displaystyle=\Bigl{(}|\dot{\bm{r}}_{2}|^{2}-z_{2}\Bigr{)}{\bm{r}}_{2}-(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2},$
$\displaystyle\widetilde{{\cal E}}_{2}$
$\displaystyle=\frac{1}{2}|\dot{\bm{r}}_{2}|^{2}-z_{2}.$

######  Lemma 1.

The following relations hold:

$\mu^{2}|{\bm{L}}_{1}|^{2}=2{\cal
E}_{1}|{\bm{c}}_{1}|^{2}+\mu^{2},\qquad\mu^{2}|\widetilde{{\bm{L}}}_{2}|^{2}=2\widetilde{{\cal
E}}_{2}|{\bm{c}}_{2}|^{2}+z_{2}^{2}|{\bm{r}}_{2}|^{2}.$ (5)

###### Proof.

From the definitions of ${\bm{L}}_{1}$ and $\widetilde{{\bm{L}}}_{2}$ we have

$\displaystyle\mu^{2}|{\bm{L}}_{1}|^{2}$
$\displaystyle=|\dot{\bm{r}}_{1}\times{\bm{c}}_{1}|^{2}-2\frac{\mu}{|{\bm{r}}_{1}|}(\dot{\bm{r}}_{1}\times{\bm{c}}_{1}\cdot{\bm{r}}_{1})+\mu^{2},$
$\displaystyle\mu^{2}|\widetilde{{\bm{L}}}_{2}|^{2}$
$\displaystyle=|\dot{\bm{r}}_{2}\times{\bm{c}}_{2}|^{2}-2z_{2}(\dot{\bm{r}}_{2}\times{\bm{c}}_{2}\cdot{\bm{r}}_{2})+z_{2}^{2}|{\bm{r}}_{2}|^{2}.$

We note that, at both epochs,

$|\dot{\bm{r}}\times{\bm{c}}|^{2}=|\dot{\bm{r}}|^{2}|{\bm{c}}|^{2},\qquad\dot{\bm{r}}\times{\bm{c}}\cdot{\bm{r}}=|{\bm{c}}|^{2},$

so that

$\displaystyle\mu^{2}|{\bm{L}}_{1}|^{2}$
$\displaystyle=|{\bm{c}}_{1}|^{2}\Bigl{(}|\dot{\bm{r}}_{1}|^{2}-2\frac{\mu}{|{\bm{r}}_{1}|}\Bigr{)}+\mu^{2},$
$\displaystyle\mu^{2}|\widetilde{{\bm{L}}}_{2}|^{2}$
$\displaystyle=|{\bm{c}}_{2}|^{2}\bigl{(}|\dot{\bm{r}}_{2}|^{2}-2z_{2}\bigr{)}+z_{2}^{2}|{\bm{r}}_{2}|^{2}.$

Relations (5) immediately follow from the definitions of ${\cal E}_{1}$ and
$\widetilde{{\cal E}}_{2}$. ∎

We consider the overdetermined polynomial system

${\bm{c}}_{1}={\bm{c}}_{2},\quad{\bm{L}}_{1}=\widetilde{{\bm{L}}}_{2},\quad{\cal
E}_{1}=\widetilde{{\cal E}}_{2},\quad z_{2}^{2}|{\bm{r}}_{2}|^{2}=\mu^{2},$
(6)

consisting of 8 equations in the 6 unknowns $\dot{\rho}_{1}$, $\xi_{1}$,
$\zeta_{1}$, $\rho_{2}$, $\dot{\rho}_{2}$, $z_{2}$. System (6) corresponds to
the equations of our OD problem.

##  4 Consistency of the equations

We will show the following property:

######  Theorem 1.

The overdetermined polynomial system (6) is generically consistent, i.e. it
always has solutions in the complex field for a generic choice of the data
${\cal P}_{1}$, ${\cal A}_{2}$, $\bm{q}_{1}$, $\dot{\bm{q}}_{1}$,
$\bm{q}_{2}$, $\dot{\bm{q}}_{2}$.

The proof makes use of the results presented in the following subsections.

We start by noting that, because of relations (5), equation

$z_{2}^{2}|{\bm{r}}_{2}|^{2}=\mu^{2}$

is a consequence of the reduced system

${\bm{c}}_{1}={\bm{c}}_{2},\quad{\bm{L}}_{1}=\widetilde{{\bm{L}}}_{2},\quad{\cal
E}_{1}=\widetilde{{\cal E}}_{2},$ (7)

of 7 equations in 6 unknowns. Therefore, to prove Theorem 1 we show that
system (7) is generically consistent.

Moreover, for a generical choice of the data ${\cal P}_{1}$, ${\cal A}_{2}$,
$\bm{q}_{1}$, $\dot{\bm{q}}_{1}$, $\bm{q}_{2}$, $\dot{\bm{q}}_{2}$, system (7)
is equivalent to

$\mathfrak{q}_{1}=\mathfrak{q}_{2}=\mathfrak{q}_{3}=\mathfrak{q}_{4}=\mathfrak{q}_{5}=\mathfrak{q}_{6}=\mathfrak{q}_{7}=0,$

where

$\begin{split}\mathfrak{q}_{1}&=(\bm{c}_{1}-\bm{c}_{2})\cdot{\bm{W}}_{12},\cr\mathfrak{q}_{2}&=(\bm{c}_{1}-\bm{c}_{2})\cdot\bm{D}_{1}\times{\bm{W}}_{12},\cr\mathfrak{q}_{3}&=(\bm{c}_{1}-\bm{c}_{2})\cdot\bm{D}_{2}\times{\bm{W}}_{12},\cr\mathfrak{q}_{4}&=\mu(\bm{L}_{1}-\tilde{\bm{L}}_{2})\cdot\bm{D}_{1},\cr\mathfrak{q}_{5}&=\mu(\bm{L}_{1}-\tilde{\bm{L}}_{2})\cdot\bm{D}_{2},\cr\mathfrak{q}_{6}&=\mu(\bm{L}_{1}-\tilde{\bm{L}}_{2})\cdot(\bm{r}_{1}\times\hat{\bm{e}}^{\rho}_{2}),\cr\mathfrak{q}_{7}&={\cal
E}_{1}-\tilde{{\cal E}}_{2},\end{split}$ (8)

with

${\bm{D}}_{j}=\bm{q}_{j}\times\hat{\bm{e}}^{\rho}_{j},\quad j=1,2,$ (9)

and

${\bm{W}}_{12}={\bm{D}}_{1}\times{\bm{D}}_{2}.$ (10)

###  4.1 The angular momentum

The angular momenta ${\bm{c}}_{1}$ and ${\bm{c}}_{2}$, written in terms of the
unknowns, become

$\displaystyle{\bm{c}}_{1}$
$\displaystyle={\bm{D}}_{1}\dot{\rho}_{1}+{\bm{N}}_{1}\xi_{1}+{\bm{O}}_{1}\zeta_{1}+{\bm{P}}_{1},$
$\displaystyle{\bm{c}}_{2}$
$\displaystyle={\bm{D}}_{2}\dot{\rho}_{2}+{\bm{E}}_{2}\rho_{2}^{2}+{\bm{F}}_{2}\rho_{2}+{\bm{G}}_{2},$

where ${\bm{D}}_{1}$, ${\bm{D}}_{2}$ are defined as in (9), and

$\displaystyle{\bm{N}}_{1}$
$\displaystyle={\bm{r}}_{1}\times\hat{\bm{e}}^{\alpha}_{1},$
$\displaystyle{\bm{E}}_{2}$
$\displaystyle=\dot{\alpha}_{2}\cos\delta_{2}\hat{\bm{e}}^{\delta}_{2}-\dot{\delta}_{2}\hat{\bm{e}}^{\alpha}_{2},$
$\displaystyle{\bm{O}}_{1}$
$\displaystyle={\bm{r}}_{1}\times\hat{\bm{e}}^{\delta}_{1},$
$\displaystyle{\bm{F}}_{2}$
$\displaystyle=\dot{\alpha}_{2}\cos\delta_{2}(\bm{q}_{2}\times\hat{\bm{e}}^{\alpha}_{2})+\dot{\delta}_{2}(\bm{q}_{2}\times\hat{\bm{e}}^{\delta}_{2})+(\hat{\bm{e}}^{\rho}_{2}\times\dot{\bm{q}}_{2}),$
$\displaystyle{\bm{P}}_{1}$
$\displaystyle={\bm{r}}_{1}\times\dot{\bm{q}}_{1},$
$\displaystyle{\bm{G}}_{2}$ $\displaystyle=\bm{q}_{2}\times\dot{\bm{q}}_{2}.$

Therefore, we can write

${\bm{c}}_{1}-{\bm{c}}_{2}={\bm{D}}_{1}\dot{\rho}_{1}-{\bm{D}}_{2}\dot{\rho}_{2}+{\bm{J}}(\xi_{1},\zeta_{1},\rho_{2}),$
(11)

with

${\bm{J}}(\xi_{1},\zeta_{1},\rho_{2})={\bm{N}}_{1}\xi_{1}+{\bm{O}}_{1}\zeta_{1}-{\bm{E}}_{2}\rho_{2}^{2}-{\bm{F}}_{2}\rho_{2}+{\bm{P}}_{1}-{\bm{G}}_{2}.$
(12)

###  4.2 Elimination of variables

Using relations (11), (12), the polynomials $\mathfrak{q}_{1}$,
$\mathfrak{q}_{2}$, $\mathfrak{q}_{3}$ defined in (8) can be written as

$\displaystyle\mathfrak{q}_{1}$
$\displaystyle=Q^{(1)}_{100}\xi_{1}+Q^{(1)}_{010}\zeta_{1}+Q^{(1)}_{002}\rho_{2}^{2}+Q^{(1)}_{001}\rho_{2}+Q^{(1)}_{000},$
$\displaystyle\mathfrak{q}_{2}$
$\displaystyle=-|{\bm{W}}_{12}|^{2}\dot{\rho}_{2}+Q^{(2)}_{100}\xi_{1}+Q^{(2)}_{010}\zeta_{1}+Q^{(2)}_{002}\rho_{2}^{2}+Q^{(2)}_{001}\rho_{2}+Q^{(2)}_{000},$
$\displaystyle\mathfrak{q}_{3}$
$\displaystyle=|{\bm{W}}_{12}|^{2}\dot{\rho}_{1}+Q^{(3)}_{100}\xi_{1}+Q^{(3)}_{010}\zeta_{1}+Q^{(3)}_{002}\rho_{2}^{2}+Q^{(3)}_{001}\rho_{2}+Q^{(3)}_{000},$

where

$\displaystyle Q^{(1)}_{100}$ $\displaystyle={\bm{N}}_{1}\cdot{\bm{W}}_{12},$
$\displaystyle Q^{(1)}_{010}$ $\displaystyle={\bm{O}}_{1}\cdot{\bm{W}}_{12},$
$\displaystyle Q^{(1)}_{002}$ $\displaystyle=-{\bm{E}}_{2}\cdot{\bm{W}}_{12},$
$\displaystyle Q^{(1)}_{001}$ $\displaystyle=-{\bm{F}}_{2}\cdot{\bm{W}}_{12},$
$\displaystyle Q^{(1)}_{000}$
$\displaystyle=({\bm{G}}_{2}-{\bm{P}}_{1})\cdot{\bm{W}}_{12},$ $\displaystyle
Q^{(2)}_{100}$
$\displaystyle={\bm{N}}_{1}\cdot{\bm{D}}_{1}\times{\bm{W}}_{12},$
$\displaystyle Q^{(2)}_{010}$
$\displaystyle={\bm{O}}_{1}\cdot{\bm{D}}_{1}\times{\bm{W}}_{12},$
$\displaystyle Q^{(2)}_{002}$
$\displaystyle=-{\bm{E}}_{2}\cdot{\bm{D}}_{1}\times{\bm{W}}_{12},$
$\displaystyle Q^{(2)}_{001}$
$\displaystyle=-{\bm{F}}_{2}\cdot{\bm{D}}_{1}\times{\bm{W}}_{12},$
$\displaystyle Q^{(2)}_{000}$
$\displaystyle=({\bm{P}}_{1}-{\bm{G}}_{2})\cdot{\bm{D}}_{1}\times{\bm{W}}_{12},$
$\displaystyle Q^{(3)}_{100}$
$\displaystyle={\bm{N}}_{1}\cdot{\bm{D}}_{2}\times{\bm{W}}_{12},$
$\displaystyle Q^{(3)}_{010}$
$\displaystyle={\bm{O}}_{1}\cdot{\bm{D}}_{2}\times{\bm{W}}_{12},$
$\displaystyle Q^{(3)}_{002}$
$\displaystyle=-{\bm{E}}_{2}\cdot{\bm{D}}_{2}\times{\bm{W}}_{12},$
$\displaystyle Q^{(3)}_{001}$
$\displaystyle=-{\bm{F}}_{2}\cdot{\bm{D}}_{2}\times{\bm{W}}_{12},$
$\displaystyle Q^{(3)}_{000}$
$\displaystyle=({\bm{P}}_{1}-{\bm{G}}_{2})\cdot{\bm{D}}_{2}\times{\bm{W}}_{12}.$

In particular, the occurence of the variables $\dot{\rho}_{1}$, $\xi_{1}$,
$\zeta_{1}$, $\dot{\rho}_{2}$ in $\mathfrak{q}_{1}$,
$\mathfrak{q}_{2}$,$\mathfrak{q}_{3}$ is at most linear. From the definitions
above, we obtain

$\displaystyle Q^{(1)}_{100}$
$\displaystyle=({\bm{N}}_{1}\times{\bm{D}}_{1})\cdot{\bm{D}}_{2}=\left[({\bm{r}}_{1}\times\hat{\bm{e}}^{\alpha}_{1})\times({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})\right]\cdot{\bm{D}}_{2}$
$\displaystyle=\left[({\bm{r}}_{1}\times\hat{\bm{e}}^{\alpha}_{1})\cdot\hat{\bm{e}}^{\rho}_{1}\right]{\bm{r}}_{1}\cdot{\bm{D}}_{2}=-({\bm{r}}_{1}\cdot\hat{\bm{e}}^{\delta}_{1})({\bm{r}}_{1}\cdot{\bm{D}}_{2}),$
$\displaystyle Q^{(1)}_{010}$
$\displaystyle=({\bm{O}}_{1}\times{\bm{D}}_{1})\cdot{\bm{D}}_{2}=\left[({\bm{r}}_{1}\times\hat{\bm{e}}^{\delta}_{1})\times({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})\right]\cdot{\bm{D}}_{2}$
$\displaystyle=\left[({\bm{r}}_{1}\times\hat{\bm{e}}^{\delta}_{1})\cdot\hat{\bm{e}}^{\rho}_{1}\right]{\bm{r}}_{1}\cdot{\bm{D}}_{2}=({\bm{r}}_{1}\cdot\hat{\bm{e}}^{\alpha}_{1})({\bm{r}}_{1}\cdot{\bm{D}}_{2}).$

Therefore, the coefficients $Q_{100}^{(1)}$ and $Q_{010}^{(1)}$ are both
vanishing iff

${\bm{r}}_{1}\cdot{\bm{D}}_{2}=0\qquad\mbox{or}\qquad{\bm{D}}_{1}=\bm{0}.$

In fact, the second condition in the alternative above is equivalent to

${\bm{r}}_{1}\cdot\hat{\bm{e}}^{\alpha}_{1}={\bm{r}}_{1}\cdot\hat{\bm{e}}^{\delta}_{1}=0.$

Thus, for a generic choice of the data, the angular momentum equations allow
us to eliminate $\dot{\rho}_{1}$, $\dot{\rho}_{2}$ and one variable between
$\xi_{1}$ and $\zeta_{1}$. Assuming $Q^{(1)}_{100}\neq 0$, we choose to
eliminate $\xi_{1}$, which can be written as

$\xi_{1}=-\frac{Q^{(1)}_{010}\zeta_{1}+\sum_{h=0}^{2}Q^{(1)}_{00h}\rho_{2}^{h}}{Q^{(1)}_{100}}.$
(13)

Substituting (13) in equations $\mathfrak{q}_{3}=\mathfrak{q}_{2}=0$ and
assuming ${\bm{W}}_{12}\neq\bm{0}$ we find

$\displaystyle\dot{\rho}_{1}$
$\displaystyle=-\frac{1}{|{\bm{W}}_{12}|^{2}Q^{(1)}_{100}}\Bigl{[}\left(Q^{(3)}_{010}Q^{(1)}_{100}-Q^{(1)}_{010}Q^{(3)}_{100}\right)\zeta_{1}+\sum_{h=0}^{2}\left(Q^{(3)}_{00h}Q^{(1)}_{100}-Q^{(1)}_{00h}Q^{(3)}_{100}\right)\rho_{2}^{h}\Bigr{]}$
(14) $\displaystyle\dot{\rho}_{2}$
$\displaystyle=\frac{1}{|{\bm{W}}_{12}|^{2}Q^{(1)}_{100}}\Bigl{[}\left(Q^{(2)}_{010}Q^{(1)}_{100}-Q^{(1)}_{010}Q^{(2)}_{100}\right)\zeta_{1}+\sum_{h=0}^{2}\left(Q^{(2)}_{00h}Q^{(1)}_{100}-Q^{(1)}_{00h}Q^{(2)}_{100}\right)\rho_{2}^{h}\Bigr{]}.$
(15)

Hence, by using relations (13), (14), and (15) we can eliminate the variables
$\xi_{1}$, $\dot{\rho}_{1}$, $\dot{\rho}_{2}$ in the generators
$\mathfrak{q}_{4}$, $\mathfrak{q}_{5}$, $\mathfrak{q}_{6}$,
$\mathfrak{q}_{7}$. The resulting polynomials, named
$\tilde{\mathfrak{q}}_{4}$, $\tilde{\mathfrak{q}}_{5}$,
$\tilde{\mathfrak{q}}_{6}$, $\tilde{\mathfrak{q}}_{7}$, can be written as

$\displaystyle\tilde{\mathfrak{q}}_{4}$
$\displaystyle=({\bm{r}}_{2}\cdot{\bm{D}}_{1})z_{2}+h_{4}(\zeta_{1},\rho_{2}),$
$\displaystyle\tilde{\mathfrak{q}}_{5}$
$\displaystyle=\tilde{\mathfrak{q}}_{5}(\zeta_{1},\rho_{2}),$
$\displaystyle\tilde{\mathfrak{q}}_{6}$
$\displaystyle=-({\bm{r}}_{1}\cdot{\bm{D}}_{2})z_{2}+h_{6}(\zeta_{1},\rho_{2}),$
$\displaystyle\tilde{\mathfrak{q}}_{7}$
$\displaystyle=z_{2}+h_{7}(\zeta_{1},\rho_{2}),$

for some bivariate polynomials $h_{4}$, $h_{6}$, $h_{7}$. For later use we
observe that

$\tilde{\mathfrak{q}}_{5}=P^{(5)}_{12}\zeta_{1}\rho_{2}^{2}+P^{(5)}_{11}\zeta_{1}\rho_{2}+P^{(5)}_{10}\zeta_{1}+P^{(5)}_{04}\rho_{2}^{4}+P^{(5)}_{03}\rho_{2}^{3}+P^{(5)}_{02}\rho_{2}^{2}+P^{(5)}_{01}\rho_{2}+P^{(5)}_{00},$

for some coefficients $P^{(5)}_{ij}$. Moreover, using equation
$\tilde{\mathfrak{q}}_{7}=0$, we can eliminate $z_{2}$ from
$\tilde{\mathfrak{q}}_{4}$ and $\tilde{\mathfrak{q}}_{6}$. In particular, with
this elimination $\tilde{\mathfrak{q}}_{6}$ becomes

$\mathfrak{p}_{6}=P^{(6)}_{20}\zeta_{1}^{2}+P^{(6)}_{12}\zeta_{1}\rho_{2}^{2}+P^{(6)}_{11}\zeta_{1}\rho_{2}+P^{(6)}_{10}\zeta_{1}+P^{(6)}_{04}\rho_{2}^{4}+P^{(6)}_{03}\rho_{2}^{3}+P^{(6)}_{02}\rho_{2}^{2}+P^{(6)}_{01}\rho_{2}+P^{(6)}_{00},$

for some coefficients $P^{(6)}_{ij}$.

###  4.3 A redundant relation

Here, we prove that equation $\mathfrak{q}_{4}=0$ can be dropped from system
(7) without losing any solution. In particular, we show:

######  Proposition 1.

If ${\bm{r}}_{1}\cdot{\bm{D}}_{2}\neq 0$, the polynomial $\mathfrak{q}_{4}$ is
generated by

$\mathfrak{q}_{1},\,\mathfrak{q}_{2},\,\mathfrak{q}_{3},\,\mathfrak{q}_{6}.$

###### Proof.

Let us consider the polynomials $\tilde{\mathfrak{q}}_{4}$,
$\tilde{\mathfrak{q}}_{6}$ obtained from $\mathfrak{q}_{4}$,
$\mathfrak{q}_{6}$ by eliminating $\dot{\rho}_{1}$, $\dot{\rho}_{2}$,
$\xi_{1}$ through relations
$\mathfrak{q}_{1}=\mathfrak{q}_{2}=\mathfrak{q}_{3}=0$. We note that

$\tilde{\mathfrak{q}}_{4}=\sum_{j=1}^{3}\mathfrak{a}_{j}\mathfrak{q}_{j}+\mathfrak{q}_{4},\qquad\tilde{\mathfrak{q}}_{6}=\sum_{j=1}^{3}\mathfrak{b}_{j}\mathfrak{q}_{j}+\mathfrak{q}_{6},$

for some polynomials $\mathfrak{a}_{j}$, $\mathfrak{b}_{j}$. We prove that
$\tilde{\mathfrak{q}}_{6}$ divides $\tilde{\mathfrak{q}}_{4}$; in particular,

$\tilde{\mathfrak{q}}_{4}=-\frac{{\bm{r}}_{2}\cdot{\bm{D}}_{1}}{{\bm{r}}_{1}\cdot{\bm{D}}_{2}}\,\tilde{\mathfrak{q}}_{6}.$
(16)

For this purpose, we first note that system
$\mathfrak{q}_{1}=\mathfrak{q}_{2}=\mathfrak{q}_{3}=0$ is equivalent to
${\bm{c}}_{1}={\bm{c}}_{2}$ and, inserting the last relation in the
expressions of $\mathfrak{q}_{4}$, $\mathfrak{q}_{6}$, after eliminating
$\dot{\rho}_{1}$, $\dot{\rho}_{2}$, $\xi_{1}$, we obtain

$\displaystyle\tilde{\mathfrak{q}}_{4}$
$\displaystyle=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})\left[{\bm{r}}_{1}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]+z_{2}({\bm{r}}_{2}\cdot{\bm{D}}_{1}),$
$\displaystyle\tilde{\mathfrak{q}}_{6}$
$\displaystyle=({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})\left[{\bm{r}}_{1}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]-z_{2}({\bm{r}}_{1}\cdot{\bm{D}}_{2}),$

where ${\bm{c}}_{1}$, ${\bm{c}}_{2}$, ${\bm{r}}_{2}$, $\dot{\bm{r}}_{1}$,
$\dot{\bm{r}}_{2}$ are meant as functions of $\zeta_{1}$, $\rho_{2}$ only.
More details about these computations can be found in Appendix A.

We have

${\bm{r}}_{2}\cdot{\bm{D}}_{1}=-({\bm{r}}_{1}\times{\bm{r}}_{2})\cdot\hat{\bm{e}}^{\rho}_{1},\qquad{\bm{r}}_{1}\cdot{\bm{D}}_{2}=({\bm{r}}_{1}\times{\bm{r}}_{2})\cdot\hat{\bm{e}}^{\rho}_{2}.$

Moreover, ${\bm{c}}_{1}={\bm{c}}_{2}$ implies that
${\bm{r}}_{1}\times{\bm{r}}_{2}$ is parallel to ${\bm{c}}_{1}$ and
${\bm{c}}_{2}$. Hence

$({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})\left[({\bm{r}}_{1}\times{\bm{r}}_{2})\cdot\hat{\bm{e}}^{\rho}_{1}\right]=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})\left[({\bm{r}}_{1}\times{\bm{r}}_{2})\cdot\hat{\bm{e}}^{\rho}_{2}\right],$

that is

$-({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})({\bm{r}}_{2}\cdot{\bm{D}}_{1})=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})({\bm{r}}_{1}\cdot{\bm{D}}_{2}).$

The latter relation immediately yields (16).

Setting

$A=-\frac{{\bm{r}}_{2}\cdot{\bm{D}}_{1}}{{\bm{r}}_{1}\cdot{\bm{D}}_{2}},$

we can write

$\mathfrak{q}_{4}=\tilde{\mathfrak{q}}_{4}-\sum_{j=1}^{3}\mathfrak{a}_{j}\mathfrak{q}_{j}=A\tilde{\mathfrak{q}}_{6}-\sum_{j=1}^{3}\mathfrak{a}_{j}\mathfrak{q}_{j}=\sum_{j=1}^{3}(A\mathfrak{b}_{j}-\mathfrak{a}_{j})\mathfrak{q}_{j}+A\mathfrak{q}_{6},$

which concludes the proof. ∎

###  4.4 The univariate polynomial

We compute a univariate polynomial $\mathfrak{u}(\rho_{2})$ of degree 8, which
is a consequence of the equations in (7). Let us consider the polynomial
system

$\mathfrak{q}_{1}=\mathfrak{q}_{2}=\mathfrak{q}_{3}=\mathfrak{q}_{5}=\mathfrak{q}_{6}=\mathfrak{q}_{7}=0$
(17)

including all the polynomials in (8) except $\mathfrak{q}_{4}$. To solve the
OD problem, we solve the polynomial system defined by (17).

Now, we compute the resultant of $\tilde{\mathfrak{q}}_{5}$ and
$\mathfrak{p}_{6}$ with respect to $\zeta_{1}$, which is denoted by
$\mathfrak{v}(\rho_{2})$. For this purpose, the terms in
$\tilde{\mathfrak{q}}_{5}$ and $\mathfrak{p}_{6}$ are grouped in the following
way:

$\displaystyle\tilde{\mathfrak{q}}_{5}$
$\displaystyle=a_{1}(\rho_{2})\zeta_{1}+a_{0}(\rho_{2}),$
$\displaystyle\mathfrak{p}_{6}$
$\displaystyle=P^{(6)}_{20}\zeta_{1}^{2}+b_{1}(\rho_{2})\zeta_{1}+b_{0}(\rho_{2}),$

where

$\displaystyle a_{1}(\rho_{2})$
$\displaystyle=P^{(5)}_{12}\rho_{2}^{2}+P^{(5)}_{11}\rho_{2}+P^{(5)}_{10},$
$\displaystyle a_{0}(\rho_{2})$
$\displaystyle=P^{(5)}_{04}\rho_{2}^{4}+P^{(5)}_{03}\rho_{2}^{3}+P^{(5)}_{02}\rho_{2}^{2}+P^{(5)}_{01}\rho_{2}+P^{(5)}_{00},$
$\displaystyle b_{1}(\rho_{2})$
$\displaystyle=P^{(6)}_{12}\rho_{2}^{2}+P^{(6)}_{11}\rho_{2}+P^{(6)}_{10},$
$\displaystyle b_{0}(\rho_{2})$
$\displaystyle=P^{(6)}_{04}\rho_{2}^{4}+P^{(6)}_{03}\rho_{2}^{3}+P^{(6)}_{02}\rho_{2}^{2}+P^{(6)}_{01}\rho_{2}+P^{(6)}_{00}.$

Therefore, we have

$\mathfrak{v}(\rho_{2})=a_{1}(\rho_{2})a_{0}(\rho_{2})b_{1}(\rho_{2})-a_{0}^{2}(\rho_{2})P^{(6)}_{20}-b_{0}(\rho_{2})a_{1}^{2}(\rho_{2}),$

which has degree 8.

We can conclude the proof of Theorem 1. In fact, $\mathfrak{v}(\rho_{2})$ has
(generically) 8 complex roots. For each of these roots, equation
$\tilde{\mathfrak{q}}_{5}=0$ gives a unique value of $\zeta_{1}$. Then,
equation $\tilde{\mathfrak{q}}_{7}=0$ and relations (13), (14), (15),
respectively, allow us to compute unique values for the remaining components
$z_{2}$, $\xi_{1}$, $\dot{\rho}_{1}$, $\dot{\rho}_{2}$ of the solutions.

##  5 An optimal property

We show that the univariate polynomial $\mathfrak{v}(\rho_{2})$ has the
minimal degree among all the univariate polynomials in the variable $\rho_{2}$
that are algebraic consequences of

${\bm{c}}_{1}-{\bm{c}}_{2},\ \mu({\bm{L}}_{1}-\widetilde{{\bm{L}}}_{2}),\
{\cal E}_{1}-\widetilde{{\cal E}}_{2}.$

More precisely, let us introduce the polynomial ideal

$I=\langle{\bm{c}}_{1}-{\bm{c}}_{2},\
\mu({\bm{L}}_{1}-\widetilde{{\bm{L}}}_{2}),\ {\cal E}_{1}-\widetilde{{\cal
E}}_{2}\rangle\subseteq\mathbb{R}[\dot{\rho}_{1},\xi_{1},\zeta_{1},\rho_{2},\dot{\rho}_{2},z_{2}].$

We will show the following result:

######  Theorem 2.

For a generic choice of the data we can find a Gröbner basis $\mathcal{G}$ in
$\mathbb{R}[\dot{\rho}_{1},\xi_{1},\zeta_{1},\rho_{2},\dot{\rho}_{2},z_{2}]$
of the ideal $I$ for the lexicographic order

$\dot{\rho}_{1}\succ\dot{\rho}_{2}\succ\xi_{1}\succ
z_{2}\succ\zeta_{1}\succ\rho_{2}.$

such that $\mathfrak{v}\in\mathcal{G}$.

###### Proof.

As a consequence of Proposition 1 we obtain

$I=\langle\mathfrak{q}_{1},\mathfrak{q}_{2},\mathfrak{q}_{3},\mathfrak{q}_{5},\mathfrak{q}_{6},\mathfrak{q}_{7}\rangle,$

where we dropped $\mathfrak{q}_{4}$ from the set of generators. Let us
consider the elimination ideal

$I^{\prime}=I\cap\mathbb{R}[\zeta_{1},\rho_{2}]=\langle\tilde{\mathfrak{q}}_{5},\mathfrak{p}_{6}\rangle.$

We compute a Gröbner basis $\mathcal{G}^{\prime}$ of the ideal $I^{\prime}$
for the lexicographic order $\zeta_{1}\succ\rho_{2}$ using the software
_Mathematica_ 111Wolfram Research, Inc. Mathematica, Version 12.1,
https://www.wolfram.com/mathematica.: in this basis we have a univariate
polynomial

$\mathfrak{u}(\rho_{2})=\sum_{i=0}^{8}u_{i}\rho_{2}^{i},$

for some constant coefficients $u_{i}$. The roots of $\mathfrak{u}$ are the
only values of the $\rho_{2}$ components of the solutions of system (17).

We note that, since the resultant $\mathfrak{v}$ is a univariate polynomial of
degree 8, $\mathfrak{v}$ is necessarily proportional to $\mathfrak{u}$, i.e.

$\mathfrak{v}=\kappa\mathfrak{u},\qquad\kappa\in\mathbb{R}\backslash\\{0\\}.$

In place of the generators $\mathfrak{q}_{1}$, $\mathfrak{q}_{2}$,
$\mathfrak{q}_{3}$, we can consider

$\displaystyle\tilde{\mathfrak{q}}_{1}$
$\displaystyle=\xi_{1}+h_{1}(\zeta_{1},\rho_{2}),$
$\displaystyle\tilde{\mathfrak{q}}_{2}$
$\displaystyle=\dot{\rho}_{2}+h_{2}(\zeta_{1},\rho_{2}),$
$\displaystyle\tilde{\mathfrak{q}}_{3}$
$\displaystyle=\dot{\rho}_{1}+h_{3}(\zeta_{1},\rho_{2}),$

where the polynomials $h_{1}$, $h_{2}$, $h_{3}$ are defined by relations (13),
(14), (15).

The set

$\mathcal{G}=\\{\tilde{\mathfrak{q}}_{1},\tilde{\mathfrak{q}}_{2},\tilde{\mathfrak{q}}_{3},\tilde{\mathfrak{q}}_{7}\\}\cup\mathcal{G}^{\prime}$

is a Gröbner basis of the ideal $I$. In fact, the leading term of every
polynomial in $I$ is divisible by the leading term of one polynomial in
$\mathcal{G}$.

∎

As a consequence of the previous theorem we obtain the following property:

######  Corollary 1.

The polynomial $\mathfrak{v}$ has the minimum degree among the univariate
polynomials in $\rho_{2}$ belonging to $I$.

##  6 Selecting the solutions

From the solutions of (7), we first discard the ones with non-real components
and the ones with $\rho_{2},z_{2}\leq 0$. From the remaining solutions we can
compute Keplerian orbital elements at epochs

$\tilde{t}_{1}=t_{1}-\frac{\rho_{1}}{c},\qquad\tilde{t}_{2}^{\,(j)}=\bar{t}_{2}-\frac{\rho_{2}^{(j)}}{c},$

where $c$ is the speed of light and $\rho_{2}^{(j)}$ is the value of
$\rho_{2}$ for the $j$-th solution.

We emphasize that the accepted solutions of the system are such that the
Keplerian elements at the two epochs are the same. In fact, from
${\bm{c}}_{1}-{\bm{c}}_{2}=\bm{0}$ we get

$i_{1}=i_{2},\qquad\Omega_{1}=\Omega_{2},\qquad
a_{1}(1-e_{1}^{2})=a_{2}(1-e_{2}^{2}),$

where $i$, $\Omega$, $a$ and $e$ denote the inclination, ascending node, semi-
major axis and eccentricity. In particular, from
$z_{2}^{2}|{\bm{r}}_{2}|^{2}-\mu^{2}=0$, which is a consequence of system (7),
we get $z_{2}=\pm\frac{\mu}{|{\bm{r}}_{2}|}$. Having discarded the solutions
with $z_{2}\leq 0$, we obtain

$\mu({\bm{L}}_{1}-\widetilde{{\bm{L}}}_{2})=\bm{0}\quad\Leftrightarrow\quad\mu({\bm{L}}_{1}-{\bm{L}}_{2})=\bm{0},$

which implies

$e_{1}=e_{2},\qquad\omega_{1}=\omega_{2}.$

Since the eccentrities $e_{1}$, $e_{2}$ are equal, then also

$a_{1}=a_{2}.$

We observe that for the other orbit determination methods using the Keplerian
integrals (Gronchi et al., 2010, 2011, 2015) the trajectories of the solutions
at the two epochs are not necessarily the same.

In the following subsections we explain how to choose or discard the remaining
solutions.

###  6.1 Selection without covariance

In case of multiple solutions, labeled with $(j)$, we make our choice
according to the following procedure. We propagate each of the computed
orbits, referring to epoch $\tilde{t}_{2}^{\,(j)}$, backward to epoch
$\tilde{t}_{1}$ in the framework of Kepler’s dynamics, and consider the norm
of the differences $\bm{\rho}_{1}^{(j)}-{\cal P}_{1}$, where

$\bm{\rho}_{1}^{(j)}={\bm{r}}^{(j)}(\tilde{t}_{1})-\bm{q}(t_{1})$

is the topocentric position of the $j$-th solution propagated to epoch
$\tilde{t}_{1}$. Note that in general we have

$|\bm{\rho}_{1}^{(j)}|\neq\rho_{1}.$

We select the solution attaining the minimum value

$\min_{j}|\bm{\rho}_{1}^{(j)}-{\cal P}_{1}|.$

###  6.2 Selection using covariance

Assume that the data

${\cal D}=({\cal P}_{1},{\cal A}_{2}),$

where

${\cal P}_{1}=(\alpha_{1},\delta_{1},\rho_{1}),\qquad{\cal
A}_{2}=(\alpha_{2},\delta_{2},\dot{\alpha}_{2},\dot{\delta}_{2}),$

have a covariance matrix

$\Gamma_{{\cal D}}=\left[\begin{array}[]{cc}\Gamma_{{\cal P}_{1}}&0\cr
0&\Gamma_{{\cal A}_{2}}\cr\end{array}\right].$

Let

${\cal S}={\cal S}({\cal D})=({\cal V}_{1}({\cal D}),{\cal R}_{2}({\cal D})),$

with

${\cal V}_{1}=(\xi_{1},\zeta_{1},\dot{\rho}_{1}),\qquad{\cal
R}_{2}=(\rho_{2},\dot{\rho}_{2})$

be a solution of

$\bm{\Phi}({\cal S};{\cal
D})=\left(\begin{array}[]{c}{\bm{c}}_{1}-{\bm{c}}_{2}\cr\left[\mu{\bm{L}}_{1}+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2}\right]\cdot\bm{q}_{2}\times{\bm{r}}_{2}\cr\left[-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\dot{\bm{r}}_{1}-(\frac{1}{2}|\dot{\bm{r}}_{2}|^{2}+{\cal
E}_{1}){\bm{r}}_{2}+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2}\right]\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\end{array}\right)=\bm{0}.$
(18)

Note that, generically, the solutions ${\cal S}$ of system (18) correspond to
the ${\cal S}$ components of the solutions of system (17).

Let $\mathbf{E}_{car}=(\mathbf{E}_{car}^{(1)},\mathbf{E}_{car}^{(2)})$ and
$\mathbf{E}_{att}=(\mathbf{E}_{att}^{(1)},\mathbf{E}_{att}^{(2)})$ be the
vectors of the Cartesian coordinates and the attributable elements at epochs
$\tilde{t}_{1}$ and $\tilde{t}_{2}$.222for simplicity, in this section we drop
the label $(j)$, referring to the possible multiple solutions. Let us
introduce the transformation
$\mathcal{T}_{att}^{car}:\mathbf{E}_{att}\to\mathbf{E}_{car}$ by (3), (4) for
both epochs, and consider the map $\bm{\Psi}$ defined by
$\bm{\Phi}=\bm{\Psi}\circ\mathcal{T}_{att}^{car}$. Equation
$\bm{\Phi}=\mathbf{0}$ is equivalent to the system

${\bm{c}}_{1}-{\bm{c}}_{2}=\bm{0},\quad\mu({\bm{L}}_{1}-\tilde{{\bm{L}}}_{2})\cdot{\bm{D}}_{2}=0,\quad\mu({\bm{L}}_{1}-\tilde{{\bm{L}}}_{2})\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})=0.$

We introduce the vector

$\widetilde{\cal V}_{1}=(\dot{\alpha}_{1},\dot{\delta}_{1},\dot{\rho}_{1}).$

The covariance matrix of the Cartesian coordinates at epoch $\tilde{t}_{1}$ is

$\Gamma_{car}^{(1)}=\frac{\partial\mathbf{E}_{car}^{(1)}}{\partial{{\cal
D}}}\Gamma_{{\cal
D}}\left[\frac{\partial\mathbf{E}_{car}^{(1)}}{\partial{{\cal
D}}}\right]^{T}\,,$

with

$\frac{\partial\mathbf{E}_{car}^{(1)}}{\partial{{\cal
D}}}=\frac{\partial\mathbf{E}_{car}^{(1)}}{\partial\mathbf{E}_{att}^{(1)}}\frac{\partial\mathbf{E}_{att}^{(1)}}{\partial{{\cal
D}}},\qquad\frac{\partial\mathbf{E}_{att}^{(1)}}{\partial{{\cal
D}}}=M\left[\begin{array}[]{c}\begin{array}[]{cc}I_{3}&O_{3\times
4}\cr\end{array}\cr\displaystyle\frac{\partial\widetilde{\cal
V}_{1}}{\partial{{\cal D}}}\cr\end{array}\right],$

where

$M=\left[\begin{array}[]{ccccccc}1&0&0&0&0&0\cr 0&1&0&0&0&0\cr 0&0&0&1&0&0\cr
0&0&0&0&1&0\cr 0&0&1&0&0&0\cr 0&0&0&0&0&1\cr\end{array}\right]$

is the matrix that exchanges the 3-rd, 4-th and 5-th lines, and

$\frac{\partial\widetilde{\cal V}_{1}}{\partial{{\cal
D}}}=\frac{\partial\widetilde{\cal V}_{1}}{\partial{\cal
V}_{1}}\frac{\partial{\cal V}_{1}}{\partial{{\cal
D}}},\qquad\frac{\partial\widetilde{\cal V}_{1}}{\partial{\cal
V}_{1}}=\left[\begin{array}[]{ccc}\frac{1}{\rho_{1}\cos\delta_{1}}&0&0\cr
0&\frac{1}{\rho_{1}}&0\cr 0&0&1\cr\end{array}\right].$

From the implicit function theorem, we have

$\frac{\partial{\widetilde{\cal S}}}{\partial{{\cal D}}}({\cal
D})=-\left[\frac{\partial\bm{\Phi}}{\partial{\widetilde{\cal
S}}}(\mathbf{E}_{att})\right]^{-1}\frac{\partial\bm{\Phi}}{\partial{{\cal
D}}}(\mathbf{E}_{att}),$ $\frac{\partial\bm{\Phi}}{\partial{\widetilde{\cal
S}}}=\left(\frac{\partial\bm{\Psi}}{\partial\mathbf{E}_{car}}\circ{\cal
T}_{att}^{car}\right)\frac{\partial{\cal
T}_{att}^{car}}{\partial{\widetilde{\cal
S}}},\qquad\frac{\partial\bm{\Phi}}{\partial{{\cal
D}}}=\left(\frac{\partial\bm{\Psi}}{\partial\mathbf{E}_{car}}\circ{\cal
T}_{att}^{car}\right)\frac{\partial{\cal T}_{att}^{car}}{\partial{{\cal D}}},$

with

$\widetilde{\cal S}=(\widetilde{\cal V}_{1},{\cal
R}_{2}),\qquad\frac{\partial\widetilde{\cal S}}{\partial{\cal
S}}=\left[\begin{array}[]{cc}\frac{\partial\widetilde{\cal
V}_{1}}{\partial{\cal V}_{1}}&O_{3\times 2}\cr O_{2\times
3}&I_{{2}}\cr\end{array}\right].$

The matrices ${\partial{\cal T}_{att}^{car}}/{\partial{\widetilde{\cal S}}}$
and ${\partial{\cal T}_{att}^{car}}/{\partial{{\cal D}}}$ are respectively
made by columns 3, 4, 6, 11, 12 and by columns 1, 2, 5, 7, 8, 9, 10 of
${\partial\mathbf{E}_{car}}/{\partial\mathbf{E}_{att}}$.

The covariance matrix of the Cartesian coordinates at epoch $\tilde{t}_{2}$ is
given by

$\Gamma_{car}^{(2)}=\frac{\partial\mathbf{E}_{car}^{(2)}}{\partial\mathcal{D}}\Gamma_{\mathcal{D}}\left[\frac{\partial\mathbf{E}_{car}^{(2)}}{\partial\mathcal{D}}\right]^{T},$

with

$\frac{\partial\mathbf{E}_{car}^{(2)}}{\partial\mathcal{D}}=\frac{\partial\mathbf{E}_{car}^{(2)}}{\partial\mathbf{E}_{att}^{(2)}}\frac{\partial\mathbf{E}_{att}^{(2)}}{\partial\mathcal{D}},\qquad\frac{\partial\mathbf{E}_{att}^{(2)}}{\partial\mathcal{D}}=\left[\begin{array}[]{c}\begin{array}[]{cc}O_{4\times
3}&I_{4}\cr\end{array}\cr\frac{\partial\mathcal{R}_{2}}{\partial\mathcal{D}}\end{array}\right]$

and

$\frac{\partial(\rho_{2},\dot{\rho}_{2})}{\partial\widetilde{\cal
S}}=\frac{\partial(\rho_{2},\dot{\rho}_{2})}{\partial\widetilde{\cal
S}}\frac{\partial\widetilde{\cal
S}}{\partial\mathcal{D}},\qquad\frac{\partial(\rho_{2},\dot{\rho}_{2})}{\partial\widetilde{\cal
S}}=\left[\begin{array}[]{ccccc}0&0&0&1&0\cr 0&0&0&0&1\cr\end{array}\right].$

For a given vector ${\bf u}\in\mathbb{R}^{3}$ we define the map

$\mathbb{R}^{3}\ni(u_{1},u_{2},u_{3})={\bf u}\mapsto\widehat{\bf
u}=\left[\begin{array}[]{ccc}0&-u_{3}&u_{2}\cr u_{3}&0&-u_{1}\cr-
u_{2}&u_{1}&0\cr\end{array}\right]\in so(3).$

Then, using $\hat{\bf u}^{T}=-\hat{\bf u}$, we have

$\frac{\partial\bm{\Psi}}{\partial\mathbf{E}_{car}}=\left[\begin{array}[]{cccc}-\widehat{{\dot{\bm{r}}}_{1}}&\widehat{{{\bm{r}}}_{1}}&\widehat{{\dot{\bm{r}}}_{2}}&-\widehat{{{\bm{r}}}_{2}}\cr\stackrel{{\scriptstyle}}{{\displaystyle\frac{\partial\Phi_{4}}{\partial{\bm{r}}_{1}}}}&\displaystyle\frac{\partial\Phi_{4}}{\partial\dot{\bm{r}}_{1}}&\displaystyle\frac{\partial\Phi_{4}}{\partial{\bm{r}}_{2}}&\displaystyle\frac{\partial\Phi_{4}}{\partial\dot{\bm{r}}_{2}}\cr\stackrel{{\scriptstyle}}{{\displaystyle\frac{\partial\Phi_{5}}{\partial{\bm{r}}_{1}}}}&\displaystyle\frac{\partial\Phi_{5}}{\partial\dot{\bm{r}}_{1}}&\displaystyle\frac{\partial\Phi_{5}}{\partial{\bm{r}}_{2}}&\displaystyle\frac{\partial\Phi_{5}}{\partial\dot{\bm{r}}_{2}}\cr\end{array}\right],$

where

$\displaystyle\Phi_{4}$
$\displaystyle=\left[\mu{\bm{L}}_{1}+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2}\right]\cdot\bm{q}_{2}\times{\bm{r}}_{2},$
$\displaystyle\Phi_{5}$
$\displaystyle=\left[-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\dot{\bm{r}}_{1}-(\frac{1}{2}|\dot{\bm{r}}_{2}|^{2}+{\cal
E}_{1}){\bm{r}}_{2}+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2}\right]\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2}),$

and

$\displaystyle\frac{\partial\Phi_{4}}{\partial{\bm{r}}_{1}}=\left(|\dot{\bm{r}}_{1}|^{2}-\frac{\mu}{|{\bm{r}}_{1}|}\right)(\bm{q}_{2}\times{\bm{r}}_{2})-\dot{\bm{r}}_{1}\left[\dot{\bm{r}}_{1}\cdot(\bm{q}_{2}\times{\bm{r}}_{2})\right]+\frac{\mu}{|{\bm{r}}_{1}|^{3}}{\bm{r}}_{1}\left[{\bm{r}}_{1}\cdot(\bm{q}_{2}\times{\bm{r}}_{2})\right],$
$\displaystyle\frac{\partial\Phi_{4}}{\partial\dot{\bm{r}}_{1}}=2\dot{\bm{r}}_{1}\left[{\bm{r}}_{1}\cdot(\bm{q}_{2}\times{\bm{r}}_{2})\right]-{\bm{r}}_{1}\left[\dot{\bm{r}}_{1}\cdot(\bm{q}_{2}\times{\bm{r}}_{2})\right]-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})(\bm{q}_{2}\times{\bm{r}}_{2}),$
$\displaystyle\frac{\partial\Phi_{4}}{\partial{\bm{r}}_{2}}=\dot{\bm{r}}_{2}\left[\dot{\bm{r}}_{2}\cdot(\bm{q}_{2}\times{\bm{r}}_{2})\right]+\left[\left(|\dot{\bm{r}}_{1}|^{2}-\frac{\mu}{|{\bm{r}}_{1}|}\right){\bm{r}}_{1}-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\dot{\bm{r}}_{1}+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2}\right]\times\bm{q}_{2},$
$\displaystyle\frac{\partial\Phi_{4}}{\partial\dot{\bm{r}}_{2}}={\bm{r}}_{2}\left[\dot{\bm{r}}_{2}\cdot(\bm{q}_{2}\times{\bm{r}}_{2})\right]+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})(\bm{q}_{2}\times{\bm{r}}_{2}),$
$\displaystyle\frac{\partial\Phi_{5}}{\partial{\bm{r}}_{1}}=-\dot{\bm{r}}_{1}\left[\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]-\frac{\mu}{|{\bm{r}}_{1}|^{3}}{\bm{r}}_{1}\left[{\bm{r}}_{2}\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]$
$\displaystyle\quad\qquad+({\bm{r}}_{2}-\bm{q}_{2})\times\left[-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\dot{\bm{r}}_{1}-\left(\frac{|\dot{\bm{r}}_{2}|^{2}}{2}-\frac{|\dot{\bm{r}}_{1}|^{2}}{2}+\frac{\mu}{|{\bm{r}}_{1}|}\right){\bm{r}}_{2}+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2}\right],$
$\displaystyle\frac{\partial\Phi_{5}}{\partial\dot{\bm{r}}_{1}}=-{\bm{r}}_{1}\left[\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\left[{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]+\dot{\bm{r}}_{1}\left[{\bm{r}}_{2}\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right],$
$\displaystyle\frac{\partial\Phi_{5}}{\partial{\bm{r}}_{2}}=-\left(\frac{|\dot{\bm{r}}_{2}|^{2}}{2}+\frac{|\dot{\bm{r}}_{1}|^{2}}{2}-\frac{\mu}{|{\bm{r}}_{1}|}\right)\left[{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]+\dot{\bm{r}}_{2}\left[\dot{\bm{r}}_{2}\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]$
$\displaystyle\quad\qquad+\left[-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\dot{\bm{r}}_{1}-\left(\frac{|\dot{\bm{r}}_{2}|^{2}}{2}+\frac{|\dot{\bm{r}}_{1}|^{2}}{2}-\frac{\mu}{|{\bm{r}}_{1}|}\right){\bm{r}}_{2}+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\dot{\bm{r}}_{2}\right]\times{\bm{r}}_{1},$
$\displaystyle\frac{\partial\Phi_{5}}{\partial\dot{\bm{r}}_{2}}=-\dot{\bm{r}}_{2}\left[{\bm{r}}_{2}\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]+{\bm{r}}_{2}\left[\dot{\bm{r}}_{2}\cdot{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right]+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\left[{\bm{r}}_{1}\times({\bm{r}}_{2}-\bm{q}_{2})\right].$

Equations (7) do not set any constraints on the mean anomalies $\ell_{1}$,
$\ell_{2}$ of the pairs of Keplerian orbits that we compute. To properly
select a solution we integrate back each orbit $\mathbf{E}_{\text{car}}^{(2)}$
computed at epoch $\tilde{t}_{2}$, to the epoch $\tilde{t}_{1}$ of the orbit
$\mathbf{E}_{\text{car}}^{(1)}$, together with its covariance matrix
$\Gamma^{(2)}$. Then, we compute a predicted position vector
$\mathcal{P}_{1,p}$ and its $3\times 3$ covariance matrix
$\Gamma_{\mathcal{P}_{1,p}}$. To define a metric to choose/discard solutions
we use a three-dimensional version of the attribution algorithm (Milani and
Gronchi, 2010), and introduce the identification penalty

$\chi_{3}^{2}=(\mathcal{P}_{1}-\mathcal{P}_{1,p})\cdot\left[C_{\mathcal{P}_{1,p}}-C_{\mathcal{P}_{1,p}}\Gamma_{0}C_{\mathcal{P}_{1,p}}\right](\mathcal{P}_{1}-\mathcal{P}_{1,p}),$

where

$C_{\mathcal{P}_{1,p}}=\Gamma_{\mathcal{P}_{1,p}}^{-1},\qquad\Gamma_{0}=C_{0}^{-1},$

with

$C_{0}=C_{\mathcal{P}_{1,p}}+C_{\mathcal{P}_{1}},\qquad
C_{\mathcal{P}_{1}}=\Gamma_{\mathcal{P}_{1}}^{-1}.$

We keep only the solutions with the lowest values of $\chi_{3}$ or, if we wish
to admit multiple solutions, we keep the ones whose value of $\chi_{3}$ is
below a certain threshold.

##  7 Numerical tests

In this section we show an application of our algorithm to asteroids
undergoing a close approch with the Earth. We consider all the near-Earth
asteroids present in the NEODyS database333http://newton.spacedys.com/neodys
to the date of December 6, 2023, and select the $1305$ asteroids with Earth
MOID444Minimum Orbit Intersection Distance, see Bowell and Muinonen (1994),
Gronchi et al. (2023) $d_{\rm min}$ smaller than $10^{-3}$ au. For each
selected asteroid we compute the points $P_{A}$ and $P_{\oplus}$, lying on the
osculating orbits of the asteroid and the Earth at epoch $t_{0}$, where the
minimum distance $d_{\rm min}$ is attained. Then, we find the closest epoch
$t_{1}>t_{0}$ when the Earth arrives at $P_{\oplus}$ and change the asteroid
phase so that it would arrive at $P_{A}$ at the same epoch $t_{1}$ in the
framework of the Sun-asteroid two-body dynamics. In this way, we try to
enhance the effect of the close approach with the Earth at $t_{1}$. A full
$n$-body propagation, from $t_{0}$ to $t_{1}$, is used to compute the
topocentric position vector ${\cal P}_{1}$ from the point of view of the Pan-
STARRS 1 telescope (et al., 2019). Continuing the propagation from $t_{1}$ to
a successive epoch $t_{2}$, we derive an attributable ${\cal A}_{2}$ by
coordinate change. As a result, the components $\alpha_{2}$, $\delta_{2}$,
$\dot{\alpha}_{2}$, $\dot{\delta}_{2}$ of ${\cal A}_{2}$ are not affected by
interpolation errors, which are usually introduced when an attributable is
obtained from the astrometric observations.

Assume we can model a close approach at epoch $t_{1}$ by an instantaneous
velocity change, like in Öpik’s theory (Öpik, 1976). Then, our algorithm can
be applied using ${\cal P}_{1},{\cal A}_{2}$ as input data, because the values
of the Keplerian integrals are conserved in $(t_{1},t_{2}]$. At $t_{1}$, when
the velocity instantaneously changes, the position vector ${\cal P}_{1}$
remains unchanged.

In practice, the close encounter is not instantaneous, but takes a finite
time: we select the epoch $t_{2}$ when the close approach is over, as
described below. In a typical close encounter, the value of Tisserand’s
parameter $T$ changes quickly, going back approximately to its pre-encounter
value when the encounter is over. In Cavallari et al. (2023, Sec. 3), for a
given value $C$ of the Jacobi constant, the authors provide a value $d_{C}$ of
the geocentric distance attaining $|\frac{dT}{dt}|\leq\varepsilon$, for a
fixed small parameter $\varepsilon$, so that we can think that the encounter
is over when the asteroid reaches that distance from the Earth. Following this
approach, we choose $t_{2}$ as the time when the value of the distance becomes
$d_{C}$.

Figure 1: Effects of the close approach on the semimajor axes (left) and on
the whole trajectories (right).

In Fig. 1 we show the effect of the close approach on these orbits. On the
left, we plot the distribution of the relative differences between the
semimajor axes $a_{1}$, $a_{2}$ obtained by the known orbits propagated at
epochs $t_{1}$ and $t_{2}$ using the software
OrbFit555http::adams.dm.unipi.it/orbfit. On the right, we plot the
distribution of the distance $\delta({\cal T}_{1},{\cal T}_{2})$ between the
propagated trajectories ${\cal T}_{1}$, ${\cal T}_{2}$, where

$\delta({\cal T}_{1},{\cal
T}_{2})=\sqrt{\frac{(a_{1}-a_{2})^{2}}{a_{1}^{2}}+(e_{1}-e_{2})^{2}+(i_{1}-i_{2})^{2}+(\Omega_{1}-\Omega_{2})^{2}+(\omega_{1}-\omega_{2})^{2}}.$
(19)

In (19) the subscripts $1$, $2$ of the Keplerian elements $a$, $e$, $i$,
$\Omega$, $\omega$ and the trajectories ${\cal T}$ refer to the epochs
$t_{1}$, $t_{2}$. In most cases the relative change in semimajor axis is
within 20%. However, there are a few cases where the change is larger, with
one extreme case (Fig 1, left) where it passes from $a_{1}\sim 2.84$ to
$a_{2}\sim 0.69$ au. The approximate values of the mean and the standard
deviation of $(a_{1}-a_{2})/a_{1}$ are about 0.012 and 0.0478, respectively.

Next, we show the performace of our algorithm using ${\cal P}_{1}$, ${\cal
A}_{2}$ as input data, which have been computed for each of these orbits by a
full $n$-body propagation with OrbFit. In case of multiple solutions, we
select the best one according to the procedure explained in Section 6.1.

In 8 cases, out of 1305, we could not obtain an orbital solution. However, by
increasing the time span $[t_{1},t_{2}]$ by $10\%$, while keeping the same
value of $t_{1}$, thus changing only ${\cal A}_{2}$, we recovered an orbital
solution in all these cases. In one case the computation of the second epoch
$t_{2}$ with Newton’s method failed. However, employing the starting guess
value for $t_{2}$, computed with a geocentric two-body dynamics, we obtained
an orbital solution also in this case.

In Fig. 2 we show the distribution of the differences between the Keplerian
elements computed with our algorithm (with subscript $c$) and the same
elements at $t_{2}$ obtained by propagation. In Fig. 2 (bottom right) we also
show the distribution of the distances $\delta({\cal T}_{2},{\cal T}_{c})$
between propagated (${\cal T}_{2}$) and computed (${\cal T}_{c}$)
trajectories, where the function $\delta$ is the same as in (19).

Figure 2: Distribution of the differences in orbital elements $a$ (top left),
$e$ (top right), $i$ (center left), $\Omega$ (center right), $\omega$ (bottom
left) between the propagated trajectory and the one computed with our
algorithm at epoch $t_{2}$. Distribution of the distances $\delta({\cal
T}_{2},{\cal T}_{c})$ between propagated and computed trajectories (bottom
right).

The results of this preliminary statistical test are satisfactory: the mean
and the standard deviation of the distributions displayed in Fig. 2 are shown
in Table 1.

| mean | std
---|---|---
$(a_{c}-a_{2})/a_{2}$ |  $-$5.0803$\times$10-4 | 0.0778
$e_{c}-e_{2}$ |  $-$0.003 | 0.0341
$i_{c}-i_{2}$ |  $-$3.8397$\times$10-4 | 0.0045
$\Omega_{c}-\Omega_{2}$ | 0.0019 | 0.0232
$\omega_{c}-\omega_{2}$ |  $-$0.0021 | 0.0244
$\delta({\cal T}_{c},{\cal T}_{2})$ | 0.0264 | 0.0877
Table 1: Mean and standard deviation of the distributions represented in Fig.
2.

Figure 3: Values of $q$ and $\upsilon$ in a log-log plot. The colors represent
the values of the distance $\delta({\cal T}_{2},{\cal T}_{c})$ between known
and computed trajectories.

Finally, we show the dependence of the results of our algorithm on the type of
close encounter. Following Cavallari et al. (2023, Sect. 3), close encounters
can be classified as deep/shallow and fast/slow, according to the values of
the minimum distance $q$ from the Earth center and the corresponding
geocentric velocity $\upsilon$. Here, we compute an approximation of $q$ and
$\upsilon$, still denoted by these symbols, using the following procedure: $q$
is the pericenter distance of the geocentric two-body orbit computed from the
Keplerian elements at $t_{1}$, and $\upsilon$ is the corresponding two-body
velocity. In Fig. 3 we plot the values of $q$ and $\upsilon$ using a color
scale representing the values of the trajectory distance $\delta({\cal
T}_{c},{\cal T}_{2})$.

For values of $q$ greater than $0.001$ au, we obtain better results (e.g.
lower values of the trajectory distance) for higher values of the velocity
$\upsilon$. For $q$ smaller than $0.001$ au, the results appear worse, and
there is no apparent correlation with velocity in the figure.

##  8 Conclusions

We have introduced a method to compute preliminary orbits with one topocentric
position vector ${\cal P}_{1}$ and a very short arc of optical observations,
from which we can derive an attributable ${\cal A}_{2}$. This method is based
on polynomial equations coming from the first integrals of Kepler’s problem,
i.e. angular momentum, energy, and Laplace-Lenz vector. Using the conservation
laws of these integrals, after introducing the auxiliary variable $z_{2}$, we
obtain a polynomial system that always has solutions, at least in the complex
field, even if ${\cal P}_{1}$, ${\cal A}_{2}$ do not correspond to the same
celestial object. There are some checks that can be performed to accept or
reject solutions.

We applied this algorithm to the computation of $1305$ NEA orbits, whose phase
was changed in order to enhance the close encounter effect with the Earth:
these preliminary results are satisfactory, see Section 7. The ideal situation
for the application of the proposed algorithm in this context would be given
by an instantaneous effect of the close encounter, which is not the case. The
selected epoch $t_{1}$ does not exactly correspond to the time of passage at
the MOID because we use a full $n$-body propagation starting from $t_{0}$, so
that at $t_{1}$ both the osculating trajectory and the time law along it may
have changed. Moreover, from this preliminary test, we saw that also varying
$t_{2}$ may affect the results. Therefore, a more detailed study is necessary
to better understand the applicability in case of close encounters and the
reliability of the computed orbits. On top of that, the sensitivity of the
algorithm to astrometric errors has still to be investigated.

Applications of the method introduced in this work to the orbit computation of
Earth satellites undergoing a maneuvre are also possible, and are worth to be
investigated in a future work.

#### Acknowledgments

The authors acknowledge the project MIUR-PRIN 20178CJA2AB “New Frontiers of
Celestial Mechanics: theory and applications”. G.F. Gronchi and G. Baù
acknowledge the project MSCA-ITN Stardust-R, Grant Agreement n. 813644 under
the H2020 research and innovation program. E. Scantamburlo acknowledges the
project “Advanced Space System Engineering to Address Broad Societal Benefits
– Starting Grant” funded by a contract between Politecnico di Torino and
Compagnia di San Paolo (CSP) 2019/2021 within the call “Attrazione e retention
di docenti di qualità”.

##  Appendix A Appendix

######  Lemma 2.

For ${\bm{c}}_{1}={\bm{c}}_{2}$, the generators $\mathfrak{q}_{4}$,
$\mathfrak{q}_{5}$, $\mathfrak{q}_{6}$ defined in (8) assume the following
form:

$\begin{split}\tilde{\mathfrak{q}}_{4}&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})\left[{\bm{r}}_{1}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]+({\bm{r}}_{2}\cdot{\bm{D}}_{1})z_{2},\cr\tilde{\mathfrak{q}}_{5}&=({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})\left[{\bm{r}}_{2}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]-({\bm{r}}_{1}\cdot{\bm{D}}_{2})\frac{\mu}{|{\bm{r}}_{1}|},\cr\tilde{\mathfrak{q}}_{6}&=({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})\left[{\bm{r}}_{1}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]-({\bm{r}}_{1}\cdot{\bm{D}}_{2})z_{2}.\end{split}$

###### Proof.

The result is given by a direct computation:

$\begin{split}\tilde{\mathfrak{q}}_{4}&=-(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\left[\dot{\bm{r}}_{1}\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})\right]-\mu\tilde{{\bm{L}}}_{2}\cdot{\bm{D}}_{1}\cr&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})-\mu\tilde{{\bm{L}}}_{2}\cdot{\bm{D}}_{1}\cr&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})-(\|{\bm{r}}_{2}\|^{2}-z_{2})\left[{\bm{r}}_{2}\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})\right]-(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\left[\dot{\bm{r}}_{2}\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})\right]\cr&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})-(\dot{\bm{r}}_{2}\times{\bm{c}}_{2})\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})+({\bm{r}}_{2}\cdot{\bm{D}}_{1})z_{2}\cr&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})-\left[\dot{\bm{r}}_{2}\times({\bm{r}}_{1}\times\dot{\bm{r}}_{1})\right]\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})+({\bm{r}}_{2}\cdot{\bm{D}}_{1})z_{2}\cr&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})-\left[{\bm{r}}_{1}(\dot{\bm{r}}_{1}\cdot\dot{\bm{r}}_{2})-\dot{\bm{r}}_{1}(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{1})\right]\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})+({\bm{r}}_{2}\cdot{\bm{D}}_{1})z_{2}\cr&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})+(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{1})\left[\dot{\bm{r}}_{1}\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{1})\right]+({\bm{r}}_{2}\cdot{\bm{D}}_{1})z_{2}\cr&=({\bm{c}}_{1}\cdot\hat{\bm{e}}^{\rho}_{1})\left[{\bm{r}}_{1}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]+({\bm{r}}_{2}\cdot{\bm{D}}_{1})z_{2},\end{split}$
$\begin{split}\tilde{\mathfrak{q}}_{5}&=\left(\dot{\bm{r}}_{1}\times{\bm{c}}_{1}-\frac{\mu}{|{\bm{r}}_{1}|}\right)\cdot({\bm{r}}_{2}\times\hat{\bm{e}}^{\rho}_{2})-({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\cr&=\left(\dot{\bm{r}}_{1}\times{\bm{c}}_{2}-\frac{\mu}{|{\bm{r}}_{1}|}\right)\cdot({\bm{r}}_{2}\times\hat{\bm{e}}^{\rho}_{2})-({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\cr&=\left((\dot{\bm{r}}_{1}\cdot\dot{\bm{r}}_{2}){\bm{r}}_{2}-({\bm{r}}_{2}\cdot\dot{\bm{r}}_{1})\dot{\bm{r}}_{2}-\frac{\mu}{|{\bm{r}}_{1}|}\right)\cdot({\bm{r}}_{2}\times\hat{\bm{e}}^{\rho}_{2})-({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\cr&=-({\bm{r}}_{2}\cdot\dot{\bm{r}}_{1})\left[\dot{\bm{r}}_{2}\cdot({\bm{r}}_{2}\times\hat{\bm{e}}^{\rho}_{2})\right]-\frac{\mu}{|{\bm{r}}_{1}|}{\bm{r}}_{1}\cdot({\bm{r}}_{2}\times\hat{\bm{e}}^{\rho}_{2})-({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{2})\cr&=({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})\left[{\bm{r}}_{2}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]-({\bm{r}}_{1}\cdot{\bm{D}}_{2})\frac{\mu}{|{\bm{r}}_{1}|},\end{split}$
$\begin{split}\tilde{\mathfrak{q}}_{6}&=-\left[(\dot{\bm{r}}_{2}\times{\bm{c}}_{2})-z_{2}{\bm{r}}_{2}\right]\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{2})+({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\cr&=-\left[(\dot{\bm{r}}_{2}\times{\bm{c}}_{1})-z_{2}{\bm{r}}_{2}\right]\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{2})+({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\cr&=-\left[(\dot{\bm{r}}_{1}\cdot\dot{\bm{r}}_{2}){\bm{r}}_{1}-(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{1})\dot{\bm{r}}_{1}-z_{2}{\bm{r}}_{2}\right]\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{2})+({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\cr&=(\dot{\bm{r}}_{2}\cdot{\bm{r}}_{1})\left[\dot{\bm{r}}_{1}\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{2})\right]+z_{2}\left[{\bm{r}}_{2}\cdot({\bm{r}}_{1}\times\hat{\bm{e}}^{\rho}_{2})\right]+({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})(\dot{\bm{r}}_{1}\cdot{\bm{r}}_{1})\cr&=({\bm{c}}_{2}\cdot\hat{\bm{e}}^{\rho}_{2})\left[{\bm{r}}_{1}\cdot(\dot{\bm{r}}_{1}-\dot{\bm{r}}_{2})\right]-({\bm{r}}_{1}\cdot{\bm{D}}_{2})z_{2}.\end{split}$

∎

## References

  * Bowell and Muinonen (1994) E. Bowell and K. Muinonen.  _Hazards due to comets and asteroids_.  Tucson: The University of Arizona Press, 1994.
  * Cavallari et al. (2023) I. Cavallari, C. Grassi, G. F. Gronchi, G. Baù, and G. B. Valsecchi.  A dynamical definition of the sphere of influence of the Earth.  _Communications in Nonlinear Science and Numerical Simulation_ , 119:107091, 2023.
  * et al. (2019) K. C. Chambers et al.  The Pan-STARRS1 surveys, 2019.
  * Gauss (1809) C. F. Gauss.  _Theoria motus corporum in sectionibus conicis solem ambientium_.  reprinted by Dover publications in 1963, 1809.
  * Gronchi et al. (2010) G. F. Gronchi, L. Dimare, and A. Milani.  Orbit determination with the two-body integrals.  _Celestial Mechanics and Dynamical Astronomy_ , 107:299–318, 2010.
  * Gronchi et al. (2011) G. F. Gronchi, D. Farnocchia, and L. Dimare.  Orbit determination with the two-body integrals. II.  _Celestial Mechanics and Dynamical Astronomy_ , 110:257–270, 2011.
  * Gronchi et al. (2015) G. F. Gronchi, G. Baù, and S. Marò.  Orbit determination with the two-body integrals. III.  _Celestial Mechanics and Dynamical Astronomy_ , 123:105–122, 2015.
  * Gronchi et al. (2017) G. F. Gronchi, G. Baù, and A. Milani.  Keplerian integrals, elimination theory and identification of very short arcs in a large database of optical observations.  _Celestial Mechanics and Dynamical Astronomy_ , 127:211–232, 2017.
  * Gronchi et al. (2023) G. F. Gronchi, G. Baù, and C. Grassi.  Revisiting the computation of the critical points of the Keplerian distance.  _Celestial Mechanics and Dynamical Astronomy_ , 135, 2023.
  * Lagrange (1783) J. L. Lagrange.  Sur le problème de la détermination des orbites des cométes d’aprés trois observations. Troisiè memémoire.  _Nouveaux mémoires de l’Académie royale des sciences et belles-lettres de Berlin_ , 1783.  Reprinted in Œuvres de Lagrange, Gauthier-Villars et fils, Paris (1869), volume 4, pp. 496–532.
  * Laplace (1780) P. S. Laplace.  Mémoire sur la détermination des orbites des comètes.  _Mémoires de l’Académie royale des sciences de Paris_ , 1780.  Reprinted in Œuvres complètes de Laplace, Gauthier-Villars et fils, Paris (1894), volume 10, pp. 93–146.
  * Milani and Gronchi (2010) A. Milani and G. F. Gronchi.  _Theory of Orbit Determination_.  Cambridge Univ. Press, 2010.
  * Milani et al. (2005) A. Milani, G. F. Gronchi, Z. Knežević, M. E. Sansaturio, and O. Arratia.  Orbit Determination with Very Short Arcs. II Identifications.  _Icarus_ , 179:350–364, 2005.
  * Öpik (1976) E. J. Öpik.  _Interplanetary encounters: close-range gravitational interactions_.  Elsevier, Amsterdam, 1976.
  * Taff and Hall (1977) L. G. Taff and D. L. Hall.  The use of angles and angular rates. I - Initial Orbit Determination.  _Celestial Mechanics_ , 16:481–488, 1977.