A&A 441, 815-818 (2005)
DOI: 10.1051/0004-6361:20053160

Orbit correction without variational equations

The orbits of Caliban, 46P/Wirtanen and 67P/Churyumov-Gerasimenko

K. Aksnes ¹ - T. Grav²

1 - Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern 0315 Oslo, Norway
2 - Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

Received 30 March 2005 / Accepted 20 June 2005

Abstract
Differential orbit correction is performed by computing an osculating orbit from observations from which the perturbations have been subtracted. No variational equations need to be integrated except to solve for non-gravitational forces. The perturbations in the observations are obtained by subtracting a two-body orbit from an integrated, perturbed orbit. The method is used to update the orbits of the Uranian satellite Caliban and of the comets 46P/Wirtanen and 67P/Churyumov-Gerasimenko with recent observations. The method is shown to agree very well with, but is considerably faster than, the conventional method based entirely on variational equations. Algorithms are given for both methods.

Key words: celestial mechanics - comets: general - planets and satellites: general

1 Introduction

This paper deals with orbit correction of solar system objects (comets, asteroids, satellites). Usually orbit correction is done by matching a perturbed orbit directly to the available observations, starting from an osculating orbit or state vector whose six elements are adjusted by least squares. This requires numerical integration of one orbit and six variational equations to derive the partial derivatives needed in the least squares solution. This method may be referred to as the coordinate perturbation method.

We present here a simpler and faster method which we shall call the observation perturbation method. The perturbations are first removed from the observations to obtain a set of unperturbed fictitious observations from which the osculating orbit is derived as a pure two-body orbit using analytic partial derivatives. The perturbations in the observations are obtained by subtraction between a numerically integrated, perturbed orbit and an analytic, unperturbed orbit in an iterative procedure. No variational equations need to be integrated except to solve for non-gravitational forces. The method was introduced by the first author who used it to compute the orbits of Jupiters's satellite Leda (Aksnes 1978) and Uranus' satellites Caliban and Sycorax (Marsden et al. 2000). The second author (Grav 2001) then applied the method to the comet 46P/Wirtanen and extended it by the inclusion of non-gravitational forces.

In the following we will give an exposition of the observation perturbation method and demonstrate its use through some examples. We will also compare the results with the corresponding ones obtained with the conventional coordinate perturbation method for which an algorithm is also given.

2 Orbit model

The equations of motion for a body of mass m attracted by a central mass m₀ and N perturbing masses m₁, m₂, ....m_N take the form

$$\ddot{\vec{r}} = -(\mu_0 + \mu)\frac{\vec{r}}{r^3} + \sum_{i... ...c{r}_i-\vec{r}}{\Delta_i^3} - \frac{\vec{r}_i}{{r_i}^3}\right)$$ (1)

where $\mu$, $\mu_0$ and $\mu_i$ are the gravitational constant times respectively the masses m, m₀ and m_i, $\vec{r}$ being the radius vector of the object whose distance from body i is $\Delta_i$. Let $\vec{s}_0=\{x_0,y_0,z_0,\dot{x_0},\dot{y_0},\dot{z_0}\}$ be the initial state vector (column vector) which we seek to adjust by fitting an orbit to observations by least squares, the O-C residuals being represented as a variation $\Delta\vec{r} = \{\Delta x,\Delta y,\Delta z\}$,

$$\Delta\vec{r} = G\Delta\vec{s}_0$$ (2)

where G is the 3 by 6 matrix

$$G = \left[ \begin{array}{cccc} \frac{\partial x}{\partial x_... ... & \frac{\partial z}{\partial \dot{z_0}} \end{array} \right] .$$ (3)

The elements of this matrix; i.e. the partial derivatives needed in the least squares adjustment, can be obtained by solving the differential equations (variational equations)

$$\ddot{G} = \frac{\partial\ddot{\vec{r}}}{\partial\vec{r}} G$$ (4)

where

$$\frac{\partial\ddot{\vec{r}}}{\partial\vec{r}}= (\mu_0+\mu)D(\vec{r}) + \sum_{i=1}^{N}\mu_iD(\vec{r}_i-\vec{r}),$$ (5)

D being the operator (Aksnes 1973)

$$D(\vec{r})= -\nabla(\frac{\vec{r}}{r^3})= \frac{3}{r^5}\left[... ...^2-r^2/3 & yz\\ zx & zy & z^2-r^2/3 \end{array} \right] \cdot$$ (6)

Initially all the elements of G and $\dot{G}$ are zero except for $G_{11}=G_{22}=G_{33}=\dot{G}_{14}=\dot{G}_{25}=\dot{G}_{36}=1$. If the observables are right ascension $\alpha$ and declination $\delta$, Eq. (2) is used to obtain

$$\left[\begin{array}{c}\Delta\alpha\cos{\delta}\\ \Delta\delta... ...& \cos{\delta} \end{array}\right]\frac{\Delta\vec{r}}{r} \cdot$$ (7)

Let $\alpha,\delta$ be the perturbed and $\alpha_0,\delta_0$ the unperturbed, computed right ascension and declination corresponding to the observed values $\alpha_{\rm obs},\delta_{\rm obs}$. By subtracting the perturbations from the observations, we get fictitious observations

$$\alpha_{0,\rm obs} = \alpha_{\rm obs}-(\alpha-\alpha_0)= \alpha_0+(\alpha_{\rm obs}-\alpha)= \alpha_0+\Delta\alpha$$

$$\delta_{0,\rm obs} = \delta_{\rm obs}-(\delta-\delta_0)= \delta_0+(\delta_{\rm obs}-\delta)= \delta_0+\Delta\delta .$$ (8)

Through the method of least squares the osculating orbit at the initial epoch is adjusted to minimize the sum of squares of the residuals $\Delta\alpha\cos{\delta_0}$ and $\Delta\delta$. This leads to the following equations of condition for the observation perturbation method,

$$\sum_{i=1}^6\cos{\delta_0}\frac{\partial\alpha_0}{\partial c_i} \Delta c_i = \Delta\alpha\cos{\delta_0}$$

$$\sum_{i=1}^6\frac{\partial\delta_0}{\partial c_i} \Delta c_i = \Delta\delta$$ (9)

where the $\Delta c$'s represent corrections to the initial state vector $\vec{s}_0$ or equivalently to orbital elements $a,e,M_0,I,\omega,\Omega$. The partial derivatives $\partial\alpha_0/\partial c_i$ and $\partial\delta_0/\partial c_i$ w.r.t. the orbital elements are derived by means of Eq. (7) and the two-body partials $\partial\vec{r}/ \partial c_i$ supplied in the Appendix. The integration and least squares solution are repeated with an improved initial orbit until convergence.

The above equations of condition also apply for the coordinate perturbation method if the unperturbed $\alpha_0,\delta_0$ are replaced by the perturbed $\alpha,\delta$. This then requires the partials in the G matrix obtained through numerical integration.

2.1 Non-gravitational perturbations

When dealing with comets, one frequently has to include forces due to outgassing from the nucleus. We here adopt a model due to Marsden et al. (1973). To Eq. (1) is added an acceleration decomposed in a radial, transverse and normal direction:

$$\Delta\ddot{\vec{r}} = g(r)\left(A_1\frac{\vec{r}}{r} + A_2\... ...\dot{\vec{r}}-\dot{r}\vec{r}}{h} + A_3\frac{\vec{h}}{h}\right)$$ (10)

where $\vec{h} = \vec{r}\times\dot{\vec{r}}$ and

$$g(r) = \alpha\left(\frac{r}{r_0}\right)^{-m}\left(1 + \left(\frac{r}{r_0}\right)^n\right)^{-k} \cdot$$ (11)

Beyond the distance r₀ the non-gravitational forces drop off very rapidly. Through curves describing the rate of vaporization of water ice from a rapidly rotating nucleus with a Bond albodo 0.1, the values r₀=2.808 AU, m=2.15, n=5.093 and k=4.6142 were deduced. The constant $\alpha=0.111262$ normalizes g(r) to unity at 1 AU, the unit for the acceleration being AU/day².

The parameters A₁, A₂ and A₃ are assumed to be solved for along with the state vector or orbital elements of the comet at some epoch. We then need the partial derivatives of the observables w.r.t. these parameters. We define a column vector

$$\vec{p}_i = \frac{\partial\vec{r}}{\partial A_i} = \left[\be... ...al A_i \\ \partial z/\partial A_i \end{array}\right], i=1,2,3$$ (12)

and derive the variational equations

$$\ddot{\vec{p}_i} = \frac{\partial\ddot{\vec{r}}}{\partial\vec... ...i + \frac{\partial\ddot{\vec{r}}}{\partial A_i}, i = 1,2,3 .$$ (13)

The partials $\partial\ddot{\vec{r}}/\partial A_i$ can be immediately deduced from Eq. (10), while the partials $\partial\ddot{\vec{r}}/\partial \vec{r}$ may be found from Eq. (5). The partials p_i are all initially zero.

2.2 Implementation

The preceding formulas have been implemented in two FORTRAN programs which we shall refer to as Model 1, the coordinate perturbation method, and Model 2, the observation perturbation method. Both programs use the same input:

1.

Epoch in Julian days.

2.

State vector ( $x,y,z,\dot{x},\dot{y},\dot{z}$) or orbital elements ( $T,\omega,\Omega,I,q,e$) at the epoch.

3.

ID's of central body and perturbing bodies ( ${\rm 1=Mercury}$, ${\rm 2=Venus}$, .... ${\rm 9=Pluto}$, ${\rm 10=Moon}$, ${\rm 11=Sun}$).

4.

Initial values (usually zero) of the non-gravitational parameters A₁,A₂,A₃ for comets.

5.

Initial and min step size in the integration and max number of iterations (normally stopped at a specified rms level).

6.

Observation file in the format used at the Minor Planet Center.

The output is the improved state vector and orbital elements and,if applicable, non-gravitational parameters. The observation residuals and their overall rms value are also given.

The integration is done with the Bulirsch-Stoer (1980) first-order method and adaptive step size control. This is a very efficient scheme of numerical integration. No additional interpolations are required to match the output to the times of observation.

In both models the state vector is integrated in accordance with Eq. (1) and, for comets, Eq. (10), and the 9 components of the non-gravitational partial derivatives p_i are obtained by integrating Eqs. (13) with the aid of Eq. (5). Model 1 requires the additional integrations of the variational Eqs. (4) for the 18 elements of G. This sums up to the integration of 30 s-order differential equations in Model 1 but only 12 in Model 2. Without non-gravitational perturbations, these numbers are reduced to 21 and 3 respectively. Model 2 thus offers a very substantial savings of integration time, especially if non-gravitational forces are excluded, although to some extent this savings will be offset by the need to compute two-body partials. Model 2 also appears to converge somewhat faster and be more stable, but this may be due to solving for corrections to the orbital elements instead of to the state vector, as in Model 1.

3 Applications to satellites and comets

We will illustrate the use of both models for one satellite and two comets.

3.1 Caliban (Uranus XVI)

Our first example is orbit improvement for the Uranian satellite Caliban based on 47 observations between June 1984 and November 1998, covering a little over 9 revolutions of the satellite which has a period of 579 days. Perturbations were included due to the Sun, Jupiter and Saturn. We started from an orbit with an average rms residual of 104 arcsec. Note that this is a rather large relative error because the orbit spans only about 17 arcmin in the sky. Model 1 diverged, while Model 2 converged with an rms residual of 0.63 arcsec after 4 integrations. With a more accurate starting orbit (rms = 37 arcsec), Model 1 converged on very nearly the same orbit, but it took 6 integrations. Thus the observation perturbation method appears to be more robust and stable than the coordinate perturbation method.

Table 1 lists the orbital elements of the satellite according to both solutions. The differences between the two sets of elements are well within the standard errors of the fits.

Table 1: Orbit of Caliban from 47 obs. 1984-1998.

3.2 46P/Wirtanen

Our second example involves the comet 46P/Wirtanen which was the original target of ESA's Rosetta mission. The orbit of this comet has been studied by several authors. Grav (2001) used the observation perturbation method to analyze observations between 1948 and 1999. He tried several models for the non-gravitational forces: (1) Marsden's model with and without a time shift $\tau$ making g(r) in Eq. (11) peak after perihelion if $\tau > 0$, (2) Sitarski's rotating nucleus model (Sitarski 1990), and (3) Sekanina's forced precession model as implemented by Krolikowska & Szutowicz (1999).

This comet is affected rather strongly by non-gravitational forces which have undergone substantial changes during the 50-year interval considered, probably mainly due to close approaches to Jupiter in 1971 and 1984. Just prior to the 1971 approach the comet's period was 6.7 years which was reduced to about 5.5 years after the 1984 approach. The same set of non-gravitational parameters gives satisfactory residuals for at most three apparitions of the comet.

Table 2 gives some results obtained by several authors for the non- gravitational parameters A₁, A₂ and A₃ for different observation intervals.

Table 2: Non-gravitational parameters in units of 10^-8 AU/day² for 46P/Wirtanen.

The best result (Grav 2001) was obtained with the use of the forced precession model with a rotating nucleus, as introduced by Krolikowska & Szutowicz (1999) for 46P/Wirtanen; see Table 3. It was necessarry to use two different values for the main parameter A: $A_{\rm I}$ before 1971 and $A_{\rm II}$ after 1971. The inclusion of a time shift $\tau$ in the second run did not have a noticeable effect on the residuals. The non- gravitational parameters are in rough agreement with those obtained by Krolikowska and Szutowicz, except that they found a much larger negative value for $\tau$ (-23.476). However, they used observations only up to 1995 and a different epoch (Jan. 27 1998), so the two sets of results are not directly comparable.

In all of Grav's calculations, the available observations were screened using the the Bielicki criterion for rejecting and weighting observations, as outlined by Bielicki & Sitarski (1991). Note that in the examples that follow (Tables 4 and 5), the observations are all given equal weight and are rejected if they produce residuals larger than 4 arcsec. The rms values of the fits can therefore be expected to be somewhat larger than if the Bielicki selection criterion had been used.

The orbital elements in Table 3 are in good agreement with the orbital elements in Table 4 derived from 555 observations between July 8 1991 and April 10 2003, using both model 1 and model 2 with the three Marsden parameters. Perturbations were included from all the nine planets and the Moon. The very good internal agreement of the two models again shows that the observation perturbation method works equally well as the coordinate perturbation method.

Table 3: Forced precession model for 46P/Wirtanen 1948-1999.

Table 4: Results for 46P/Wirtanen from 555 obs. 1991-2003.

3.3 67P/Churyumov-Gerasimenko

Our last example is orbit computation for the comet 67P/Churyumov- Gerasimenko which has a period of 6.6 yr. This comet is of special interest since the Rosetta spacecraft is scheduled to rendezvous with the comet in 2014. Table 5 gives the orbital elements and the non-gravitational parameters obtained for this comet by means of model 1 and model 2 based on 1301 observations between Sep. 9, 1969 and June 20, 2004. Again all nine planets and the Moon were included as perturbers (inclusion of the Moon reduced the rms residual by 0.2 arcsec). For comparison, we also show the corresponding results obtained by Krolikowska (2003) based on 1169 observations between Sep. 9, 1969 and March 13, 2003. The fact that the 1.16 arcsec rms residual obtained by Krolikowska is smaller than the 1.57 arcsec in model 1 and model 2, is probably due to Krolikowska's use of the Bielicki observation selection criterion. Even so, the three sets of orbital elements and non-gravitational parameters are in excellent agreement.

The small values for A₁, A₂ and A₃ show that this comet is much less affected by non-gravitational forces than is 46P/Wirtanen. For a discussion of the orbital evolution of this comet, the reader is referred to the article by Krolikowska (2003).

Table 5: Results for 67P/Churyumov-Gerasimenko.

4 Conclusions

Our results show that the observation perturbation method is fully compatible with the conventional coordinate perturbation method. The former method is simpler and requires less computer time; it also seems to have better convergence when applied to satellites. One might suspect that this method will handle very strongly perturbed orbits less well than the conventional method, but this is not borne out in the case of 46P/Wirtanen which had two close approahes to Jupiter, although both methods had problems satisfying the observations near those times.

Acknowledgements

We wish to thank Gareth V. Williams, Associate Director of the Minor Planet Center, Cambridge, Mass., for providing the observations on which this paper is based.

References

• Aksnes, K. 1973, Cel. Mech., 8, 99 [NASA ADS] (In the text)
• Aksnes, K. 1978, AJ, 83, 1249 [NASA ADS] [CrossRef] (In the text)
• Bielicki, M., & Sitarski, G. 1991, Acta Astron., 41, 309 [NASA ADS] (In the text)
• Grav, T. 2001, The Orbital Motion of Periodic Comet 46P/Wirtanen, Cand. Scient. (MS) thesis, University of Oslo (In the text)
• Krolikowska, M., & Sitarski, G. 1996, A&A, 310, 992 [NASA ADS]
• Krolikowska, M., & Szutowicz, S. L. 1999, A&A, 343, 997 [NASA ADS] (In the text)
• Krolikowska, M. 2003, Acta Astron., 53, 195 [NASA ADS] (In the text)
• Marsden, B. G., Sekanina, Z., & Yeomans, D. K. 1973, AJ, 78, 211 [NASA ADS] [CrossRef] (In the text)
• Marsden, B. G. 1983, Minor Planet Circular, 8274
• Marsden, B. G., Williams, G. V., & Aksnes, K. 2000, in Minor Bodies in the Outer Solar System (Berlin: Springer Verlag), 187 (In the text)
• Sitarski, G. 1990, Acta Astron., 44, 91 [NASA ADS] (In the text)
• Stoer. J., & Bulirsch, R. 1980, in Introduction to Numerical Analysis (New-York: Springer Verlag), Chap. 7 (In the text)

 

  
Online Material

Appendix: Two-body partials

Let $a,e,M_0,I,\omega,\Omega$ be the orbital elements, referred to the ecliptic, of an object in an elliptic orbit. The object's Earth-equatorial coordinates x,y,z centered on the primary may then be expressed as

$$\left[\begin{array}{c} x\\ y\\ z\\ \end{array}\right] = \left... ... \left[\begin{array}{c} \bar{x}\\ \bar{y}\\ \end{array}\right]$$ (14)

where $\bar{x}$ and $\bar{y}$ are components of the radius vector rin the orbital system,

$$\bar{x}= a(\cos{E}-e),\; \bar{y}=a\sqrt{1-e^2}\sin{E},\; r=a(1-e\cos{E})$$ (15)

Here E is given from Kepler's equation,

$$E-e\sin{E}=M=M_0+n(t-t_0),\; n=\sqrt{\frac{\mu}{a^3}},$$ (16)

t₀ being the epoch of the orbital elements. The elements of the rotation matrix are

$$\cos{\Omega}\cos{\omega} - \sin{\Omega}\cos{I}\sin{\omega}$$                                               R₁₁ =

$$-\cos{\Omega}\sin{\omega} - \sin{\Omega}\cos{I}\cos{\omega}$$ R₁₂ =

$$(\sin{\Omega}\cos{\omega}+\cos{\Omega}\cos{I}\sin{\omega}) \cos{\epsilon} - \sin{\omega}\sin{I}\sin{\epsilon}$$ R₂₁ =

$$(-\sin{\Omega}\sin{\omega}+\cos{\Omega}\cos{I}\cos{\omega}) \cos{\epsilon} - \cos{\omega}\sin{I}\sin{\epsilon}$$ R₂₂ =

$$\sin{I}\sin{\omega}\cos{\epsilon}+(\cos{\omega}\sin{\Omega} +\sin{\omega}\cos{\Omega}\cos{I})\sin{\epsilon}$$ R₃₁ =

$$\sin{I}\cos{\omega}\cos{\epsilon}+(-\sin{\omega}\sin{\Omega} +\cos{\omega}\cos{\Omega}\cos{I})\sin{\epsilon}$$ R₃₂ = (17)

where $\epsilon=23^\circ43928$ (J2000) is the obliquity of the ecliptic. We need the partial derivatives of x,y,z w.r.t. $a,e,M_0,I,\omega,\Omega$:

$$\frac{\partial x}{\partial a} R_{11}\frac{\partial\bar{x}}{\partial a} + R_{12}\frac{\partial\bar{y}}{\partial a}$$ =

$$\vdots$$ (18)

$$\frac{\partial z}{\partial \Omega} \frac{\partial R_{31}} {\partial \Omega}\bar{x} + \frac{\partial R_{32}}{\partial \Omega}\bar{y}$$ =

The partials of the R's w.r.t. $I,\omega,\Omega$ are readily derived, while the required partials of $\bar{x}$ and $\bar{y}$ are found to be

$$\frac{\partial\bar{x}}{\partial a} \frac{\bar{x}}{a}+ \frac{3a}{2r}n(t-t_0)\sin{E}$$ =

$$\frac{\partial\bar{x}}{\partial e} -a - \frac{\bar{y}^2}{r(1-e^2)}$$ =

$$\frac{\partial\bar{x}}{\partial M_0} -\frac{a\bar{y}}{r\sqrt{1-e^2}}$$ = (19)

$$\frac{\partial\bar{y}}{\partial a} \frac{\bar{y}}{a} - \frac{3a}{2r}n(t-t_0)\sqrt{1-e^2}\cos{E}$$ =

$$\frac{\partial\bar{y}}{\partial e} \frac{\bar{x}\bar{y}}{r(1-e^2)}$$ =

$$\frac{\partial\bar{y}}{\partial M_0} \frac{a}{r}\sqrt{1-e^2}(\bar{x}+ae)$$ =