Generalization of a method by Mossotti for initial orbit determination

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.


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.

The original method
Mossotti's method [17] leads to a set of linear equations for the components of c ⊕ − c, where c ⊕ and 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.

Units and preliminary definitions
In order to simplify the notation we use the rescaled time θ, defined by where m ⊕ and m are the masses of the Earth and the Sun, and g = Gm , where G is Newton's gravitational constant. We also use κ to denote Gauss's constant √ g. In the following we assume that θ ≈ κt, neglecting the constant m ⊕ /m . Take three of the four observations, at epochs t 1 < t 2 < t 3 , and set where the superscript t stands for transposition. We write c and c ⊕ for the orbital angular momenta of the asteroid and the Earth, respectively, and introduce their unit vectorŝ where c = |c|, c ⊕ = |c ⊕ |, being |x| the Euclidean norm of a vector x.
We also write r i and q i for the heliocentric positions of the asteroid and the Earth at the three epochs t i (i = 1, 2, 3), use ρ i = r i − q i for the geocentric position of the asteroid, and introduce the unit vectorsq Finally, we denote the parameters 1 of the orbits of the asteroid and the Earth by p, p ⊕ . They are defined by

Geometric relations
With the purpose of writing Mossotti's equations in a compact form, let us introduce the matrices where (x 1 | x 2 | x 3 ) is the matrix whose columns are the vectors x j . The corresponding adjugate matrices are We recall the following property, which holds for any square matrix M : where I is the identity matrix. The rank of the matrices Q and R is 2, since each triplet {q 1 , q 2 , q 3 } and {r 1 , r 2 , r 3 } is made by coplanar vectors and we assume that see Figure 1. This implies Since the angular momenta c, c ⊕ are respectively orthogonal to the orbital planes of the asteroid and the Earth, we have for i = 1, 2, 3. Then, we introduce the vectors Sun Earth's plane of motion Asteroid's plane of motion Figure 1: Geometry of the three observations.
Moreover, recalling that R = Q + P , we get Since 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 By (8), (9) and the definitions of the τ k , T k given in (5) we obtain where ε 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 Writing and using (4), we have Therefore, introducing the coefficients and simplifying ρ i we can write relations (12) as

Combining geometry of observations with two-body dynamics
Since the time intervals θ 12 , θ 23 are small compared to the orbital period, 2 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 θ 23 , θ 31 , θ 12 centering the expansion at the intermediate epoch.

Mossotti's equations for c ⊕ − 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 where the subscripts 1, 2 of γ and ϕ refer to the two triplets of observations. Set and assume w = 0. Then the general solution of (23) has the form with λ ∈ R, giving the direction of c ⊕ − c.
Remark 3. If γ 1 −ϕ 1 and γ 2 −ϕ 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 λ we proceed as follow.
Inserting the expression (24) of the general solution in (25) and (4) with i = 2 we obtain, respectively, and Substituting the expression (26) of ρ 2 in (27) yields a quadratic equation in λ, which is here the only unknown: This equation can be written as Substituting these expressions in (24)

Topocentric method
We first introduce some notation. Let us define q ⊕ , p 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 where ρ, ρ geo are the topocentric and geocentric positions of the asteroid, respectively (see Figure 2).

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 ρ and q: here ρ represents the topocentric position of the asteroid, and q gives the heliocentric position of the observer. Moreover, we introduce the matrices We recall the geometrical relations for i = 1, 2, 3. We also introduce the quantities which are the same as in (6), and define the matrix Like in the geocentric case, we have the relations where (i, j, k) is varied so that all the 6 permutations of the set {1, 2, 3} are considered, and C ki denotes the element of the k-th row and i-th column of the matrix C. Subtracting (33) from (34) we obtain and following the same procedure as in the geocentric case we can write (35) as where the coefficients α ik are defined as in (13), with the new interpretation for ρ i and q i . Using the same approximations as in (15), we also note that where q ⊕,2 = |q ⊕,2 | and The presence of the term adj(P )P obs θ in (37) prevents us from making the same simplification that allowed to express the α ik as functions of known quantities. Noting that adj(P )P obs = O(p obs ), where p obs = |p obs |, we neglect this term and obtain det(P )δ θ 1 6 As a consequence, the expressions for α 13 and α 31 given in (21) can still be used, with the new interpretation for the vectors ρ and q.

A linear equation involving c ⊕ − c
Choosing (i, j, k) = (1, 2, 3), (3, 2, 1) in (36) and eliminating κρ 2 from the resulting equations, we obtain with Like in the geocentric case, defining the vectors we can write equation (38) as with In this way, we can write a linear system of three equations for c ⊕ − c using only three observations. However, we expect that the matrix of this system is ill-conditioned because Q = Q ⊕ + O(p obs ), so that the vectors q 1 , q 2 , q 3 are almost coplanar.

Equations of the topocentric method
Following Section 2.5, if we consider two different choices of the three observations, we obtain the system Set and assume w = 0. Then the general solution of (41) takes the form where λ ∈ R, and g is a particular solution of (41), e.g. the one fulfilling g · w = 0. In order to constrain the values of λ we proceed as follow. Note that we can write (36) with (i, j, k) = (1, 2, 3) as Inserting the general solution (42) into (43) and (30) with i = 2 we obtain and (c ⊕ − λw − g) ·ρ 2 ρ 2 = (λw + g) · q 2 − c ⊕ · p obs,2 .
which can be written as Equation (47) can be compared with equation (28). It is worth noting that in the topocentric formulation we do not have the solution λ = 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 q 1 , q 2 , q 3 , q 4 are almost coplanar.

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 θ = κ ⊕ t, with κ ⊕ = √ Gm ⊕ . Here the vector q represents the geocentric positions of the observer, r and ρ = r − q give the geocentric and topocentric positions of the debris, and c = r ×ṙ is its orbital angular momentum. We consider the orthogonal decomposition where e 3 is the unit vector of the Earth rotation axis, see Figure 3. Moreover, we introduce the vector where the last equality holds because p obs is constant. We can write a quadratic equation analogous to (47) simply by substituting the vectors c ⊕ , q ⊕ with c obs , q ⊥ , and the parameter p ⊕ with |q ⊥ |.

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 M geo and the latter by M top .
In Sections 5.1, 5.2 we compare M geo with M top 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 M geo and M top are equally spaced in time. The time interval ∆t between two consecutive observations is varied in the two intervals: The comparison is based on the computation of the angular momentum vector c and the following related quantities: its magnitude c = |c| and directionĉ = c/c, the longitude of the (ascending) node Ω, 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 M geo , M top , respectively. Note that while M geo always produces one solution for c, the method M top can give two solutions. If this is the case, we select the one for which |c t − c top | is smaller. In Section 5.3 the method M top 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 ∆t is not constant.

Tests without astrometric errors
We consider simulated observations of the asteroid Vesta without astrometric error. The time ∆t between two consecutive observations is varied in the two intervals I 1 , I 2 specified in (48). Given ∆t, we select 10 5 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 website 4 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 (c t ). Finally, we apply M geo , M top to each set of four observations. Let us first consider the case of short arcs of observation. For each ∆t selected in I 1 we compute the differences between the true values of the inclination (i t ), longitude of the node (Ω t ), magnitude of the angular momentum vector (c t ) of Vesta, and the values obtained by either M geo or M top . Some relevant statistical quantities related to these errors are shown in Figure 4 as functions of ∆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 M top improves as ∆t increases, and it stabilizes when the time interval is about 30 minutes. It is remarkable that for ∆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. We then allow ∆t to take values in the interval I 2 , which corresponds to wider arcs of observation. For each ∆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 ∆t. M top continues to show better results and a smoother behavior for values of ∆t smaller than 30 days, even if the improvement over M geo is less pronounced with respect to that shown in Figure 4. In such interval, the geocentric method is much more sensitive to ∆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 ∆t ≈ 3 weeks, and M top obtains in 75% of the solutions an error smaller than 0.2% and 0.03% in the vectors c andĉ, 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 ∆t grows. The solutions with both methods deteriorate for ∆t > 80 days since they rely on Taylor's expansions with respect to ∆t.
In Table 1 we report the inclination, longitude of the node, and magnitude of the angular momentum vector obtained by M geo and M top from observations of Vesta without astrometric error, considering different values of ∆t with the same epoch for the first observation. We also show for M geo 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 ∆t is large enough, the values of i and c of this spurious solution are close to the wrong solution from M top (also labeled by "2"). The same behavior is not observed for Ω because the Earth orbital inclination is small and therefore small variations in i can cause large deviations in Ω. For ∆t = 30 minutes (0.02 days) the solution from M geo with label "1" is very close to the wrong solution from M top : in this case only M top gives values of i, Ω, c close to the true ones.

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.   arcsec in the values of right ascension and declination, which is typical of modern asteroid surveys like Pan-STARRS1. Neither M geo nor M top gives reliable results for ∆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, M top is better than M geo for any considered ∆t. A smoother behavior of M top is observed for ∆t smaller than 30 days. Both methods perform best for ∆t ≈ 3 weeks. In

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 (M top ) 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 M top 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 M top 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. ∼ 97%) and M top provided solutions for 1395 objects (i.e. ∼ 91%). However, we had 59 occurrences of a negative discriminant of equation (47) with M top , and for some of them we were able to While Gauss's method can yield three different solutions for the same set of observations and M top 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 M top for 81. Moreover, both of them failed in 30 of these cases. In 12 cases M top 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 M top 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 ∼ 77 seconds, while M top ∼ 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-criterion 5 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).

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  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.
methods show an optimal behavior for a time separation ∆t between two consecutive obser-