# Generalization of a method by Mossotti for initial orbit determination

Giovanni F. Gronchi  Dipartimento di Matematica, Università di Pisa, Italy
Giulio Baù  Dipartimento di Matematica, Università di Pisa, Italy  Óscar
Rodríguez  Dipartimento di Matematica, Università di Pisa, Italy  Robert
Jedicke  Institute for Astronomy, University of Hawai‘i, USA  Joachim Moeyens
Vera C. Rubin Observatory, University of Washington, USA

###### Abstract

Here we revisit an initial orbit determination method introduced by O. F.
Mossotti employing four geocentric sky-plane observations and a linear
equation to compute the angular momentum of the observed body. We then extend
the method to topocentric observations, yielding a quadratic equation for the
angular momentum. The performance of the two versions are compared through
numerical tests with synthetic asteroid data using different time intervals
between consecutive observations and different astrometric errors. We also
show a comparison test with Gauss’s method using simulated observations with
the expected cadence of the VRO-LSST telescope.

##  1 Introduction

In 1816 Ottaviano F. Mossotti introduced a method for initial orbit
determination of a solar system body employing four optical observations, e.g.
the values of right ascension and declination. Assuming geocentric
observations, Mossotti’s method allows to write linear equations for the
computation of the orbital angular momentum [17]. Then the orbit can be
reconstructed, e.g. by Gibbs’ method [10]. This procedure has the advantage to
avoid the computation of the roots of the eight degree polynomial appearing in
the classical methods by Laplace [13], Lagrange [12], and Gauss [5], which
need only three observations but can give rise to multiple solutions. A review
of the methods by Laplace and Gauss together with a geometric interpretation
of the occurrence of multiple solutions can be found in [7], [14].

Mossotti’s work was appreciated by Gauss himself, see [6]. This method has
been reviewed in [1], where the authors state that the computation of the
solution can be seriously affected by the observational errors due to the
terms that are neglected in the employed approximation.

In this work we recall Mossotti’s original method and show that it is possible
to define a topocentric version which leads to a quadratic equation for the
angular momentum. This generalization of the method turns out to be suitable
for orbit determination of Earth satellites too. We investigate the
performance of the methods and their sensitivity to observational errors by
some numerical tests with simulated data: we compare the original geocentric
method with this topocentric version using different time intervals between
the observations and different astrometric errors. We also show a comparison
test with Gauss’s method using simulated observations with the expected
cadence of the VRO-LSST telescope.

##  2 The original method

Mossotti’s method [17] leads to a set of linear equations for the components
of ${\bm{c}}_{\oplus}-{\bm{c}}$, where ${\bm{c}}_{\oplus}$ and ${\bm{c}}$ are
the angular momenta of the Earth and a solar system body, respectively. The
observations are supposed to be made from the center of the Earth. We assume
that the observed body is an asteroid, moving along an elliptic Keplerian
trajectory with the Sun as the center of force. We also assume that the total
observational arc is covered in a much shorter time than the orbital period.
This method is also suitable to be used with hyperbolic or parabolic orbits.

Here we illustrate all the formulae which are necessary for a numerical
implementation following Mossotti’s paper steps [17]. However, in the early
XIXth century Linear Algebra had not been developed yet, and several formulae
in [17] can be written and derived in a shorter way. The original formulae can
be recovered using Table 3 in Appendix A.

###  2.1 Units and preliminary definitions

In order to simplify the notation we use the rescaled time $\theta$, defined
by

$\theta=t\sqrt{g\left(1+\frac{m_{\oplus}}{m_{\odot}}\right)},$

where $m_{\oplus}$ and $m_{\odot}$ are the masses of the Earth and the Sun,
and $g=Gm_{\odot}$, where $G$ is Newton’s gravitational constant. We also use
$\kappa$ to denote Gauss’s constant $\sqrt{g}$. In the following we assume
that

$\theta\approx\kappa t,$

neglecting the constant ${m_{\oplus}}/{m_{\odot}}$.

Take three of the four observations, at epochs $t_{1}<t_{2}<t_{3}$, and set

$\theta_{ij}=\kappa(t_{j}-t_{i}),$
$\bm{\theta}=(\theta_{23},\,\theta_{31},\,\theta_{12})^{t},$

where the superscript $t$ stands for transposition. We write ${\bm{c}}$ and
${\bm{c}}_{\oplus}$ for the orbital angular momenta of the asteroid and the
Earth, respectively, and introduce their unit vectors

$\hat{{\bm{c}}}=\frac{{\bm{c}}}{c},\qquad\hat{{\bm{c}}}_{\oplus}=\frac{{\bm{c}}_{\oplus}}{c_{\oplus}},$

where $c=|{\bm{c}}|$, $c_{\oplus}=|{\bm{c}}_{\oplus}|$, being $|\bm{x}|$ the
Euclidean norm of a vector $\bm{x}$.

We also write $\bm{r}_{i}$ and $\bm{q}_{i}$ for the heliocentric positions of
the asteroid and the Earth at the three epochs $t_{i}$ ($i=1,2,3$), use
${\bm{\rho}}_{i}=\bm{r}_{i}-\bm{q}_{i}$ for the geocentric position of the
asteroid, and introduce the unit vectors

$\hat{\bm{q}}_{i}=\frac{\bm{q}_{i}}{q_{i}},\qquad{\hat{\bm{\rho}}}_{i}=\frac{{\bm{\rho}}_{i}}{\rho_{i}},\qquad
i=1,2,3,$

where $q_{i}=|\bm{q}_{i}|$, $\rho_{i}=|{\bm{\rho}}_{i}|$. Finally, we denote
the parameters111in [17] these are called semiparametri. of the orbits of the
asteroid and the Earth by $p$, $p_{\oplus}$. They are defined by

$p\kappa^{2}=c^{2},\qquad p_{\oplus}\kappa^{2}=c_{\oplus}^{2}.$

###  2.2 Geometric relations

With the purpose of writing Mossotti’s equations in a compact form, let us
introduce the matrices

$P=(\bm{\rho}_{1}\,|\,\bm{\rho}_{2}\,|\,\bm{\rho}_{3}),\qquad
Q=(\bm{q}_{1}\,|\,\bm{q}_{2}\,|\,\bm{q}_{3}),\qquad
R=P+Q=(\bm{r}_{1}\,|\,\bm{r}_{2}\,|\,\bm{r}_{3}),$ (1)

where $(\bm{x}_{1}\,|\,\bm{x}_{2}\,|\,\bm{x}_{3})$ is the matrix whose columns
are the vectors $\bm{x}_{j}$. The corresponding adjugate matrices are

$\begin{split}\text{adj}(P)&=\bigl{(}\bm{\rho}_{2}\times\bm{\rho}_{3}\,|\,\bm{\rho}_{3}\times\bm{\rho}_{1}\,|\,\bm{\rho}_{1}\times\bm{\rho}_{2}\bigr{)}^{t},\cr\text{adj}(Q)&=\bigl{(}\bm{q}_{2}\times\bm{q}_{3}\,|\,\bm{q}_{3}\times\bm{q}_{1}\,|\,\bm{q}_{1}\times\bm{q}_{2}\bigr{)}^{t},\cr\text{adj}(R)&=\bigl{(}\bm{r}_{2}\times\bm{r}_{3}\,|\,\bm{r}_{3}\times\bm{r}_{1}\,|\,\bm{r}_{1}\times\bm{r}_{2}\bigr{)}^{t}.\cr\end{split}$

We recall the following property, which holds for any square matrix $M$:

$M\,\text{adj}(M)=\text{adj}(M)\,M=\det(M)\,I,$ (2)

where $I$ is the identity matrix.

Figure 1: Geometry of the three observations.

The rank of the matrices $Q$ and $R$ is 2, since each triplet
$\\{\bm{q}_{1},\bm{q}_{2},\bm{q}_{3}\\}$ and
$\\{\bm{r}_{1},\bm{r}_{2},\bm{r}_{3}\\}$ is made by coplanar vectors and we
assume that

$\bm{q}_{i}\times\bm{q}_{j}\neq\bm{0},\qquad\bm{r}_{i}\times\bm{r}_{j}\neq\bm{0},\qquad
1\leq i<j\leq 3,$

see Figure 1. This implies

$Q\,\text{adj}(Q)=0,\qquad R\,\text{adj}(R)=0.$ (3)

Since the angular momenta ${\bm{c}}$, ${\bm{c}}_{\oplus}$ are respectively
orthogonal to the orbital planes of the asteroid and the Earth, we have

${\bm{c}}\cdot\bm{r}_{i}={\bm{c}}\cdot(\bm{q}_{i}+{\bm{\rho}}_{i})=0,\quad{\bm{c}}_{\oplus}\cdot\bm{q}_{i}=0,$

which lead us to

${\bm{c}}\cdot{\bm{\rho}}_{i}=({\bm{c}}_{\oplus}-{\bm{c}})\cdot\bm{q}_{i},$
(4)

for $i=1,2,3$. Then, we introduce the vectors

$\bm{\tau}=(\tau_{1},\,\tau_{2},\,\tau_{3})^{t},\qquad\bm{T}=(T_{1},\,T_{2},\,T_{3})^{t},$

with

$\tau_{1}\sqrt{p}={\bm{r}}_{2}\times{\bm{r}}_{3}\cdot\hat{{\bm{c}}},\hskip
28.45274ptT_{1}\sqrt{p_{\oplus}}=\bm{q}_{2}\times\bm{q}_{3}\cdot\hat{{\bm{c}}}_{\oplus},$
(5a)
$\tau_{2}\sqrt{p}={\bm{r}}_{3}\times{\bm{r}}_{1}\cdot\hat{{\bm{c}}},\hskip
28.45274ptT_{2}\sqrt{p_{\oplus}}=\bm{q}_{3}\times\bm{q}_{1}\cdot\hat{{\bm{c}}}_{\oplus},$
(5b)
$\tau_{3}\sqrt{p}={\bm{r}}_{1}\times{\bm{r}}_{2}\cdot\hat{{\bm{c}}},\hskip
28.45274ptT_{3}\sqrt{p_{\oplus}}=\bm{q}_{1}\times\bm{q}_{2}\cdot\hat{{\bm{c}}}_{\oplus}.$
(5c)

Noting that

$\bm{\tau}=\frac{1}{\sqrt{p}}\text{adj}(R)\hat{\bm{c}},\qquad\bm{T}=\frac{1}{\sqrt{p_{\oplus}}}\text{adj}(Q)\hat{\bm{c}}_{\oplus},$
(6)

by (3) we have

$R\,\bm{\tau}=\bm{0},\qquad Q\,\bm{T}=\bm{0}.$

Moreover, recalling that $R=Q+P$, we get

$P\bm{\tau}=Q(\bm{T}-\bm{\tau}).$ (7)

Since
$(\bm{r}_{i}\times\bm{r}_{j})\times{\bm{c}}=(\bm{q}_{i}\times\bm{q}_{j})\times{\bm{c}}_{\oplus}={\bm{0}}$,
we have

$\displaystyle\left(\bm{r}_{i}\times\bm{r}_{j}\cdot\hat{\bm{c}}\right)\left(\hat{\bm{c}}\cdot\bm{\rho}_{i}\right)$
$\displaystyle=\bm{r}_{i}\times\bm{r}_{j}\cdot\bm{\rho}_{i}=\bm{q}_{i}\times\bm{q}_{j}\cdot\bm{\rho}_{i}-\bm{\rho}_{i}\times\bm{\rho}_{j}\cdot\bm{q}_{i},$
(8)
$\displaystyle(\bm{q}_{i}\times\bm{q}_{j}\cdot\hat{{\bm{c}}}_{\oplus})(\hat{{\bm{c}}}_{\oplus}\cdot\bm{\rho}_{i})$
$\displaystyle=\bm{q}_{i}\times\bm{q}_{j}\cdot\bm{\rho}_{i}.$ (9)

######  Remark 1.

If the orbits of the Earth and the asteroid are almost coplanar, then the
equations (8), (9) are close to be degenerate. This singularity is a common
feature of orbit determination methods, where the value of the geodesic
curvature of the observed arc plays an important role, see [14, Chap. 9].

The last equality in (8) is obtained by observing that

$\bm{r}_{i}\times\bm{r}_{j}\cdot\bm{\rho}_{i}=\bm{q}_{i}\times\bm{r}_{j}\cdot\bm{\rho}_{i}=\bm{q}_{i}\times\bm{q}_{j}\cdot\bm{\rho}_{i}+\bm{q}_{i}\times\bm{\rho}_{j}\cdot\bm{\rho}_{i}.$

By (8), (9) and the definitions of the $\tau_{k}$, $T_{k}$ given in (5c) we
obtain

$\displaystyle\varepsilon_{ijk}\frac{\tau_{k}}{\kappa}\bm{c}\cdot\bm{\rho}_{i}$
$\displaystyle=\bm{q}_{i}\times\bm{q}_{j}\cdot\bm{\rho}_{i}-\bm{\rho}_{i}\times\bm{\rho}_{j}\cdot\bm{q}_{i},$
(10)
$\displaystyle\varepsilon_{ijk}\frac{T_{k}}{\kappa}\bm{c}_{\oplus}\cdot\bm{\rho}_{i}$
$\displaystyle=\bm{q}_{i}\times\bm{q}_{j}\cdot\bm{\rho}_{i},$ (11)

where $\varepsilon_{ijk}$ denotes the Levi-Civita symbol, and the indexes
$i,j,k$ vary so that all the 6 permutations of the set $\\{1,2,3\\}$ can be
considered.

Subtracting (10) from (11) we obtain

$\varepsilon_{ijk}\frac{1}{\kappa}(T_{k}{\bm{c}}_{\oplus}-\tau_{k}{\bm{c}})\cdot\bm{\rho}_{i}=\bm{\rho}_{i}\times\bm{\rho}_{j}\cdot\bm{q}_{i}.$
(12)

Writing

$T_{k}\bm{c}_{\oplus}-\tau_{k}\bm{c}=T_{k}\left(\bm{c}_{\oplus}-\bm{c}\right)+\left(T_{k}-\tau_{k}\right)\bm{c},$

and using (4), we have

$\displaystyle\left[T_{k}\left(\bm{c}_{\oplus}-\bm{c}\right)+\left(T_{k}-\tau_{k}\right)\bm{c}\right]\cdot\bm{\rho}_{i}$
$\displaystyle=T_{k}\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\bm{\rho}_{i}+\left(T_{k}-\tau_{k}\right)\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\bm{q}_{i}$
$\displaystyle=\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\left[\hat{\bm{\rho}}_{i}+\left(1-\frac{\tau_{k}}{T_{k}}\right)\frac{q_{i}}{\rho_{i}}\hat{\bm{q}}_{i}\right]T_{k}\rho_{i}.$

Therefore, introducing the coefficients

$\alpha_{ik}=\left(1-\frac{\tau_{k}}{T_{k}}\right)\frac{q_{i}}{\rho_{i}},$
(13)

and simplifying $\rho_{i}$ we can write relations (12) as

$\varepsilon_{ijk}\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\left(\hat{\bm{\rho}}_{i}+\alpha_{ik}\hat{\bm{q}}_{i}\right)T_{k}=\kappa\rho_{j}\hat{\bm{\rho}}_{i}\times\hat{\bm{\rho}}_{j}\cdot\bm{q}_{i}.$
(14)

###  2.3 Combining geometry of observations with two-body dynamics

Since the time intervals $\theta_{12}$, $\theta_{23}$ are small compared to
the orbital period,222 the observations of asteroids at Mossotti’s epoch were
no more than one per night, and the time interval between two of them covered
a few days. we can consider Taylor’s expansions of the position of the
asteroid and the center of the Earth in their orbital planes as power series
of $\theta_{23}$, $\theta_{31}$, $\theta_{12}$ centering the expansion at the
intermediate epoch. Then, we neglect the terms depending on the powers of
$\theta_{ij}$ greater than 3. In this way we obtain

$\bm{\tau}\simeq\bm{\theta}-\frac{1}{6r_{2}^{3}}\bm{\theta}^{3},\quad\quad\bm{T}\simeq\bm{\theta}-\frac{1}{6q_{2}^{3}}\bm{\theta}^{3},$
(15)

where $r_{2}=|\bm{r}_{2}|$ and

$\bm{\theta}^{3}=\bm{\theta}\odot\bm{\theta}\odot\bm{\theta}=\left(\theta_{23}^{3},\,\theta_{31}^{3},\,\theta_{12}^{3}\right)^{t},$

with $\odot$ denoting the Hadamard product.333 if
$\bm{a}=(a_{1},a_{2},a_{3})^{t}$ and $\bm{b}=(b_{1},b_{2},b_{3})^{t}$, then
$\bm{a}\odot\bm{b}=(a_{1}b_{1},a_{2}b_{2},a_{3}b_{3})^{t}$.

Let us define

$\bm{u}=(u_{1},\,u_{2},\,u_{3})^{t}=\text{adj}(\hat{P})Q{\bm{\theta}^{3}},$

with
$\hat{P}=(\hat{\bm{\rho}}_{1}\,|\,\hat{\bm{\rho}}_{2}\,|\,\hat{\bm{\rho}}_{3})$.
Multiplying relation (7) on the left by $\text{adj}(\hat{P})$, and taking into
account the approximations (15), we get

$\displaystyle\text{adj}(\hat{P})P\bm{\theta}$
$\displaystyle\simeq\text{adj}(\hat{P})P\bm{\tau}=\text{adj}(\hat{P})Q(\bm{T}-\bm{\tau})$
(16)
$\displaystyle\simeq\text{adj}(\hat{P})Q\bm{\theta}^{3}\left(\frac{1}{6r_{2}^{3}}-\frac{1}{6q_{2}^{3}}\right)=\frac{1}{6}\left(\frac{1}{r_{2}^{3}}-\frac{1}{q_{2}^{3}}\right)\bm{u}.$

Recalling (2), and noting that

$\text{adj}(\hat{P})P=\det(\hat{P})\textrm{diag}\\{\rho_{1},\rho_{2},\rho_{3}\\},$

relation (16) yields

$\det(\hat{P})\bm{\delta}\odot\bm{\theta}\simeq\frac{1}{6}\left(\frac{1}{r_{2}^{3}}-\frac{1}{q_{2}^{3}}\right)\bm{u},$
(17)

with

$\bm{\delta}=(\rho_{1},\rho_{2},\rho_{3})^{t}.$

###  2.4 A linear equation involving ${\bm{c}}_{\oplus}-{\bm{c}}$

Choosing $(i,j,k)=(1,2,3),(3,2,1)$ in (14) and eliminating $\kappa\rho_{2}$
from the two resulting equations, we obtain

$\begin{split}&\left(\hat{\bm{\rho}}_{1}\times\hat{\bm{\rho}}_{2}\cdot\bm{q}_{1}\right)\left[\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\left(\hat{\bm{\rho}}_{3}+\alpha_{31}\hat{\bm{q}}_{3}\right)\right]T_{1}=\\\
&\left(\hat{\bm{\rho}}_{2}\times\hat{\bm{\rho}}_{3}\cdot\bm{q}_{3}\right)\left[\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\left(\hat{\bm{\rho}}_{1}+\alpha_{13}\hat{\bm{q}}_{1}\right)\right]T_{3}.\end{split}$
(18)

In this equation the only unknowns different from
$\left(\bm{c}_{\oplus}-\bm{c}\right)$ are the coefficients $\alpha_{13}$ and
$\alpha_{31}$. Using the approximations of $\bm{\tau}$ and $\bm{T}$ given by
(15) in (13) we have

$\alpha_{13}=\frac{q_{1}}{\rho_{1}}\frac{T_{3}-\tau_{3}}{T_{3}}\simeq\frac{q_{1}}{\rho_{1}T_{3}}\frac{\theta_{12}^{3}}{6}\left(\frac{1}{r_{2}^{3}}-\frac{1}{q_{2}^{3}}\right)\simeq\frac{q_{1}}{\rho_{1}}\frac{\theta_{12}^{2}}{6}\left(\frac{1}{r_{2}^{3}}-\frac{1}{q_{2}^{3}}\right),$
(19)

where we used $\frac{1}{T_{3}}=\frac{1}{\theta_{12}}(1+O(\theta_{12}))$. In a
similar way we obtain

$\alpha_{31}\simeq\frac{q_{3}}{\rho_{3}}\frac{\theta_{23}^{2}}{6}\left(\frac{1}{r_{2}^{3}}-\frac{1}{q_{2}^{3}}\right).$
(20)

Inserting the approximations of $\rho_{1},\rho_{3}$ given by (17) into (19),
(20) we can express $\alpha_{13}$, $\alpha_{31}$ with known quantities:

$\alpha_{13}\simeq\frac{\det(\hat{P})q_{1}{\theta_{12}^{2}\theta_{23}}}{u_{1}},\qquad\alpha_{31}\simeq\frac{\det(\hat{P})q_{3}{\theta_{23}^{2}\theta_{12}}}{u_{3}}.$
(21)

Defining the coefficients

$a_{1}=\frac{q_{2}\left(\hat{\bm{\rho}}_{2}\times\hat{\bm{\rho}}_{3}\cdot\bm{q}_{3}\right)}{T_{1}\sqrt{p_{\oplus}}}=\frac{\hat{\bm{\rho}}_{2}\times\hat{\bm{\rho}}_{3}\cdot\hat{\bm{q}}_{3}}{\hat{\bm{q}}_{2}\times\hat{\bm{q}}_{3}\cdot\hat{\bm{c}}_{\oplus}},\qquad
a_{3}=\frac{q_{2}\left(\hat{\bm{\rho}}_{1}\times\hat{\bm{\rho}}_{2}\cdot\bm{q}_{1}\right)}{T_{3}\sqrt{p_{\oplus}}}=\frac{\hat{\bm{\rho}}_{1}\times\hat{\bm{\rho}}_{2}\cdot\hat{\bm{q}}_{1}}{\hat{\bm{q}}_{1}\times\hat{\bm{q}}_{2}\cdot\hat{\bm{c}}_{\oplus}},$

and the vectors

$\bm{\gamma}=a_{1}\left(\hat{\bm{\rho}}_{1}+\alpha_{13}\hat{\bm{q}}_{1}\right),\qquad\bm{\varphi}=a_{3}\left(\hat{\bm{\rho}}_{3}+\alpha_{31}\hat{\bm{q}}_{3}\right),$

equation (18) becomes

$(\bm{\gamma}-\bm{\varphi})\cdot\left(\bm{c}_{\oplus}-\bm{c}\right)=0.$ (22)

######  Remark 2.

Choosing the index pair $\\{(2,3,1),(1,3,2)\\}$ or $\\{(3,1,2),(2,1,3)\\}$ in
(14) we can obtain two additional equations analogous to (22). However, with
the employed approximation of $\bm{T}$ and $\bm{\tau}$, the three equations
are just the same (see [17, Sect. 28]).

###  2.5 Mossotti’s equations for ${\bm{c}}_{\oplus}-{\bm{c}}$

With the aim of writing two independent linear equations, all the four
observations are used. If we consider two different choices of the three
observations, out of the available four, we obtain the system

$\left\\{\begin{aligned}
&(\bm{\gamma}_{1}-\bm{\varphi}_{1})\cdot\left(\bm{c}_{\oplus}-\bm{c}\right)=0\\\
&(\bm{\gamma}_{2}-\bm{\varphi}_{2})\cdot\left(\bm{c}_{\oplus}-\bm{c}\right)=0,\end{aligned}\right.$
(23)

where the subscripts $1$, $2$ of $\bm{\gamma}$ and $\bm{\varphi}$ refer to the
two triplets of observations. Set

$\bm{w}=(\bm{\gamma}_{1}-\bm{\varphi}_{1})\times(\bm{\gamma}_{2}-\bm{\varphi}_{2})$

and assume $\bm{w}\neq\bm{0}$. Then the general solution of (23) has the form

$\bm{c}_{\oplus}-\bm{c}=\lambda\bm{w},$ (24)

with $\lambda\in\mathbb{R}$, giving the direction of $\bm{c}_{\oplus}-\bm{c}$.

######  Remark 3.

If $\bm{\gamma}_{1}-\bm{\varphi}_{1}$ and $\bm{\gamma}_{2}-\bm{\varphi}_{2}$
are almost parallel, then system (23) is almost degenerate: we can try to
avoid this singularity by choosing other triplets of observations.

In order to constrain the values of $\lambda$ we proceed as follow. Choosing
$(i,j,k)=(1,2,3)$ in equation (14), we have

$\frac{a_{3}c_{\oplus}\rho_{2}}{q_{2}}=(\bm{c}_{\oplus}-\bm{c})\cdot\left(\hat{\bm{\rho}}_{1}+\alpha_{13}\hat{\bm{q}}_{1}\right).$
(25)

Inserting the expression (24) of the general solution in (25) and (4) with
$i=2$ we obtain, respectively,

$\rho_{2}=\lambda\frac{q_{2}}{a_{1}a_{3}c_{\oplus}}\bm{w}\cdot\bm{\gamma}_{1},$
(26)

and

$\lambda\bm{w}\cdot\bm{q}_{2}=(\bm{c}_{\oplus}-\lambda\bm{w})\cdot\bm{\rho}_{2}=(\bm{c}_{\oplus}-\lambda\bm{w})\cdot\hat{\bm{\rho}}_{2}\rho_{2}.$
(27)

Substituting the expression (26) of $\rho_{2}$ in (27) yields a quadratic
equation in $\lambda$, which is here the only unknown:

$\lambda\bm{w}\cdot\bm{q}_{2}=\lambda\frac{q_{2}}{a_{1}a_{3}c_{\oplus}}\left(\bm{w}\cdot\bm{\gamma}_{1}\right)\big{[}(\bm{c}_{\oplus}-\lambda\bm{w})\cdot\hat{\bm{\rho}}_{2}\big{]}.$

This equation can be written as

$(\bm{w}\cdot\bm{\gamma}_{1})(\bm{w}\cdot\hat{\bm{\rho}}_{2})\lambda^{2}+\big{[}a_{1}a_{3}c_{\oplus}(\bm{w}\cdot\hat{\bm{q}}_{2})-(\bm{w}\cdot\bm{\gamma}_{1})(\bm{c}_{\oplus}\cdot\hat{\bm{\rho}}_{2})\big{]}\lambda=0,$
(28)

whose solutions are

$\lambda=0,$ (29a)
$\lambda=\frac{(\bm{w}\cdot\bm{\gamma}_{1})(\bm{c}_{\oplus}\cdot\hat{\bm{\rho}}_{2})-a_{1}a_{3}c_{\oplus}(\bm{w}\cdot\hat{\bm{q}}_{2})}{(\bm{w}\cdot\bm{\gamma}_{1})(\bm{w}\cdot\hat{\bm{\rho}}_{2})}.$
(29b)

Substituting these expressions in (24) gives two possible values of the
angular momentum ${\bm{c}}$. The solution (29a) yields
${\bm{c}}_{\oplus}={\bm{c}}$, and is usually discarded, so that (29b) is
regarded as the only solution. In this way the equations of Mossotti’s method
can be considered linear.

##  3 Topocentric method

We first introduce some notation. Let us define $\bm{q}_{\oplus}$,
$\bm{p}_{\rm obs}$ as the heliocentric position of the Earth center, and the
geocentric position of the observer, respectively. The heliocentric positions
of the observer and asteroid are

$\displaystyle\bm{q}=\bm{q}_{\oplus}+\bm{p}_{\rm obs},$
$\displaystyle{\bm{r}}=\bm{q}+{\bm{\rho}}=\bm{q}_{\oplus}+{\bm{\rho}}_{\rm
geo},$

where ${\bm{\rho}}$, ${\bm{\rho}}_{\rm geo}$ are the topocentric and
geocentric positions of the asteroid, respectively (see Figure 2).

SunEarthAsteroid$\bm{q}_{\oplus}$$\bm{q}$$\bm{p}_{\rm
obs}$$\bm{r}$${\bm{\rho}}$${\bm{\rho}}_{\rm geo}$ Figure 2: Geocentric and
topocentric point of view.

###  3.1 Geometric relations

As in the geocentric case, we select three observations of the asteroid out of
the available four, and introduce the matrices $P$, $Q$, $R$ as in (1), but
with a different interpretation for the vectors $\bm{\rho}$ and $\bm{q}$: here
$\bm{\rho}$ represents the topocentric position of the asteroid, and $\bm{q}$
gives the heliocentric position of the observer. Moreover, we introduce the
matrices

$Q_{\oplus}={\bigl{(}\bm{q}_{\oplus,1}\,|\,\bm{q}_{\oplus,2}\,|\,\bm{q}_{\oplus,3}\bigr{)}},\qquad
P_{\text{obs}}={\bigl{(}\bm{p}_{\text{obs},1}\,|\,\bm{p}_{\text{obs},2}\,|\,\bm{p}_{\text{obs},3}\bigr{)}}.$

We recall the geometrical relations

${\bm{c}}\cdot\bm{r}_{i}={\bm{c}}\cdot(\bm{q}_{i}+{\bm{\rho}}_{i})=0,\quad{\bm{c}}_{\oplus}\cdot\bm{q}_{\oplus,i}={\bm{c}}_{\oplus}\cdot(\bm{q}_{i}-\bm{p}_{\text{obs},i})=0,$

that lead us to

${\bm{c}}\cdot{\bm{\rho}}_{i}=({\bm{c}}_{\oplus}-{\bm{c}})\cdot\bm{q}_{i}-{\bm{c}}_{\oplus}\cdot\bm{p}_{\text{obs},i},$
(30)

for $i=1,2,3$. We also introduce the quantities

$\bm{\tau}=\frac{1}{\sqrt{p}}\text{adj}(R)\hat{\bm{c}},\qquad\bm{T}=\frac{1}{\sqrt{p_{\oplus}}}\text{adj}(Q_{\oplus})\hat{\bm{c}}_{\oplus},$
(31)

which are the same as in (6), and define the matrix

$\mathcal{C}=\left[\text{adj}(Q)-\text{adj}(Q_{\oplus})\right]\hat{P}.$ (32)

Like in the geocentric case, we have the relations

$\displaystyle\varepsilon_{ijk}\frac{\tau_{k}}{\kappa}\bm{c}\cdot\bm{\rho}_{i}$
$\displaystyle=\bm{q}_{i}\times\bm{q}_{j}\cdot\bm{\rho}_{i}-\bm{\rho}_{i}\times\bm{\rho}_{j}\cdot\bm{q}_{i},$
(33)
$\displaystyle\varepsilon_{ijk}\frac{T_{k}}{\kappa}\bm{c}_{\oplus}\cdot\bm{\rho}_{i}+{\cal
C}_{ki}\rho_{i}$ $\displaystyle=\bm{q}_{i}\times\bm{q}_{j}\cdot\bm{\rho}_{i},$
(34)

where $(i,j,k)$ is varied so that all the 6 permutations of the set
$\\{1,2,3\\}$ are considered, and $\mathcal{C}_{ki}$ denotes the element of
the $k$-th row and $i$-th column of the matrix $\mathcal{C}$.

Subtracting (33) from (34) we obtain

$\varepsilon_{ijk}\frac{1}{\kappa}(T_{k}{\bm{c}}_{\oplus}-\tau_{k}{\bm{c}})\cdot\bm{\rho}_{i}+{\cal
C}_{ki}\rho_{i}=\bm{\rho}_{i}\times\bm{\rho}_{j}\cdot\bm{q}_{i},$ (35)

and following the same procedure as in the geocentric case we can write (35)
as

$\varepsilon_{ijk}\left[\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\left(\hat{\bm{\rho}}_{i}+\alpha_{ik}\hat{\bm{q}}_{i}\right)+\frac{\kappa\mathcal{C}_{ki}}{T_{k}}-\frac{\alpha_{ik}}{q_{i}}\bm{c}_{\oplus}\cdot\bm{p}_{\text{obs},i}\right]T_{k}=\kappa\rho_{j}\hat{\bm{\rho}}_{i}\times\hat{\bm{\rho}}_{j}\cdot\bm{q}_{i},$
(36)

where the coefficients $\alpha_{ik}$ are defined as in (13), with the new
interpretation for $\rho_{i}$ and $q_{i}$.

Using the same approximations as in (15), we also note that

$\displaystyle\text{adj}(\hat{P})P\bm{\theta}$
$\displaystyle\simeq\text{adj}(\hat{P})P\bm{\tau}=\text{adj}(\hat{P})Q_{\oplus}(\bm{T}-\bm{\tau})-\text{adj}(\hat{P})P_{\text{obs}}\bm{\tau}$
(37)
$\displaystyle\simeq\frac{1}{6}\bigg{(}\frac{1}{r_{2}^{3}}-\frac{1}{q_{\oplus,2}^{3}}\bigg{)}\bm{u}-\text{adj}(\hat{P})P_{\text{obs}}\bm{\theta},$

where $q_{\oplus,2}=|\bm{q}_{\oplus,2}|$ and

$\bm{u}=(u_{1},u_{2},u_{3})^{t}=\text{adj}(\hat{P})Q_{\oplus}\bm{\theta}^{3}.$

The presence of the term $\text{adj}(\hat{P})P_{\text{obs}}\bm{\theta}$ in
(37) prevents us from making the same simplification that allowed to express
the $\alpha_{ik}$ as functions of known quantities. Noting that

$\text{adj}(\hat{P})P_{\text{obs}}=\mathcal{O}(p_{\text{obs}}),$

where $p_{\text{obs}}=|\bm{p}_{\text{obs}}|$, we neglect this term and obtain

$\det(\hat{P})\bm{\delta}\odot\bm{\theta}\simeq\frac{1}{6}\left(\frac{1}{r_{2}^{3}}-\frac{1}{q_{\oplus,2}^{3}}\right)\bm{u},$

with

$\bm{\delta}=(\rho_{1},\rho_{2},\rho_{3})^{t}.$

As a consequence, the expressions for $\alpha_{13}$ and $\alpha_{31}$ given in
(21) can still be used, with the new interpretation for the vectors
$\bm{\rho}$ and $\bm{q}$.

###  3.2 A linear equation involving ${\bm{c}}_{\oplus}-{\bm{c}}$

Choosing $(i,j,k)=(1,2,3),(3,2,1)$ in (36) and eliminating $\kappa\rho_{2}$
from the resulting equations, we obtain

$\displaystyle
a_{3}\left[\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\left(\hat{\bm{\rho}}_{3}+\alpha_{31}\hat{\bm{q}}_{3}\right)+\frac{\kappa\mathcal{C}_{13}}{T_{1}}-\frac{\alpha_{31}}{q_{3}}\bm{c}_{\oplus}\cdot\bm{p}_{\text{obs},3}\right]=$
(38) $\displaystyle
a_{1}\left[(\bm{c}_{\oplus}-\bm{c})\cdot\left(\hat{\bm{\rho}}_{1}+\alpha_{13}\hat{\bm{q}}_{1}\right)+\frac{\kappa\mathcal{C}_{31}}{T_{3}}-\frac{\alpha_{13}}{q_{1}}\bm{c}_{\oplus}\cdot\bm{p}_{\text{obs},1}\right],$

with

$a_{1}=\frac{\left(\hat{\bm{\rho}}_{2}\times\hat{\bm{\rho}}_{3}\cdot\bm{q}_{3}\right){q}_{2}}{T_{1}\sqrt{p_{\oplus}}}=\frac{\left(\hat{\bm{\rho}}_{2}\times\hat{\bm{\rho}}_{3}\cdot\bm{q}_{3}\right){q}_{2}}{\bm{q}_{\oplus,2}\times\bm{q}_{\oplus,3}\cdot\hat{\bm{c}}_{\oplus}},\quad
a_{3}=\frac{\left(\hat{\bm{\rho}}_{1}\times\hat{\bm{\rho}}_{2}\cdot\bm{q}_{1}\right){q}_{2}}{T_{3}\sqrt{p_{\oplus}}}=\frac{\left(\hat{\bm{\rho}}_{1}\times\hat{\bm{\rho}}_{2}\cdot\bm{q}_{1}\right){q}_{2}}{\bm{q}_{\oplus,1}\times\bm{q}_{\oplus,2}\cdot\hat{\bm{c}}_{\oplus}}.$
(39)

Like in the geocentric case, defining the vectors

$\bm{\gamma}=a_{1}\left(\hat{\bm{\rho}}_{1}+\alpha_{13}\hat{\bm{q}}_{1}\right),\qquad\bm{\varphi}=a_{3}\left(\hat{\bm{\rho}}_{3}+\alpha_{31}\hat{\bm{q}}_{3}\right),$

we can write equation (38) as

$(\bm{\gamma}-\bm{\varphi})\cdot\left(\bm{c}_{\oplus}-\bm{c}\right)=D,$ (40)

with

$\displaystyle D$
$\displaystyle=\kappa\left(\frac{a_{3}\mathcal{C}_{13}}{T_{1}}-\frac{a_{1}\mathcal{C}_{31}}{T_{3}}\right)+\left(\frac{a_{1}\alpha_{13}}{q_{1}}\bm{p}_{\text{obs},1}-\frac{a_{3}\alpha_{31}}{q_{3}}\bm{p}_{\text{obs},3}\right)\cdot\bm{c}_{\oplus}.$

######  Remark 4.

We can add to (40) two equations choosing the index pairs
$\\{(2,3,1),(1,3,2)\\}$ and $\\{(3,1,2),(2,1,3)\\}$ in (36). In this way, we
can write a linear system of three equations for ${\bm{c}}_{\oplus}-{\bm{c}}$
using only three observations. However, we expect that the matrix of this
system is ill-conditioned because $Q=Q_{\oplus}+\mathcal{O}(p_{\rm obs})$, so
that the vectors $\bm{q}_{1},\bm{q}_{2},\bm{q}_{3}$ are almost coplanar.

###  3.3 Equations of the topocentric method

Following Section 2.5, if we consider two different choices of the three
observations, we obtain the system

$\left\\{\begin{aligned}
&(\bm{\gamma}_{1}-\bm{\varphi}_{1})\cdot\left(\bm{c}_{\oplus}-\bm{c}\right)=D_{1}\\\
&(\bm{\gamma}_{2}-\bm{\varphi}_{2})\cdot\left(\bm{c}_{\oplus}-\bm{c}\right)=D_{2}.\end{aligned}\right.$
(41)

Set

$\bm{w}=(\bm{\gamma}_{1}-\bm{\varphi}_{1})\times(\bm{\gamma}_{2}-\bm{\varphi}_{2}),$

and assume $\bm{w}\neq\bm{0}$. Then the general solution of (41) takes the
form

$\bm{c}_{\oplus}-\bm{c}=\lambda\bm{w}+\bm{g},$ (42)

where $\lambda\in\mathbb{R}$, and $\bm{g}$ is a particular solution of (41),
e.g. the one fulfilling $\bm{g}\cdot\bm{w}=0$.

In order to constrain the values of $\lambda$ we proceed as follow. Note that
we can write (36) with $(i,j,k)=(1,2,3)$ as

$\frac{a_{3}c_{\oplus}\rho_{2}}{q_{2}}=\left(\bm{c}_{\oplus}-\bm{c}\right)\cdot\left(\hat{\bm{\rho}}_{1}+\alpha_{13}\hat{\bm{q}}_{1}\right)+\frac{\kappa\mathcal{C}_{31}}{T_{3}}-\frac{\alpha_{13}}{q_{1}}\bm{c}_{\oplus}\cdot\bm{p}_{\text{obs},1}.$
(43)

Inserting the general solution (42) into (43) and (30) with $i=2$ we obtain

${\rho}_{2}=\frac{1}{b}(\lambda\bm{w}+\bm{g})\cdot\bm{\gamma}_{1}+f,$ (44)

where

$b=\frac{a_{1}a_{3}c_{\oplus}}{q_{2}},\qquad
f=\frac{q_{2}}{a_{3}c_{\oplus}}\left(\frac{\kappa\mathcal{C}_{31}}{T_{3}}-\frac{\alpha_{13}}{q_{1}}\bm{c}_{\oplus}\cdot\bm{p}_{\text{obs},1}\right),$
(45)

and

$\left({\bm{c}}_{\oplus}-\lambda\bm{w}-\bm{g}\right)\cdot\hat{{\bm{\rho}}}_{2}\,\rho_{2}=(\lambda\bm{w}+\bm{g})\cdot\bm{q}_{2}-{\bm{c}}_{\oplus}\cdot\bm{p}_{\text{obs},2}.$
(46)

Substituting the expression (44) of $\rho_{2}$ in (46) we get

$\left({\bm{c}}_{\oplus}-\lambda\bm{w}-\bm{g}\right)\cdot\hat{{\bm{\rho}}}_{2}\left[\frac{1}{b}(\lambda\bm{w}+\bm{g})\cdot\bm{\gamma}_{1}+f\right]=(\lambda\bm{w}+\bm{g})\cdot\bm{q}_{2}-{\bm{c}}_{\oplus}\cdot\bm{p}_{\text{obs},2},$

which can be written as

$\begin{split}&(\bm{w}\cdot{\bm{\gamma}}_{1})(\bm{w}\cdot\hat{\bm{\rho}}_{2})\lambda^{2}+\bigl{[}(\bm{w}\cdot{\bm{q}}_{2})b-(\bm{w}\cdot\bm{\gamma}_{1})(\bm{c}_{\oplus}-\bm{g})\cdot\hat{\bm{\rho}}_{2}+(\bm{w}\cdot{\hat{\bm{\rho}}}_{2})(\bm{g}\cdot\bm{\gamma}_{1}+bf)\bigr{]}\lambda\\\\[4.30554pt]
&+b\left(\bm{g}\cdot{\bm{q}_{2}}-\bm{c}_{\oplus}\cdot\bm{p}_{\rm{obs},2}\right)-\left(\bm{g}\cdot\bm{\gamma}_{1}+bf\right)(\bm{c}_{\oplus}-\bm{g})\cdot\hat{\bm{\rho}}_{2}=0.\end{split}$
(47)

Equation (47) can be compared with equation (28). It is worth noting that in
the topocentric formulation we do not have the solution $\lambda=0$ as in
Mossotti’s original method, so that this formulation leads to a quadratic
equation.

######  Remark 5.

In the topocentric case we could add a third linear equation to system (41) by
choosing three different triplets of observations, out of the available four.
However, we expect that also in this case the system is ill-conditioned,
because the vectors $\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},\bm{q}_{4}$ are almost
coplanar.

##  4 Mossotti’s method for space debris

We can follow the same scheme introduced in Section 3 for the computation of
the orbits of space debris, assuming that the Earth is spherical and rotates
with uniform angular velocity. In this case we use the rescaled time

$\theta=\kappa_{\oplus}t,$

with $\kappa_{\oplus}=\sqrt{Gm_{\oplus}}$.

Figure 3: Geometry of observations of space debris.

Here the vector $\bm{q}$ represents the geocentric positions of the observer,
$\bm{r}$ and ${\bm{\rho}}=\bm{r}-\bm{q}$ give the geocentric and topocentric
positions of the debris, and ${\bm{c}}={\bm{r}}\times\dot{{\bm{r}}}$ is its
orbital angular momentum. We consider the orthogonal decomposition

$\bm{q}=\bm{q}^{\perp}+\bm{p}_{\rm obs},$

with

$\bm{q}^{\perp}=({\bm{e}}_{3}\times\bm{q})\times{\bm{e}}_{3},\qquad\bm{p}_{\rm
obs}=(\bm{q}\cdot{\bm{e}}_{3}){\bm{e}}_{3},$

where ${\bm{e}}_{3}$ is the unit vector of the Earth rotation axis, see Figure
3. Moreover, we introduce the vector

${\bm{c}}_{\rm
obs}=\bm{q}^{\perp}\times\frac{d}{dt}{\bm{q}}^{\perp}=\bm{q}^{\perp}\times\dot{\bm{q}},$

where the last equality holds because $\bm{p}_{\rm obs}$ is constant.

We can write a quadratic equation analogous to (47) simply by substituting the
vectors ${\bm{c}}_{\oplus}$, $\bm{q}_{\oplus}$ with ${\bm{c}}_{\rm obs}$,
$\bm{q}^{\perp}$, and the parameter $p_{\oplus}$ with $|\bm{q}^{\perp}|$.

##  5 Numerical tests

In this section we test the performance of Mossotti’s original geocentric
method (see Section 2) and its topocentric version introduced in Section 3. In
the following, we denote the former by Mgeo and the latter by Mtop.

In Sections 5.1, 5.2 we compare Mgeo with Mtop using simulated observations
(right ascension and declination) computed for the site of the Pan-STARRS1
telescope, mount Haleakala, Hawaii, without taking into account observability
conditions, i.e. the asteroids are not necessarily visible in the night sky.
Moreover, we assume that the four observations given in input to Mgeo and Mtop
are equally spaced in time. The time interval $\Delta t$ between two
consecutive observations is varied in the two intervals:

$I_{1}=[15,200]\text{ minutes,}\qquad I_{2}=[0.25,100]\text{ days.}$ (48)

The comparison is based on the computation of the angular momentum vector
${\bf c}$ and the following related quantities: its magnitude $c=|{\bf c}|$
and direction $\hat{{\bf c}}={\bf c}/c$, the longitude of the (ascending) node
$\Omega$, and the inclination $i$. We denote by $x_{t}$ the _true_ value of
the quantity $x$, and by $x_{geo}$, $x_{top}$ the values of $x$ computed by
Mgeo, Mtop, respectively. Note that while Mgeo always produces one solution
for ${\bf c}$, the method Mtop can give two solutions. If this is the case, we
select the one for which $|{\bf c}_{t}-{\bf c}_{top}|$ is smaller.

In Section 5.3 the method Mtop is compared to Gauss’s method for initial orbit
determination. Synthetic data have been obtained that take into account the
observability conditions and the expected real cadence of the observations
from the Vera C. Rubin Observatory, which is currently under construction in
Chile. With respect to the previous tests, we remark that for any set of four
observations the time interval $\Delta t$ is not constant.

###  5.1 Tests without astrometric errors

We consider simulated observations of the asteroid Vesta without astrometric
error. The time $\Delta t$ between two consecutive observations is varied in
the two intervals $I_{1}$, $I_{2}$ specified in (48). Given $\Delta t$, we
select 105 different initial epochs in a random way, and for each of them we
generate four observations. For this purpose, we use the orbit of Vesta at the
epoch 59200 MJD from the AstDyS-2
website444https://newton.spacedys.com/astdys/, last access February 13, 2021.
and propagate it to the desired epochs assuming Keplerian motion. The angular
momentum vector defined by the Keplerian orbit is assumed to be the true
solution (${\bf c}_{t}$). Finally, we apply Mgeo, Mtop to each set of four
observations.

Let us first consider the case of short arcs of observation. For each $\Delta
t$ selected in $I_{1}$ we compute the differences between the true values of
the inclination ($i_{t}$), longitude of the node ($\Omega_{t}$), magnitude of
the angular momentum vector ($c_{t}$) of Vesta, and the values obtained by
either Mgeo or Mtop. Some relevant statistical quantities related to these
errors are shown in Figure 4 as functions of $\Delta t$. We note that the
topocentric version of Mossotti’s method provides much better results than the
original method. In particular, we observe that the performance of Mtop
improves as $\Delta t$ increases, and it stabilizes when the time interval is
about 30 minutes. It is remarkable that for $\Delta t$ larger than 25 minutes
the error in the inclination is smaller than 0.01 degrees for the solutions
that fall within the 1st and 3rd quartile and it is smaller than 0.1 degrees
for the solutions that fall within the 5th and 95th percentile.

Figure 4: Statistics related to the errors in the inclination $i$, longitude
of the node $\Omega$, and magnitude of the angular momentum vector $c$,
obtained by Mgeo (_left_) and Mtop (_right_) as functions of $\Delta t$ in the
case of observations of Vesta without astrometric error. In particular, we
show the median (black line), the 1st and 3rd quartiles (blue and red lines),
and the 5th and 95th percentiles (gray lines). The angles are in degrees and
the $\Delta t$ step is 1.5 seconds. Figure 5: Statistics related to the errors
in the angular momentum vector ${\bf c}$ and its direction $\hat{\bf c}$,
obtained by Mgeo and Mtop as functions of $\Delta t$ in the case of
observations of Vesta without astrometric error. In particular, we show the
median (Q2), the 1st and 3rd quartiles (Q1, Q3), and the 95th percentile
(P95). Here the $\Delta t$ step is 15 minutes.

We then allow $\Delta t$ to take values in the interval $I_{2}$, which
corresponds to wider arcs of observation. For each $\Delta t$ selected in
$I_{2}$ we compute the errors in the angular momentum vector and its
direction. Statistical quantities related to these errors are shown in Figure
5 as functions of $\Delta t$. Mtop continues to show better results and a
smoother behavior for values of $\Delta t$ smaller than 30 days, even if the
improvement over Mgeo is less pronounced with respect to that shown in Figure
4. In such interval, the geocentric method is much more sensitive to $\Delta
t$ and large oscillations having a period of one day appear. This is due to
not accounting for the topocentric position of the observer. We also notice
that both methods exhibit an optimal performance for $\Delta t\approx$ 3
weeks, and Mtop obtains in $75\%$ of the solutions an error smaller than
$0.2\%$ and $0.03\%$ in the vectors ${\bf c}$ and $\hat{{\bf c}}$,
respectively.

If we consider time intervals longer than 30 days, the performance of the two
methods is almost comparable, which is expected because the terms introduced
in the topocentric version become smaller as $\Delta t$ grows. The solutions
with both methods deteriorate for $\Delta t>80$ days since they rely on
Taylor’s expansions with respect to $\Delta t$.

${\Delta t}$ (days) | 0.02 | 1 | 10 | 50 | 100 | True value
---|---|---|---|---|---|---
$i_{1}$ (geo) | 17.80160 | 7.06410 | 7.20588 | 7.09387 | 7.05670 | 7.14165
$i_{2}$ (geo) | 0.00315 | 0.00297 | 0.00303 | 0.00436 | 0.00278
$i_{1}$ (top) | 7.14611 | 7.06413 | 7.20572 | 7.09407 | 7.05700
$i_{2}$ (top) | 17.79875 | 0.00156 | 0.00328 | 0.00305 | 0.00881
$\Omega_{1}$ (geo) | 150.28159 | 103.18144 | 104.36892 | 103.20691 | 105.07597 | 103.80838
$\Omega_{2}$ (geo) | 140.61411 | 146.13681 |  $-$143.70328 | 160.84846 | 179.83017
$\Omega_{1}$ (top) | 103.83824 | 103.17391 | 104.35977 | 103.19931 | 105.07223
$\Omega_{2}$ (top) | 150.27855 | 167.05122 |  $-$112.30086 |  $-$174.23000 |  $-$66.15466
$c_{1}$ (geo) | 1125.88148 | 0.02600 | 0.02650 | 0.02694 | 0.03463 | 0.02633
$c_{2}$ (geo) | 0.01721 | 0.01721 | 0.01719 | 0.01720 | 0.01719
$c_{1}$ (top) | 0.02628 | 0.02600 | 0.02650 | 0.02694 | 0.03463
$c_{2}$ (top) | 1124.55534 | 0.01720 | 0.01719 | 0.01720 | 0.01718
Table 1: Values of the inclination ($i$), longitude of the node ($\Omega$),
magnitude of the angular momentum vector ($c$), obtained by Mtop and Mgeo,
with observations of Vesta not affected by astrometric error and with
different time intervals between consecutive observations. The epoch of the
true orbit is 59200 MJD, the epoch of the first observation is $\sim$
59200.012 MJD (due to aberration correction), and the angles are in degrees.

In Table 1 we report the inclination, longitude of the node, and magnitude of
the angular momentum vector obtained by Mgeo and Mtop from observations of
Vesta without astrometric error, considering different values of $\Delta t$
with the same epoch for the first observation. We also show for Mgeo with the
label “2” the solution of Mossotti’s original method which is always
discarded, i.e. the one for which the angular momenta of the asteroid and the
Earth are equal. We observe that if $\Delta t$ is large enough, the values of
$i$ and $c$ of this spurious solution are close to the _wrong_ solution from
Mtop (also labeled by “2”). The same behavior is not observed for $\Omega$
because the Earth orbital inclination is small and therefore small variations
in $i$ can cause large deviations in $\Omega$. For $\Delta t=30$ minutes
($0.02$ days) the solution from Mgeo with label “1” is very close to the wrong
solution from Mtop: in this case only Mtop gives values of $i$, $\Omega$, $c$
close to the true ones.

###  5.2 Tests with astrometric errors

The same numerical tests described in the previous section for the asteroid
Vesta are carried out by introducing an astrometric error with zero mean and
standard deviation (rms) of 0.1 arcsec in the values of right ascension and
declination, which is typical of modern asteroid surveys like Pan-STARRS1.
Neither Mgeo nor Mtop gives reliable results for $\Delta t$ between 15 and 200
minutes. Indeed, determining a preliminary orbit from a single short arc is a
challenging task, sometimes impossible without considering infinitely many
solutions [16], [15]. However, the computation of a preliminary orbit can be
performed by linking together two or more short arcs (e.g. [8], [9]). On the
other hand, both methods yield satisfactory results for time intervals between
two consecutive observations larger than 8 hours with only a slight
degradation of their performance with the introduction of astrometric error
(compare Figures 5 and 6). As in the case without astrometric error, Mtop is
better than Mgeo for any considered $\Delta t$. A smoother behavior of Mtop is
observed for $\Delta t$ smaller than 30 days. Both methods perform best for
$\Delta t\approx$ 3 weeks.

Figure 6: Same as Figure 5 but the observations of Vesta are affected by
astrometric error with zero mean and rms of $0.1$ arcsec. Here the $\Delta t$
step is 15 minutes.

In our third test we use Mgeo and Mtop with synthetic observations of the
546077 numbered asteroids known to the date of October 19, 2020. Their orbital
elements at the epoch 59200 MJD (from AstDyS-2) provide the true orbit and are
used to simulate observations from the Pan-STARRS1 telescope site. The epoch
of the first observation is obtained from 59200 MJD considering aberration
correction; then the subsequent three observations are simulated by Keplerian
propagation. The time $\Delta t$ is varied in the interval $I_{2}$ given in
(48), and the standard deviation of the astrometric error spans from 0 to 1
arcsec. Figures 7, 8 show that for time intervals between 20 and 40 days both
Mgeo and Mtop give good results. A closer look at this range of values of
$\Delta t$ for Mgeo reveal the same oscillations displayed in Figures 5 and 6.
On the other hand, Mtop smooths out such oscillations. The best performance of
both Mgeo and Mtop is reached again for $\Delta t\approx$ 3 weeks.

Figure 7: Performance of Mossotti’s original method (_left_) and its
topocentric version (_right_) with observations of 546077 numbered asteroids.
For each asteroid we compute the difference between the true angular momentum
vector ${\bf c}_{t}$ at the epoch 59200 MJD, and the same quantity ${\bf c}$
obtained by either Mgeo or Mtop on a uniform grid of values of $\Delta t$ and
rms of the astrometric error in the intervals $[0.25,\,100]$ days and
$[0,\,1]$ arcsec. Colored representations of the median and the 95th
percentile related to the error $|{\bf c}_{t}-{\bf c}|/c_{t}$ are displayed.
The step of $\Delta t$ is 6 hours, the step of rms is $0.002$ arcsec.

Figure 8: Same as Figure 7 for the error in the direction of the angular
momentum vector.

###  5.3 Tests on synthetic survey data without astrometric errors

Our final test introduces more realism into the simulations in an attempt to
assess the performance of the topocentric version of Mossotti’s method (Mtop)
on synthetic data that uses a realistic cadence and accounts for actual
observability of the asteroids. With the motivation that Mossotti’s initial
orbit determination method could be applied to linking observations of unknown
asteroids in contemporary and future asteroid surveys, we applied Mtop to
synthetic observations from the Vera Rubin Observatory’s (VRO) Legacy Survey
of Space and Time (LSST) [11]. They have developed a high-fidelity survey
scheduler that will be employed in final operations but is currently being
used to simulate and optimize the survey strategy [2, 3, 18]. We used a single
survey simulation for one month of surveying that did not include any
astrometric error. Then we extracted the first four synthetic detections of
all the detected numbered NEOs, Trojans, Centaurs, and TNOs, but only a small
subset of the detected main belt objects so that they would not dominate our
results. The epochs of observation were all within about 30 days and the
astrometric positions were generated with a full $n$-body integration. This
process produced a set of four detections of 1535 objects distributed
throughout the solar system. We then processed all the detections with both
Mtop and our implementation of Gauss’s method, by limiting our search to
bounded orbits only.

In our sample Gauss’s method was able to produce orbits for 1493 objects (i.e.
$\sim 97$%) and Mtop provided solutions for 1395 objects (i.e. $\sim 91$%).
However, we had 59 occurrences of a negative discriminant of equation (47)
with Mtop, and for some of them we were able to recover an acceptable orbit by
setting the discriminant equal to zero. Doing so increases the number of
solutions for Mtop to 1454 (i.e. $\sim 95$%).

While Gauss’s method can yield three different solutions for the same set of
observations and Mtop can yield two, on this set of data they had multiple
solutions for about 50% and 45% of the objects, respectively. Gauss’s method
did not produce any orbit for 42 objects and Mtop for 81. Moreover, both of
them failed in 30 of these cases. In 12 cases Mtop was able to obtain at least
one orbit while Gauss’s method was not, and for some cases it found an
acceptable orbit.

The primary benefit of Mtop is that in our limited testing it appears to be
about $6$ times faster than Gauss’s method. We repeated the orbit computation
1000 times for each of the 1535 objects using an Intel Xeon processor, with
base clock 3.30 GHz: Gauss’s algorithm took $\sim 77$ seconds, while Mtop
$\sim 13$ seconds.

Gauss’s technique provides better solutions for objects throughout the solar
system (Figure 9 and Table 2). When comparing the derived orbital elements
with their actual values we used the derived orbit solution that had the
lowest 5-element $D$-criterion555The $D$-criterion quantifies the difference
between two orbits using all the orbital elements except for the mean anomaly.
relative to the actual orbit [4]. A simple visual comparison of the results
suggests that Gauss’s method is more likely to produce good orbital elements
and less likely to yield wildly different values. Quantitative orbital element
comparisons confirm this impression (Table 2).

Figure 9: All the panels present the derived orbital parameter on the $y$-axis versus the actual parameter value on the $x$-axis using error-free data from a VRO-LSST simulation. The top panels are for Gauss’s method and the bottom panels use the topocentric version of Mossotti’s method. From left to right the panels compare $\log_{10}(a/[au])$, where $a$ is the semi-major axis, eccentricity, and inclination. The red reference line in each figure has a slope of 1. | Gauss | Mossotti
---|---|---
orbital element | $\bar{\Delta}$ | rms | $\bar{\Delta}$ | rms
semi-major axis (au) | 0.13 | 3.44 | 0.28 | 3.18
eccentricity | 0.01 | 0.07 | 0.03 | 0.12
inclination (deg) |  $-$0.02 | 0.27 | 0.05 | 2.31
long. node (deg) |  $-$0.09 | 1.83 | 0.10 | 4.58
arg. perihelion (deg) |  $-$0.38 | 8.51 | 0.10 | 19.17
Table 2: Mean and rms of the difference between the actual and derived orbital
elements for both Gauss’s method and the topocentric version of Mossotti’s
method. For the angular elements we eliminated some outliers to better
illustrate the difference between the methods in the vast majority of cases.

##  6 Conclusions

In this paper we have revisited Mossotti’s orbit determination method working
with four geocentric observations of a celestial body, and extended it to the
case of topocentric observations. While Mossotti’s method yields linear
equations for the components of the angular momentum vector, the topocentric
version leads to a quadratic equation. Numerical simulations with synthetic
observations both without and with astrometric error show that the topocentric
method improves the original one. Considering all the numbered asteroids, and
generating for each of them four observations equally spaced in time, we find
that both these methods show an optimal behavior for a time separation $\Delta
t$ between two consecutive observations of about 3 weeks. Finally, we compare
the new method with Gauss’s method using synthetic observations without
astrometrical error that reproduce the expected scheduling of the Vera Rubin
Observatory’s (VRO) Legacy Survey of Space and Time (LSST), characterized by
an average $\Delta t$ of about 4 days. Gauss’s method provides good orbits for
a larger number of objects than the topocentric version of Mossotti’s method,
which, on the other hand, is faster.

##  7 Acknowledgments

We thank Dr. Lynne Jones, Dr. Siegfried Eggl, Dr. Sam Cornwall, and Dr. Mario
Jurić of the University of Washington (WA) for assistance in identifying,
accessing, and understanding the appropriate VRO/LSST simulations. We also
thank the anonymous referees for their useful comments. GFG and GB acknowledge
the project MIUR-PRIN 20178CJA2B titled “New frontiers of Celestial Mechanics:
theory and applications”. GFG, GB and OR have been partially supported by the
MSCA-ITN Stardust-R, Grant Agreement n. 813644 under the H2020 research and
innovation program.

##  Appendix A Comparison with Mossotti’s original paper

Mossotti |  here |  |  Mossotti |  here
---|---|---|---|---
$(x,\,y,\,z)^{t}$ |  $\bm{r}$ |  |  $(Q^{\prime\prime\prime},\,-Q^{\prime\prime},\,Q^{\prime})^{t}$ |  $\displaystyle\frac{1}{\theta_{12}^{2}}\Bigl{(}\frac{u_{1}}{\theta_{23}},\frac{u_{2}}{\theta_{31}},\frac{u_{3}}{\theta_{12}}\Bigr{)}^{t}$
$(X,\,Y,\,Z)^{t}$ |  $\bm{q}$ |  |  $\left(\begin{matrix}B_{1}^{\prime\prime\prime}&B_{2}^{\prime\prime\prime}&B_{3}^{\prime\prime\prime}\\\ -B_{1}^{\prime\prime}&-B_{2}^{\prime\prime}&-B_{3}^{\prime\prime}\\\ B_{1}^{\prime}&B_{2}^{\prime}&B_{3}^{\prime}\\\ \end{matrix}\right)$ |
$(\mu,\,\nu,\,\omega)^{t}$ |  $\hat{\bm{q}}$ |  |  |  $\text{adj}(\hat{Q})\hat{P}$
$R,\,D$ |  $q$ |  |  |
$(m,\,n,\,o)^{t}$ |  $\hat{\bm{\rho}}$ |  |  $\hskip 7.11317pt\left(\begin{matrix}b_{1}^{\prime\prime\prime}&b_{2}^{\prime\prime\prime}&b_{3}^{\prime\prime\prime}\\\ -b_{1}^{\prime\prime}&-b_{2}^{\prime\prime}&-b_{3}^{\prime\prime}\\\ b_{1}^{\prime}&b_{2}^{\prime}&b_{3}^{\prime}\\\ \end{matrix}\right)$ |
$\delta$ |  $\rho$ |  |  |  $\text{adj}(\hat{P})\hat{Q}$
$\kappa(C^{\prime\prime\prime},\,-C^{\prime\prime},\,C^{\prime})^{t}$ |  $\bm{c}_{\oplus}$ |  |  |
$\kappa(c^{\prime\prime\prime},\,-c^{\prime\prime},\,c^{\prime})^{t}$ |  $\bm{c}$ |  |  $(\psi^{\prime},\,\psi^{\prime\prime},\,\psi^{\prime\prime\prime})$ |  $\frac{1}{\kappa}\bm{c}^{t}\hat{P}$
$(T^{\prime\prime\prime},\,-T^{\prime\prime},\,T^{\prime})^{t}$ |  $\bm{T}$ |  |  $(\chi^{\prime},\,\chi^{\prime\prime},\,\chi^{\prime\prime\prime})$ |  $\frac{1}{\kappa}\bm{c}_{\oplus}^{t}\hat{P}$
$(\tau^{\prime\prime\prime},\,-\tau^{\prime\prime},\,\tau^{\prime})^{t}$ |  $\bm{\tau}$ |  |  $(\varphi^{\prime},\,\varphi^{\prime\prime},\,\varphi^{\prime\prime\prime})$ |  $\frac{1}{\kappa}(\bm{c}_{\oplus}-\bm{c})^{t}\hat{Q}$
$(\theta^{\prime\prime\prime},\,-\theta^{\prime\prime},\,\theta^{\prime})^{t}$ |  $\bm{\theta}$ |  |  $(\Phi^{\prime},\,\Phi^{\prime\prime},\,\Phi^{\prime\prime\prime})$ |  $\frac{1}{\kappa}(\bm{c}_{\oplus}-\bm{c})^{t}\hat{P}$
|  |  |  |
Table 3: Correspondence table between Mossotti’s and our notation.

## References

  * [1] A. Celletti and G. Pinzari.  Four classical methods for determining planetary elliptic elements: a comparison.  Celestial Mechanics and Dynamical Astronomy, 93:1–52, 2005.
  * [2] A. J. Connolly et al.  An end-to-end simulation framework for the Large Synoptic Survey Telescope.  In George Z. Angeli and Philippe Dierickx, editors, Modeling, Systems Engineering, and Project Management for Astronomy VI, volume 9150 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, page 915014, August 2014.
  * [3] F. Delgado and M. A. Reuter.  The LSST Scheduler from design to construction.  In Alison B. Peck, Robert L. Seaman, and Chris R. Benn, editors, Observatory Operations: Strategies, Processes, and Systems VI, volume 9910 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, page 991013, July 2016.
  * [4] J. D. Drummond.  A test of comet and meteor shower associations.  Icarus, 45(3):545–553, March 1981.
  * [5] C. F. Gauss.  Theoria motus corporum in sectionibus conicis solem ambientium.  Reprinted by Dover publications in 1963, 1809.
  * [6] C. F. Gauss.  Werke, volume VI.  Available from Gallica, 1874.
  * [7] G. F. Gronchi.  Multiple solutions in preliminary orbit determination from three observations.  Celestial Mechanics and Dynamical Astronomy, 103/4:301–326, 2009\.
  * [8] G. F. Gronchi, G. Baù, and S. Marò.  Orbit determination with the two-body integrals. III.  Celestial Mechanics and Dynamical Astronomy, 123/2:105–122, 2015\.
  * [9] 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/2:211–232, 2017\.
  * [10] S. Herrick.  Astrodynamics. Vol. 1. Van Nostrand Reinhold, 1976.
  * [11] Ž. Ivezić et al.  LSST: From Science Drivers to Reference Design and Anticipated Data Products.  The Astrophysical Journal, 873(2), 2019.
  * [12] 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.
  * [13] 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.
  * [14] A. Milani and G. F. Gronchi.  Theory of Orbit Determination.  Cambridge Univ. Press, 2010.
  * [15] A. Milani, G. F. Gronchi, M. De’michieli Vitturi, and Z. Knežević.  Orbit determination with very short arcs. I admissible regions.  Celestial Mechanics and Dynamical Astronomy, 90(1-2):57–85, September 2004.
  * [16] A. Milani, M. E. Sansaturio, G. Tommei, O. Arratia, and S. R. Chesley.  Multiple solutions for asteroid orbits: Computational procedure and applications.  Astronomy & Astrophysics, 431:729–746, February 2005.
  * [17] O. F. Mossotti.  Nuova analisi del problema di determinare le orbite dei corpi celesti (1816-1818).  Domus Galileana, Pisa, 1942.  http://www.mat.uniroma2.it/%7Esimca/Testi/
mossotti-nuo-ana.pdf.

  * [18] E. Naghib, P. Yoachim, R. J. Vanderbei, A. J. Connolly, and R. L. Jones.  A Framework for Telescope Schedulers: With Applications to the Large Synoptic Survey Telescope.  The Astronomical Journal, 157(4), 2019.