Determination of Orbital Elements and Ephemerides using the Geocentric Laplace's Method

Abstract

This paper presents a methodology for Initial Orbit Determination (IOD) based on a modification of the Laplace's geocentric method. The orbital elements for Near-Earth asteroids (1864) Daedalus, 2003 GW, 2019 JA8, a Hungaria-type asteroid (4690) Strasbourg, and the asteroids of the Main Belt (1738) Oosterhoff, (2717) Tellervo, (1568) Aisleen and (2235) Vittore were calculated. Input data observations from the Minor Planet Center MPC database and Astronomical Observatory of the Technological University of Pereira (OAUTP; MPC code W63) were used. These observations cover observation arcs of less than 22 days. The orbital errors, in terms of shape and orientation for the estimated orbits of the asteroids, were calculated. The shape error was less than 53 × 10^-3 AU, except for the asteroid 2019 JA8. On the other hand, errors in orientation were less than 0.1 rad, except for (4690) Strasbourg. Additionally, we estimated ephemerides for all bodies for up to two months. When compared with actual ephemerides, the errors found allowed us to conclude that these bodies can be recovered in a field of vision of 95' × 72' (OAUTP field). This shows that Laplace's method, though simple, may still be useful in the IOD study, especially for observatories that initiate programs of minor bodies observation.

1. INTRODUCTION

Progress in the acquisition of astronomical images has dramatically increased the available astrometric information of minor bodies of the Solar System. The way to calculate orbits in the last decades has changed (Bowell 1999), so now it is possible to use all this information. This process involves three fundamental stages: an Initial Orbit Determination (IOD) (Marsden 1985), a least squares fitting (Vallado & McClain 2001) and quality control (Milani et al. 2008). Additionally, there are challenges raised by the appearance of digital astrometry. Finding a solution to this problem is of great importance for humanity, in order to discover and characterize, in the least amount of time, minor bodies that imply an imminent risk of collision with the Earth (Aristova et al. 2018).

The pioneers to offer a solution to the problem of IOD were Laplace (Curtis 2014) and Gauss (Gauss 1809; Gauss & Davis 2004; Hwang 2009); however, these methods were conceived just over two centuries ago. Because the acquisition of an astronomical image took a great deal of time, these classical methods were elaborated for that context. At present, a large amount of data can be obtained on an observation night, so these methods must be reconsidered (Klokacheva 1991; Sokolskaya 1997) or replaced by alternate ones. For that purpose, comprehensive reviews of the Gauss & Laplace methods have been performed (Milani et al. 2008; Mirtorabi 2014), or completely alternative methods have been investigated, to calculate a finite set of preliminary orbits for solar system bodies, using the two-body integrals (Gronchi 2004; Gronchi et al. 2010).

These methods are adjusted more adequately to the large amount of astrometric information currently available; namely, they use all possible data (unlike the Gauss method that uses just three (Danby 1962). The use of all possible data permits astronomers to perform a good initial orbit es-timation for the celestial body, even when the observational data covers a very short section of the orbit, that is VSA (a very short arc). This definition is used when the sequence of observations belongs to the same object, or if it can be adjusted to a smooth curve (usually a low-degree polynomial Milani et al. 2005).

On the other hand, in studies such as Quijano et al. (2010), orbital elements of the asteroid 2003QO104 were determined using, principally, a linear interpolation by Lagrange method. However, this methodology does not calculate the temporal derivatives of both coordinates. In other words, an attributable is not computed, and therefore it is not possible to determine either the geodesic curvature or the proper motion, which are essential constants for the development of the method. The term attributable was first introduced by Milani et al. (2001) in order to use it to find identifications of asteroids with known orbits, such identification is called attribution. This definition is very important at the present time in the development of methods for orbital determination, for example, methods developed by (Klokacheva 1991; Sokolskaya 1997) do not use this type of approach.

As mentioned earlier, due to a large number of observations available, it is necessary to resort to new methods that allow the monitoring of minor bodies of the Solar System. For this purpose, modern methods can be used based on the invariants of the motion (Villarraga & Quintero-Salazar 2016) or those that estimate the admissible region in which minor bodies can be found (Milani et al. 2004; Mosquera & Salazar 2019).

However, to give a more straightforward solution when there are VSA observations, it is possible to resort to variations of classical methods, such as Laplace (Milani et al. 2008). In this method, computational implementation is simpler compared to modern methods. Through this method, it is possible to know the orbital parameters of minor bodies. This latter is especially important to observatories that are initiating minor bodies observation programs.

Therefore, the purpose of this paper is to discuss a modification of the geocentric Laplace's method that allows for the use of all available astrometric data of the body in question. Moreover, we show how from this data it is possible to calculate the geodesic curvature value of the trajectory and how we use it to find the topocentric distance. Subsequently, we proceed to find the initial values for the problem of the two bodies, to determine the orbital parameters that describe the trajectory of the object.

Finally, the ephemerides for the subsequent re-observation of the minor bodies under study are estimated. Input data, OAUTP observations [MPC code W63 (Villarraga et al. 2017)], among others, available in the MPC database were used. To validate the results, we made a comparison with the NASA JPL Web-Interface service, both for orbital parameters (obtaining orbital errors for shape and orientation), and for the estimated ephemerides (when determining if the object is recoverable or not from a field of vision of 95’ × 72’).

This article is structured as follows: first, the methodology used in this study to recover minor bodies is described in detail. Next, the proposed methodology is validated by comparing the obtained results with those available in the NASA JPL Web-interface service. Finally, the main conclusions are presented, in which the reliability of the modification of the geocentric Laplace's method, by recovering minor bodies from the Solar System is demonstrated. This is useful even when the input data observations represent short time intervals.

2. MATERIALS AND METHODS

In this section, a modification of the geocentric Laplace's method to estimate a preliminary orbit is presented. The flowchart that illustrates this process is shown in Fig. 1.

Fig. 1. Flowchart of the developed algorithm for Initial Orbit Determination (IOD) using the modification of the geocentric Laplace's method.

Initially, a minor body \(\text{B}\), belonging to the Solar System that moves around the Sun with heliocentric position r is considered. This is observed from Earth ε which has a vector radius R, known for a determined moment of certain time (Fig. 2). The central idea is to find the state vector that describes this trajectory. To achieve this, a least squares fitting to the registered observations is made, and from this, an attributable is obtained. In the next stage, using the same fitting, the state vector of the Earth was estimated.

Fig. 2. Asteroid movement around the Sun and observed from Earth.

Furthermore, the methodology involves the calculation of the proper motion, the orthonormal base adapted to the apparent movement of \(\text{B}\) on the celestial sphere, and the constants of the geodesic curvature and along-track acceleration. Then, the dynamic and geometric equations were applied to find two additional constants that completed the set of constants necessary to estimate the magnitude of the heliocentric position of \(\text{B}\). Ultimately, we made a selection of the allowable solutions to find the state vector; with this, it was possible to estimate both the orbital elements and ephemerides. In the following section, methodology is explained in detail.

2.1 Computing attributable

One of the input parameters of the algorithm is the set of astrometric data obtained from different observation nights (Fig. 1, Block 1). Generally these data are expressed as follows:

\(\left(t_{i},\alpha _{i},\delta _{i}\right)i = 1,…,m\)       (1)

Where t is the observation time, α the right ascension, δ the declination, where m ≥ 3 is the total number of registered observations.

Then, a polynomial model fitting was performed for an average time t from the t_i (Milani & Gronchi 2010), of the form:

\(\xi \left(t\right) = \xi \left(\overset{¯}{t}\right) + \overset{˙}{\xi }\left(\overset{¯}{t}\right)\left(t−\overset{¯}{t}\right) + \frac{1}{2}\overset{¨}{\xi }\left(\overset{¯}{t}\right)\left(t−\overset{¯}{t}\right)^{2}\)       (2)

Where ξ represents the fitting for α and δ, respectively. From these quantities, the attributable vector \(A = \left(\alpha ,\overset{˙}{\alpha },\overset{¨}{\alpha },\delta ,\overset{˙}{\delta },\overset{¨}{\delta }\right)\) which is associated with a covariance matrix that yields an estimate of errors (Milani et al. 2001; Fig. 1, Block 2), was calculated.

At this point, the information contained in the attributable A does not include either the value of the topocentric distance ρ or the radial velocity \(\overset{˙}{\rho }\) of the body under study. Then, the next step consists in finding the value of such quantities that allow resolving the IOD, determining the orbital elements of a minor body from angular information only.

2.2 Poincaré interpolation

The next step is to estimate the state vector of the Earth (R, \(\overset{˙}{R}\)) , that is, the heliocentric position of the center of the Earth (Fig. 1, Block 3). In order to obtain more precise results, Poincaré (Poincaré 1906; Milani et al. 2008) suggested that these quantities should be calculated consistently with the same interpolation used to calculate the attributable vector and not using exact formulas. Therefore, the position of the Earth for t_i was estimated, and finally, (R, \(\overset{˙}{R}\)) was estimated using the eq. (2). In this case, a geocentric approach was implemented (Milani et al. 2007; Curtis 2014).

2.3 Computing proper motion

Both the position and heliocentric velocity of \(\text{B}\) are given by:

\(\mathrm{r}=\mathrm{R}+\rho \hat{\rho} \quad \dot{\mathrm{r}}=\dot{\mathrm{R}}+\dot{\rho} \hat{\rho}+\rho \dot{\alpha} \hat{\rho}_{\alpha}+\rho \dot{\delta} \hat{\rho}_{\delta}\)       (3)

Where \(\hat{\rho}_{\alpha}=\frac{\partial \hat{\rho}}{\partial \alpha} \text { and } \hat{\rho}_{\delta}=\frac{\partial \hat{\rho}}{\partial \delta}\).

Being:

\(\hat{\rho}=(\cos \alpha \cos \delta, \sin \alpha \cos \delta, \sin \delta)\)

\(\hat{\rho}_{\alpha}=(-\sin \alpha \cos \delta, \cos \alpha \cos \delta, 0)\)

\(\hat{\rho}_{\delta}=(-\cos \alpha \sin \delta,-\sin \alpha \sin \delta, \cos \delta)\)

In the previous set of equations, it is proven that:\(\hat{\rho} \cdot \hat{\rho}_{\alpha}=\hat{\rho} \cdot \hat{\rho}_{\delta}=\hat{\rho}_{\alpha} \cdot \hat{\rho}_{\delta}=0 \cdot\|\hat{\rho}\|=\left\|\hat{\rho}_{\delta}\right\|=1,\left\|\hat{\rho}_{\alpha}\right\|=\cos \delta\) . The next step of the method is to compute these vectors (Fig. 1, Block 4). Laplace's method is a technique used to approximate the values of ρ and \(\overset{˙}{\rho }\), so that eq. (3) can be determined. Being the topocentric position of \(\rho=\rho \hat{\rho} \text { of } \mathbf{B}\), topocentric velocity is giving by:

\(\frac{d \rho}{d t}=\dot{\rho} \hat{\rho}+\rho \frac{d \hat{\rho}}{d t}\)       (4)

Using an orthonormal base adapted to the apparent trajectory of \(\text{B}\) in the celestial sphere, that is, the image of \(\hat{\rho}(t)\) it follows that (Danby 1962; Milani & Gronchi 2010):

\(\boldsymbol{V}=\frac{d \hat{\rho}}{d t}=\frac{d \hat{\rho}}{d s} \frac{d s}{d t}=\eta \hat{V}\)       (5)

Where s is arc length parameter.

Then, proper motion \(\eta = \left\|\text{v}\right\|\) was established. Its value is one of the most important quantities that was estimated using the modification of the geocentric Laplace's method (Fig. 1, Block 5):

\(\eta=\left\|\frac{d \hat{\rho}}{d t}\right\|=\sqrt{\dot{\alpha}^{2} \cos ^{2} \delta+\dot{\delta}^{2}}\)       (6)

Likewise, the unit vector \(\hat{\mathrm{V}}\) was set, so that:

\(\eta=\frac{d s}{d t^{\prime}} \quad \hat{V}=\frac{d \hat{\rho}}{d s}\)       (7)

2.4 Computing vectors based on the attributable

Now, \(\hat{\mathrm{V}}\) has unit magnitude and \(\frac{d \hat{\rho}}{d s}\) is a perpendicular vector to \(\hat{\mathrm{V}}\), so that \(\hat{v} \cdot \hat{\rho}=0\). This allows to define a threedimensional basis if the following vector is added:

\(\hat{n}=\hat{\rho} \times \hat{V}\)       (8)

The mathematical notation of the derivative with respect to s is established as prime; and according to the above relationships, the next properties are satisfied:

\(\begin{gathered} \hat{V}^{\prime} \cdot \hat{\rho}=\frac{d}{d s}(\hat{V} \cdot \hat{\rho})-\hat{V} \cdot \hat{\rho}^{\prime}=-1 \\ \hat{V}^{\prime} \cdot \hat{V}=\frac{1}{2} \frac{d}{d s}\|\hat{V}\|^{2}=0 \end{gathered}\)       (9)

Now, the next derivatives are introduced:

\(\begin{aligned} &\alpha^{\prime}=\dot{\alpha} \frac{1}{\eta}, \quad \alpha^{\prime \prime}=\frac{1}{\eta^{3}}(\ddot{\alpha} \eta-\dot{\alpha} \dot{\eta}) \\ &\delta^{\prime}=\dot{\delta} \frac{1}{\eta}, \quad \delta^{\prime \prime}=\frac{1}{\eta^{3}}(\ddot{\delta} \eta-\dot{\delta} \dot{\eta}) \end{aligned}\)

From the second expression in (7), and considering the previous derivatives for the unit vectorv \(\hat{\mathrm{V}}\), it follows:

\(\hat{V}=\alpha^{\prime} \frac{\partial \hat{\rho}}{\partial \alpha}+\delta^{\prime} \frac{\partial \hat{\rho}}{\partial \delta}=\frac{1}{\eta}\left(\dot{\alpha} \hat{\rho}_{\alpha}+\dot{\delta} \hat{\rho}_{\delta}\right)\)       (10)

To unit vector \(\hat{n}\), the previous expression is replaced in (8), then, the following expression was obtained:

\(\hat{n}=\alpha^{\prime} \cos \delta \hat{\rho}_{\delta}-\frac{\delta^{\prime}}{\cos \delta} \hat{\rho}_{\alpha}=\frac{1}{\eta}\left(\dot{\alpha} \cos \delta \hat{\rho}_{\delta}-\frac{\dot{\delta}}{\cos \delta} \hat{\rho}_{\alpha}\right)\)       (11)

Since eqs. (10) and (11) are in terms of attributable values, the next step is to calculate these values (Fig. 1, Block 6).

2.5 Computing the geodesic curvature and along-track acceleration

According to the properties of (9) vector \(\hat{\mathrm{V}}\) can be expressed as:

\(\hat{V}^{\prime}=-\hat{\rho}+\kappa \hat{n}\)       (12)

The scalar quantity denoted by κ, is the geodesic curvature, which measures the deviation of the path of a great circle (a geodesic on the sphere). This value is unknown and hence is one of quantities to be resolved.

From the eq. (5) the derivative with respect to time is obtained, and then, replaced into (12), so that:

\(\frac{d^{2} \hat{\rho}}{d t^{2}}=\dot{\eta} \hat{V}+\eta^{2}(-\hat{\rho}+\kappa \hat{n})\)       (13)

In this equation, there is one additional scalar quantity \(\overset{˙}{\eta }\), which was estimated from the attributable as showed latter. To find this scalar quantity the following relationships are considered:

\(\begin{gathered} \hat{\rho}_{\alpha \alpha}=(-\cos \delta \cos \alpha,-\cos \delta \sin \alpha, 0) \\ \hat{\rho}_{a \delta}=(-\sin \delta \sin \alpha,-\sin \delta \cos \alpha, 0) \\ \hat{\rho}_{\delta \delta}=(-\cos \delta \cos \alpha,-\cos \delta \sin \alpha,-\sin \delta) \end{gathered}\)

Where \(\hat{\rho}_{a \alpha}=\partial\left(\hat{\rho}_{\alpha}\right) / \partial \alpha, \hat{\rho}_{\alpha \delta}=\partial\left(\hat{\rho}_{\delta}\right) / \partial \alpha, \hat{\rho}_{s s}=\partial\left(\hat{\rho}_{\alpha}\right) / \partial \delta\), and also the following properties are established:

\(\begin{gathered} \hat{\rho}_{a \alpha} \cdot \hat{\rho}_{a}=\hat{\rho}_{s \delta} \cdot \hat{\rho}_{a}=\hat{\rho}_{a \delta} \cdot \hat{\rho}_{s}=\hat{\rho}_{\delta \delta} \cdot \hat{\rho}_{\delta}=0 \\ \hat{\rho}_{\alpha \delta} \cdot \hat{\rho}_{a}=\sin \delta \cos \delta=\hat{\rho}_{a x} \cdot \hat{\rho}_{\delta}=\sin \delta \cos \delta \end{gathered}\)       (14)

In order to find κ, scalar product with \(\hat{n}\) is calculated to both sides of eq. (12):

\(\tilde{\kappa}=\hat{V}^{\prime} \cdot \hat{n}\)       (15)

Then eq. (10) is derived with respect to s. Simplifying it yields:

\(\hat{V}^{\prime}=\left(\alpha^{\prime \prime} \hat{\rho}_{\alpha}+\delta^{\prime \prime} \hat{\rho}_{\delta}\right)+\left(\alpha^{\prime 2} \hat{\rho}_{\alpha \alpha}+2 \alpha^{\prime} \delta^{\prime} \hat{\rho}_{\alpha \delta}+\delta^{\prime 2} \hat{\rho}_{\delta \delta}\right)\)

Finally, the previous expression is substituted into (15) and, making the necessary simplifications and substitutions involving (14), an expression for κ is obtained:

\(\kappa=\frac{1}{\eta^{3}}\left[\cos \delta(\dot{\alpha} \ddot{\delta}+\dot{\delta} \ddot{\alpha})+\dot{\alpha} \sin \delta\left(\dot{\delta}^{2}+\eta^{2}\right)\right]\)       (16)

Therefore, from eq. (16), it was possible to compute the geodesic curvature based on the variables of the attributable, which is, as previously mentioned, one of the quantities to be estimated in this method (Fig. 1, Block 7). Otherwise, to find the value of \(\overset{˙}{\eta }\), the scalar product on both sides of the eq. (13) with \(\hat{V}\) was performed. Thus we have:

\(\dot{\eta}=\frac{d^{2} \hat{\rho}}{d t^{2}} \cdot \hat{V}\)       (17)

Finally, using both eqs. (17) and (6), the value for alongtrack acceleration was obtained:

\(\dot{\eta}=\frac{1}{\eta}\left(\ddot{\alpha} \dot{\alpha} \cos ^{2} \delta-\dot{\alpha}^{2} \dot{\delta} \cos \delta \sin \delta+\ddot{\delta} \dot{\delta}\right)\)       (18)

Then, we proceeded to estimate this value numerically (Fig. 1, Block 8).

2.6 Computing constants through dynamical and geometric equations

By replacing (5) into eq. (4), the next expression is obtained:

\(\dot{\rho}=\dot{\rho} \hat{\rho}+\rho \eta \hat{v}\)

Deriving this expression again with respect to time and using (12), yields:

\(\ddot{\rho}=\left(\ddot{\rho}-\rho \eta^{2}\right) \hat{\rho}+(2 \dot{\rho} \eta+\rho \dot{\eta}) \hat{V}+\rho \eta^{2} \kappa \hat{n}\)       (19)

The two-body dynamic model was assumed in order to obtain preliminary orbits. In this model, the minor body is only affected gravitationally by the Sun; and therefore, no other external force alters its trajectory. Using Newton's Law of gravitation, the heliocentric accelerations of the Earth and the studied body are obtained, respectively, as:

\(\ddot{R}=-\frac{\mu}{R^{3}} R, \ddot{r}=-\frac{\mu}{r^{3}} r\)       (20)

From eq. (3) it is know that:

\(\ddot{\rho}=\ddot{r}-\ddot{R}=\mu\left(\frac{1}{R^{3}}-\frac{1}{r^{3}}\right) R-\frac{\mu}{r^{3}} \rho\)       (21)

A scalar product with \(\hat{V}\) in both eqs. (21) and (19) was performed. The result was equalized, and then the terms were organized as follows:

\(2 \dot{\rho} \eta+\rho \dot{\eta}=\mu\left(\frac{1}{R^{3}}-\frac{1}{r^{3}}\right)(R \cdot \hat{V})\)       (22)

This equation shows that if the values of ρ and r are known, the values for velocity \(\overset{˙}{\rho }\) can be found. Then, a scalar product with \(\hat{n}\); on both eq. (21) and (19) is performed, this yields:

\(\rho \eta^{2} \kappa=\mu\left(\frac{1}{R^{3}}-\frac{1}{r^{3}}\right)(R \cdot \hat{n})\)       (23)

This expression is organized, and the result is multiplied on both sides by R³, so it is reduced to:

\(C \frac{\rho}{R}=1-\frac{R^{3}}{r^{3}} \)       (24)

Where

\(C=\frac{\eta^{2} \kappa R^{3}}{\mu(\hat{R} \cdot \hat{n})}\)

Eq. (24) is known in literature as dynamical equation (Milani & Gronchi 2010).

2.7 Computing the magnitude of the heliocentric position

At this point we proceeded to perform the numerical estimation of the variables of interest \(r,\overset{˙}{r},\rho ,\overset{˙}{\rho }\). In principle, there is a geometric relationship between R, r, and ρ called the geometric equation (Danby 1962):

\(r^{2}=R^{2}+\rho^{2}+2 R \rho \cos \varepsilon\)       (25)

Where \(\cos \varepsilon=\hat{R} \cdot \hat{\rho}\). Up to this point, the constant C and cos ε are determined numerically (Fig. 1, Block 9) since the variables on which they depend are known quantities. From eq. (24), then:

\(\rho^{2}=\frac{R^{2}}{C^{2}}\left(1-\frac{2 R^{3}}{r^{3}}+\frac{R^{6}}{r^{6}}\right)\)       (26)

By replacing the previous expression into (25) and multiplying it by C²r⁶, the following expression is obtained:

\(\begin{aligned} P(r)=& C^{2} r^{8}-R^{2} r^{6}\left(C^{2}+2 C \cos \varepsilon+1\right) \\ &+2 R^{5} r^{3}(1+C \cos \varepsilon)-R^{8}=0 \end{aligned}\)       (27)

The above expression is a polynomial equation of the 8^th degree in r. By solving it, the value of r is found (Fig. 1; Block 10).

2.8 Determination of State Vector, Orbital Elements and Ephemerides

In order to find a preliminary orbit and subsequently propagate ephemerides, the following Matlab functions were implemented:

• Solutions for r that were either imaginary or negative were filtered, as well as the trivial solution r = R (Milani et al. 2008; Fig. 1, Block 11).

• Eq. (25) for ρ with each filtered value of r was solved (Fig. 1, Block 12). Then, the invalid solutions for ρ were filtered (Fig. 1, Block 13).

• By using the eq. (22), the corresponding value for radial velocity was obtained \(\overset{˙}{\rho }\) (Fig. 1, Block 14).

• Using the already-known magnitudes, the state vector (3) was determined (Fig. 1, Block 15).

• The orbit to generate ephemerides using the f and g series was propagated (Danby 1962; Fig. 1, Block 16).

• The correction for aberration was applied (Fig. 1, Blocks 16 to 18); then, the state vector was converted to orbital elements.

• Finally, the orbital elements and the ephemerides of the minor body were obtained (Fig. 1, Block 19).

2.9 Analysis of Errors

One way to carry out the error analysis is to compare all the estimated orbital parameters with the nominal ones; however, this is quite difficult (Schaeperkoetter 2011). An alternative is to reduce the number of parameters to be compared, considering that in any orbit the orbital parameters can be grouped into two groups, some associated with shape and others associated with orientation. The first contains information regarding the orbit shape on orbit plane: the semi major axis and eccentricity (a, e), while the second group contains information in respect to the relative orbit orientation: inclination, longitude of ascending node, and argument of perigee (i, Ω, ω). These error parameters are referred to as error in shape d, and orbit orientation Φ, respectively. This technique consists in reducing the number of parameters to be compared from 5 to 2 (Mortari et al. 2006). This technique has the advantage of making the graphical interpretation of the error easier, using a polar representation. In our case, the latter method is useful because we are using several sample asteroids, and with this analysis, it is possible to get a more global idea of the results of the method.

If the nominal parameters are (a, e, i, Ω, ω, θ) and the estimates are (a^* , e^* , i^*, Ω^* , ω^*, θ^*) , where θ is true anomaly, the error in the shape d is calculated as:

\(d=\sqrt{\left(a-a^{*}\right)^{2}+\left(b-b^{*}\right)^{2}}\)       (28)

Where \(b = a\sqrt{1−e^{2}}\) is the minor semi-axis. Orbit shape is related with the major and minor semi-axis which are dimensionally consistent parameters. Using both parameters the orbit shape can be represented as a point in the a − b plane. Therefore, the orbit error shape can be simply described by the eq. (28).

For orbit orientation, error Φ is calculated as:

\(\cos \Phi=\frac{1}{2}\left(t r\left[C C^{*^{T}}\right]-1\right)\)       (29)

Where

\(C=R_{3}(\omega+\theta) R_{1}(i) R_{3}(\Omega)=\left[\begin{array}{lll} \hat{r} & \hat{h} \times \hat{r} & \hat{h} \end{array}\right]^{T}\)       (30)

 

Being R₁ and R₃ the rotation matrices around the first and third coordinate axis, respectively. Mathematically, angle Φ is the main angle of corrective matrix between matrices C C and C^*. This can be interpreted as the angle between the directions of the estimated angular momentum and the real one; if this angle is small, the estimated orbit resembles the real one in terms of its position and orientation in three-dimensional space.

In the description of an orbit error, there are two parts: an error in the orientation Φ, which is an angle, and another in the position d, which is a distance. Hence, a polar representation is used to describe the error of the orbit through the point (Mortari et al. 2006):

\(E = de^{/ \Phi }\)       (31)

In this representation, the error will be smaller, the closer the point is to the origin of coordinates, and smaller the angle of inclination for the positive x-axis.

3. RESULTS

The developed algorithms were implemented as Matlab functions. The set of asteroids listed in Table 1 was used to validate the proposed methodology. As input data, observations reported to MPC (https://www.minorplanetcenter. net/db_search) and by OAUTP (MPC code W63) were used. Because of in the OAUTP we carry out observation programs in minor bodies, we decided to test our methodology with both NEA and MBA asteroids. From these asteroids, we had already made observation reports, and we complemented with those recorded by other observatories. These observations cover the observation intervals shown in the second column of Table 1 (in all cases, the observational interval used is a VSA). In the web site (https://www.utp. edu.co/observatorioastronomico/astrometria/orbital-elements- and-recovery-using-laplace), full data tables of the observations used in this study are available.

Table 1. Samples of the asteroid recovery. The number of input data N, interval of the initial observations Δt, the orbit type, the proper motion η, and the ρ value are shown.

The Poincaré interpolation and the calculation of the attributable vector were performed to apply the modification of the geocentric Laplace's method. This process resulted in a state vector that was used to determine orbital elements of each of the studied asteroids. The results are shown in Table 2.

Table 2. Orbital elements: Semimajor axis (a), eccentricity (e), inclination (i), the argument of perihelion (Ω), ascending node (ω), time of perihelion passage (T

Otherwise, Fig. 3 shows the orbits in the heliocentric ecliptic system, as well as its projections on the axes of the NEA asteroids and Hungarian-type asteroid. Fig. 4 shows the same for MBA-type asteroids. The animations of the trajectories obtained in this study are available at https:// www.utp.edu.co/observatorioastronomico/astrometria/ orbital-elements-and-recovery-using-laplace.

Fig. 3. Orbits of the NEA and Hungaria-type asteroids. Plots are in a heliocentric ecliptic system.

Fig. 4. Orbits of the MBA-type asteroids. Plots are in a heliocentric ecliptic system.

To analyze the obtained results, the errors in shape and orientation of the orbits for the set of 8 asteroids were estimated. The orbital elements data taken from the NASA JPL Web-interface service (https://ssd.jpl.nasa.gov/horizons. cgi#results) were assumed to be true. The results are shown graphically in Fig. 5.

Fig. 5. Estimated orbital errors for each asteroid: (a) The error in shape, and (b) the error in orientation.

As can be seen, the errors in shape are less than 53 × 10^–3 AU, except for the 2019 JA8 asteroid (Fig. 5a). These results are comparable with the real model of the orbits, which is intrinsically linked to the use of a two-body model (Gronchi et al. 2010). The error orientation Φ is less than 0.1 rad, except for the (4690) Strasbourg (Fig. 5b).

The errors in each of the estimated orbits in a polar representation are shown in Fig. 6. It is observed that the orbits are of good quality even when the dynamic model used only considers two bodies. Fig. 6a illustrates the object that shows the highest estimated orbital error is 2019 JA8. This may be because this asteroid was at an average distance from the observer of less than 0.35 AU, with a significant proper motion (Table 1). Under such conditions, Laplace's method should consider the position of the observer on Earth (which implies the reformulation of the entire mathematical model). Additionally, due to the proximity of the asteroid to the Earth, the gravitational force influence of the planet should be considered when calculating. That is, a three-body model should be implemented. On the other hand, in thcase of (4690) Strasbourg, the orientation error was greater than the other asteroids. This may be because this one was tested with a shorter time interval. However, as we clarified in section 2.9, seeing the error in a polar representation is important. Despite the aforementioned challenges, the estimated orbit for both 2019 JA8 and (4690) Strasbourg are of enough quality to recover this object again.

Fig. 6. Orbit error complex representation for the asteroids. (a) All used asteroids, and (b) Zooming into the plot.

In the rest of the cases, when zooming into the plot (Fig. 6b), the errors get close enough to the origin, and with minimal angles (the latter also for 2019 JA8). In this way, it is observed that the estimated orbits were of good quality.

However, despite the proximity between the estimated orbital parameters in respect to the real ones, the effective utility for the IOD method is determined by the capacity it has to produce ephemerides with enough precision so that the body can be re observed after initial observations. The purpose of our algorithm is to establish an IOD with the available information of the object. The IOD allows for the establishment of an ephemeris to perform the recovery of that object. This recovery allows for the gathering of a great number of observations. The more observations, the better the estimation of the orbit, which is the least-squares principle procedure input (Danby 1962). In our case, we tested the recovery simulation, calculating the ephemeris for 15, 30, 45, and 60 days after the last observation for each of the asteroids presented in Table 3. As it is evident in Table 3, it was possible to recover all the objects. The results are compared to those provided in the NASA JPL Web-interface service. During this time, it would be possible to collect a significant amount of additional observations that would improve the IOD.

Table 3. Ephemerides for sample asteroids. Date shows time in UT for the ephemerides at 15, 30, 45 and 60 days after the last observation entered as input. The symbols α and δ refer to the coordinates in right ascension and declination, respectively.

To simulate a possible recovery of each of the objects, the error in right ascension α and declination δ, were calculated as \(\left(\Delta \alpha ,\Delta \delta \right) = \left|\alpha _{v}−\alpha _{a},\delta _{v}−\delta _{a}\right|\).

Where (α_v, δ_v) are the NASA JPL coordinates and (α_a, δ_a) are obtained in this study by propagating the state vector obtained by using modification of the geocentric Laplace's method to the dates indicated in Table 3.

Columns 4 and 5 of Table 3 show the errors calculated for each of the coordinates. For example, the highest error is recorded at 41.94 arcmin in right ascension, for the 2019 JA8 asteroid, during a period of two months after its last observation used as input. As can be seen, the highest recorded error in decline was 2003 GW with 31.44 arcmin. At the mean time of the observations, both asteroids were at a very close distance to Earth (as can be seen in Table 1) and therefore their mean motion are greater than the others. This explains the reason why major errors occurred. In both cases, the obtained values make the re-observation very likely for these dates, while considering the estimation of the initial state vector was made from a geocentric approach.

The boundary conditions for the modification of the geocentric Laplace's method are determined by the convergence of the method. In this case, the convergence is subject to the fact that it is a VSA; therefore, it would not be possible to apply it to data that have too short observation intervals (TSA, too short arcs). In this case, it would be necessary resort to other methods developed by Mosquera & Salazar (2019). Geocentric Laplace's method would not be optimal for objects that get too close to Earth at the time of their observation (NEAs at their closest approach); moreover, it would not be suitable for studying objects that are at very large distances (that present a very low proper motion) as in the case of observing trans Neptunian objects.

Likewise, the Figs. 7 and 8 show the ephemerides obtained from the NASA JPL (purple dots) compared to those estimated in this study (pink squares). The field of vision given by the instrumental assembly of the OAUTP is 95´ × 72´ (blue dotted-line box).

Fig. 7. Calculated ephemerides for dates after the last recorded input for NEA and Hungaria-type asteroids: from row one, the results are shown for: (1864) Daedalus, 2003 GW (4690) Strasbourg and 2019 JA8, respectively. Blue dotted-line box is the field of vision given by the instrumental assembly of the OAUTP.

Fig. 8. Calculated ephemerides for dates after the last recorded input for MBA-type asteroids: from row one, the results are shown for: (1738) Oosterhoff (2717) Tellervo (1568) Aisleen y (2235) Vittore, respectively. Blue dotted-line box is the field of vision given by the instrumental assembly of the OAUTP.

According to Figs. 7 and 8, all the studied asteroids would be recoverable from the OAUTP within the field of observation, making the recovery procedure possible as well as a subsequent refinement of the initial calculation of the orbit. The most significant error between the calculated ephemerides and the actual position was for 2019 JA8 asteroid. This result was expected, given the error in the polar representation (Fig. (a)). However, despite this, the object can still be recovered within the visual field.

4. CONCLUSIONS

It was demonstrated that the modification of the geocentric Laplace's method, which considers an input given by a set of astrometric measurements of a minor body (generally more than three observations), permitted the estimation of values for quantities such as proper motion and geodesic curvature; with which the algorithm's equations were solved. From the data used, the one that covers the longest time interval is approximately 22 days. It was possible to verify that, despite the short time intervals, the methodology allows for the production of initial orbits of enough quality to subsequently re-observe the bodies.

The errors of the estimated orbit concerning those provided by the HORIZONS Web-Interface service of NASA's JPL were estimated. These errors were grouped in shape and orientation for each of the orbits so that the errors in shape were less than 53 × 10^–3 AU, except for the asteroid 2019 JA8, and in terms of orientation, the errors were less than 0.1 rad, except for (4690) Strasbourg. Both sets of errors are expected values when using the two-body problem as a dynamic model, i.e. orbital shape and orientation are expected values for Keplerian type orbital parameter.

Additionally, the ephemeris for each of the bodies studied, up to two months after the last observation used as input, was estimated. These results were compared with the NASA JPL Web-interface service. In all cases, it is shown that they are recoverable within the field of vision of the instrumentation of the OAUTP (95´ × 72´). The highest error obtained in α was 41.94´ and 31.44´ for δ, respectively for (2019 JA8) and (2003 GW) asteroids. Therefore, the proposed methodology shows the importance of modification of the geocentric Laplace's method . The method is relatively simple to implement, and the results allow for the recovery of minor bodies from the Solar System, even when the input data observations represent very short arcs.

This study represents the second step in the minor bodies observation program of the Solar System implemented in the Astronomical Observatory of the Technological University of Pereira, after obtaining the MPC code (W63). The codes written in Matlab (https://www.utp.edu.co/observatorioastronomico/ astrometria/orbital-elements-and-recovery-using-laplace) can be useful for observatories that begin their process of studying the IOD of asteroids.