A&A 397, 353-359 (2003)
DOI: 10.1051/0004-6361:20021518

An upgraded theory for Helene, Telesto, and Calypso

P. Oberti¹ - A. Vienne²

1 - Observatoire de la Côte d'Azur, BP 4229, 06304 Nice Cedex 4, France
2 - Institut de mécanique céleste et de calcul des éphémérides, UMR 8028 du CNRS and Université de Lille, 1 impasse de l'Observatoire, 59000 Lille, France

Received 1 August 2002 / Accepted 10 October 2002

Abstract
A tridimensional model including the perturbation due to Saturn's oblateness provides accurate solutions for the motions of Tethys and Dione's Lagrangian satellites Telesto, Calypso, and Helene. Expanded around the Lagrangian points and compared to numerical simulations to check their reliabilities, the solutions are fitted to almost twenty years of data thanks to recently published observations. User-friendly compact series give Telesto, Calypso, and Helene's positions and velocities in the mean ecliptic and equinox of J2000.

Key words: planets and satellites: individual: Saturn - planets and satellites: general - celestial mechanics

1 Introduction

Recently published observations of Tethys and Dione's Lagrangian satellites Telesto, Calypso, and Helene (Veiga & Vieira Martins 2000; Veiga et al. 2002), give an opportunity to revisit and improve the theory of their motions to face an almost 2-decade observation time span.

The mass ratios Tethys/Saturn ( $m/M=1.20\times10^{-6}$) and Dione/Saturn ( $m/M=1.85\times10^{-6}$) ensure the stability of the Lagrangian points L₄and L₅ in the Restricted Circular 3-Body Problem. Due to the 1:1-resonance with the companion, a satellite around these positions shows librations in its motion resulting in a tadpole-shaped orbit in a rotating reference frame. Expansions around the equilibrium positions provided accurate solutions in the bidimensional case (Deprit & Rom 1970).

Pointed out as a significant perturbation by numerical simulations, Saturn's oblateness (through the first zonal coefficient of Saturn's potential developed in spherical harmonics) was then introduced in the model by addition of perturbing terms in the second and following orders of the previously expanded non-perturbed Hamiltonian (Oberti 1990).

This new theory, now developed in a tridimensional context, deals with the perturbation at the uppermost level of the expansion, the quadratic part of the Hamiltonian (generally called H₀), improving the convergence of the solution and offering an explicit generalization suitable for other studies.

Moreover, the model could prove useful for the analysis of the CASSINI space probe's accurate observations of these faint satellites that should be obtained during its nearing exploration of the Saturnian system.

2 Analytical model

For the sake of convenience, computations are carried out in dimensionless units (Szebehely 1967). The distance between the primaries, and the sum of their masses (called $\mu=m/(M+m)$ for the second primary and $1-\mu$ for Saturn), are assumed to be unity. The reference frame is a barycentric rotating frame with a rotational period of $2\pi/k$ (in dimensionless time $\tau$), where the value of k will be computed in order to insure a consistent model involving fixed primaries. The xy-plane contains the primaries, the x-axis being located along the primaries and positively oriented from Saturn toward the second primary.

The theory will be expanded around L₄. Dealing with the symmetrical point L₅ requires only slight modifications indicated in Sect. 7.

With $\gamma=1-2\mu$, L₄ is located at the point $\left({\gamma\over2},{\sqrt{3}\over2},0\right)$.

Let x, y, z be the satellite's Cartesian coordinates and P_x, P_y, P_z their conjugate variables in the barycentric rotating frame. The Hamiltonian of the motion is:

$$\begin{eqnarray*}H={1\over2}\left(P_x^2+P_y^2+P_z^2\right)-k\left(xP_y-yP_x\righ... ...er r^2}\left(1-{3z^2\over r^2}\right)\right] -{\mu\over\Delta}, \end{eqnarray*}$$

where:

$$r=\sqrt{\left(x+\mu\right)^2+y^2+z^2},$$

$$\Delta=\sqrt{\left(x+\mu-1\right)^2+y^2+z^2},$$

and $J={1\over2}J_2r_{\rm e}^2$ stands for half the first Saturnian zonal harmonic coefficient J₂ (Campbell & Anderson 1989) multiplied by the square of Saturn's equatorial radius ( $r_{\rm e}=60~268$ km), appropriately scaled. The value of J, depending then on the selected second primary, is $3.4088\times10^{-4}$ for Tethys and $2.0780\times10^{-4}$for Dione.

L₄ is no longer an equilibrium position (Sharma & Subba Rao 1976). The conditions:

$$x={\gamma\over2},\quad y={\sqrt{3}\over2},\quad z=0,\quad \dot x=0,\quad\dot y=0,\quad\dot z=0,$$

lead to:

$$\begin{eqnarray*}\lefteqn{\dot P_x=\left(k^2-1\right){\gamma\over2} -{3J\over2}... ...ver2}\left({1+\gamma\over2}\right),} \\ \lefteqn{\dot P_z=0,} \end{eqnarray*}$$

where dots stand for derivation with respect to the dimensionless time. Canceling the accelerations would require two different values of k:

$$k_1=\sqrt{1+3J\left({1+\gamma\over2}\right)},\quad k_2=\sqrt{1+3J\left({1+\gamma\over2\gamma}\right)}\cdot$$

The value of k is in fact computed in order to immobilize the second primary, also perturbed by Saturn's oblateness. The point $\left(1,0\right)$ (or any point on the unit circle) turns into an equilibrium position for the Hamiltonian:

$$K={1\over2}\left(P_x^2+P_y^2\right)-k\left(xP_y-yP_x\right) -{1\over r}\left(1+{J\over r^2}\right),$$

where:

$$r=\sqrt{x^2+y^2},$$

when:

$$k=\sqrt{1+3J}.$$

With k₂>k>k₁, L₄ becomes an almost-equilibrium position:

$$\dot P_x=-{3J\over2}\mu,\quad \dot P_y={3J\sqrt{3}\over2}\mu,\quad\dot P_z=0.$$

A first canonical transformation defines a coordinate system in the vicinity of L₄:

$$\begin{eqnarray*}\lefteqn{x={\gamma\over2}+C\xi-S\eta,\quad P_x=-k{\sqrt{3}\ove... ...er2}+S P_\xi+C P_\eta,}\\ \lefteqn{z=\zeta,\quad P_z=P_\zeta,} \end{eqnarray*}$$

where S and C are the sine and cosine of a rotation devoted to removing a $\xi\eta$-term in the expression of the zero-order Hamiltonian below. This kind of coordinate system is generally referred to as the normal frame at an equilibrium position.

Applying the transformation turns the Hamiltonian into:

$$\begin{eqnarray*}H={1\over2}\left(P_\xi^2+P_\eta^2+P_\zeta^2\right) -k\left(\xi... ...1 -{3~\zeta^2\over r^2}\right)\right]-{{1-\gamma}\over2\Delta}, \end{eqnarray*}$$

where:

$$r=\sqrt{\left(\xi+S_1\right)^2+\left(\eta+C_1\right)^2+\zeta^2},$$

$$\Delta=\sqrt{\left(\xi+S_2\right)^2+\left(\eta+C_2\right)^2+\zeta^2},$$

$$S_1={S\sqrt{3}+C\over2},\quad C_1={C\sqrt{3}-S\over2},$$

$$S_2={S\sqrt{3}-C\over2},\quad C_2={C\sqrt{3}+S\over2}\cdot$$

3 Local solutions

Using, for example, the software Mathematica, we expand the functions of the coordinates in the vicinity of $\xi=\eta=\zeta=0$:

$$\begin{eqnarray*}{1\over r}&=&1-S_1\xi-C_1\eta+{1\over2}\left(3S_1^2-1\right)\xi... ...1\over2}\left(3~C_1^2-1\right)\eta^2 -{1\over2}~\zeta^2+\cdots, \end{eqnarray*}$$

$$\begin{eqnarray*}{1\over r^3}&=&1-3S_1\xi-3C_1\eta+{3\over2}\left(5S_1^2-1\right... ...3\over2}\left(5~C_1^2-1\right)\eta^2 -{3\over2}~\zeta^2+\cdots, \end{eqnarray*}$$

$$\begin{eqnarray*}{1\over r^5}&=&1-5S_1\xi-5C_1\eta+{5\over2}\left(7S_1^2-1\right... ...5\over2}\left(7~C_1^2-1\right)\eta^2- {5\over2}~\zeta^2+\cdots, \end{eqnarray*}$$

and similar expressions for $\Delta$ with the second index. Then, scaling with respect to the total power p+q+r of the monomials  $\xi^p\eta^q\zeta^r$, the Hamiltonian can be expanded up to an arbitrary order n.

Neglecting the remainder:

$$H_{\left(n\right)}=H_0+H_1+{1\over2}H_2+\cdots+{1\over n!}H_n,$$

H₀ containing the terms where $p+q+r\le2$.

Defining:

$$\alpha={3\over4}\left[1+5J\left({1+\gamma\over2}\right)\right... ...a={3\over4}\left[\gamma+5J\left({1+\gamma\over2}\right)\right],$$

the substitutions into H₀ lead to:

$$\begin{eqnarray*}H_0&=&{1\over2}\left(P_\xi^2+P_\eta^2+P_\zeta^2\right) -k\left... ...ft[-1-9J\left({1+\gamma\over2}\right)\right]{\zeta^2\over2}\cdot \end{eqnarray*}$$

The terms in $\xi$ and $\eta$, resulting from the fact that L₄ is not exactly an equilibrium position, are small because $J\left(1-\gamma\right)/2=J\mu$. At the precision required by the observations, they can safely be ignored, transforming this way the Lagrangian point L₄ in an actual equilibrium position for the theory.

Removing the $\xi\eta$-term leads to:

$$S=-\sqrt{\delta-1\over2~\delta},\quad C=\sqrt{\delta+1\over2... ...lta},\quad \delta=\sqrt{1+3\left({\beta\over\alpha}\right)^2},$$

and H₀ finally reads:

$$\begin{eqnarray*}H_0&=&{1\over2}\left(P_\xi^2+P_\eta^2+P_\zeta^2\right)\\ &&-k\... ... P_\xi\right) -{a\over2}\xi^2-{b\over2}\eta^2-{c\over2}\zeta^2, \end{eqnarray*}$$

where:

$$a={1\over2}\left[1-{3\over2}~\delta+3J\left({1 +\gamma\over2}\right)\left(3-{5\over2}~\delta\right)\right],$$

$$b={1\over2}\left[1+{3\over2}~\delta+3J\left({1 +\gamma\over2}\right)\left(3+{5\over2}~\delta\right)\right],$$

$$c=-1-9J\left({1+\gamma\over2}\right)\cdot$$

Then, let us focus on the linear system deriving from the quadratic Hamiltonian H₀. The eigenvalues are computed solving:

$$\lambda^4+\left(2k^2-a-b\right)\lambda^2+ab+k^2a+k^2b+k^4=0,$$

$$\lambda^2-c=0,$$

leading to imaginary solutions $\pm~\rm i~\omega_1,\pm~i~\omega_2, \pm~i~\omega_3$. The zero-order frequencies of the motion are:

$$\omega_1=\sqrt{2k^2-a-b+d\over2}\sim1,$$

$$\omega_2=\sqrt{2k^2-a-b-d\over2}\sim\sqrt{\mu},$$

$$\omega_3=\sqrt{-c}\sim1,\quad d=\sqrt{\left(a-b\right)^2-8k^2\left(a+b\right)},$$

and H₀ can be changed into:

$$H_0={1\over2}\left(P_1^2+\omega_1^2Q_1^2\right) -{1\over2}\l... ...a_2^2Q_2^2\right) +{1\over2}\left(P_3^2+\omega_3^2Q_3^2\right)$$

by a canonical transformation of the form:

$$\xi=a_{13}P_1+a_{14}P_2,\quad P_\xi=a_{31}Q_1+a_{32}Q_2,$$

$$\eta=a_{21}Q_1+a_{22}Q_2,\quad P_\eta=a_{43}P_1+a_{44}P_2,$$

$$\zeta=Q_3,\quad P_\zeta=P_3.$$

The required expression for H₀ and the canonicity of the transformation provide ten relations leading to:

$$a_{13}={1\over\sqrt{v^2-2kv-a}},\quad a_{14}={1\over\sqrt{a+2ku-u^2}},$$

$$a_{21}={-\omega_1\over\sqrt{u^2+2ku-b}},\quad a_{22}={\omega_2\over\sqrt{b-2kv-v^2}},$$

$$a_{31}=ua_{21},\quad a_{32}=va_{22},\quad a_{43}=va_{13},\quad a_{44}=ua_{14},$$

where:

$$u={b-a+d\over4k},\quad v={b-a-d\over4k}\cdot$$

Finally, the canonical transformation:

$$Q_1=\sqrt{2I_1\over\omega_1}\sin\varphi_1,\quad P_1=\sqrt{2I_1\omega_1}\cos\varphi_1,$$

$$Q_2=\sqrt{2I_2\over\omega_2}\sin\varphi_2,\quad P_2=\sqrt{2I_2\omega_2}\cos\varphi_2,$$

$$Q_3=\sqrt{2I_3\over\omega_3}\sin\varphi_3,\quad P_3=\sqrt{2I_3\omega_3}\cos\varphi_3$$

leads to:

$$H_0=I_1\omega_1-I_2\omega_2+I_3\omega_3.$$

The zero-order Hamiltonian H₀ represents a triple harmonic oscillator under its normalized form, providing zero-order solutions in the vicinity of the equilibrium position. The resulting motion is composed of two elliptic motions in the $\xi\eta$-plane, one with a short period ( $2\pi/\omega_1$) and one with a long period ( $2\pi/\omega_2$), completed by an oscillation along the $\zeta$-axis with a short period ( $2\pi/\omega_3$).

4 Asymptotic expansions

Normalizing H₀, the $(I,\varphi)$-variables are suitable for performing a nth-order Lie transformation $(I,\varphi)\rightarrow(I^\prime,\varphi^\prime)$ on the approximation H_(n) of the Hamiltonian of the motion (Deprit 1969). First, the last two canonical transformations are applied to the components H₁, H₂, ..., H_n for dealing with $(I,\varphi)$-variables. Then, the Lie algorithm is carried out using an algebraic processor from the University of Namur (Moons 1991). It results in a nth-order normalized Hamiltonian depending only on the variables $(I^\prime)$, providing nth-order solutions in the vicinity of the equilibrium position.

The Lie algorithm was carried out up to the fourth order mainly to accurately compute the frequencies of the motion (they are only updated at even steps). It gives the theory a better reliability for prediction. Non-significant solution terms (with respect to the statistics of the solution) will be removed after the fit to the observations.

The development process is a bit less straightforward than in the case of the geosynchronous satellite because of extra terms that are not canceled by the orientation of the normal frame (Oberti 1994). To increase the convergence, each even generator has to be completed by polynomials in the action variables devoted to removing solution terms having the unfortunate propensity to grow too fast.

5 Some checks

The value of the theory in the dimensionless space, at a selected time $\tau$, depends on six internal parameters: three amplitudes I⁰₁, I⁰₂, I⁰₃, and three phases $\varphi^0_1$, $\varphi^0_2$, $\varphi^0_3$ (dropping the primes to save writing). When going back to the real world, another two (external) parameters can be fitted for fine tuning: the companion's longitude $\theta_0$ at $\tau=0$ (its initial value is provided by the companion's theory) and the scaling factor r₀ between the theoretical space and the actual one (its initial value is computed from the companion's observed period).

With an expansion around the equilibrium position, the accuracy of the theory strongly depends on the satellite's libration amplitude, and its ability to correctly represent the motion increases with the development order. To check the inner limits of the theory (due to the finite expansion of the Hamiltonian), the solution was fitted by least squares (internal parameters only) to a numerical integration of the same Hamiltonian system over twenty years.

With initial values causing libration amplitudes of 10, 20, 30, 40, and 50 degrees, the accuracy is of order 10^-6, 10^-4, 10^-3, 10^-2, and $5\times10^{-2}$ arcsec. Helene's libration amplitude being of order 30 degrees, and Telesto and Calypso's much smaller, the theory qualifies for the job regarding this test.

The theory assumes a non-Keplerian circular motion for the second primary and doesn't take into account the neighboring satellites. To check the outer limits of the theory (due to this more complex environment), the solution was fitted by least squares to simulated observations generated at each observation time by a numerical integration of the studied satellite. The numerical model, including the first six Saturn's satellites (Vienne & Duriez 1995; Duriez & Vienne 1997) was previously fitted to the actual observations.

When only the six internal parameters are fitted, the precision strongly depends on the initial time chosen for the series through the induced companion's initial position. Its average value is of order 0.15 arcsec for the three satellites. Allowing the fit of the companion's initial longitude removes this dependency, and the precision is about 0.03 arcsec for Helene, 0.05 arcsec for Telesto, and 0.08 arcsec for Calypso. The effect of the last parameter (the scaling factor) is barely visible. Despite a greater libration amplitude, Helene's motion seems to be a bit more accurately modeled, possibly because the motion of Tethys is less regular than Dione's.

Here again, the theory pass the test for the three satellites.

The last test shows the relevance of the model: without or with Saturn's oblateness (I_J=0 or 1) and bidimensional or tridimensional (I_D=2 or 3). To remove any side effects, the companion's longitude and the scaling factor were fitted in every case. The resulting precisions are summarized in Table 1, in arcsec, for the different values of I_J/I_D. Comparatively to the basic model (0/2), only the combination of the two extensions seems to make the difference...

   

Table 1: Relevance of the model I_J/I_D (arcsec).

Model	Helene	Telesto	Calypso
0/2	0.20	0.60	0.75
0/3	0.15	0.40	0.65
1/2	0.10	0.55	0.60
1/3	0.03	0.05	0.08

6 Fits to the observations

They are performed using the observations from Reitsema (1981a, 1981b) obtained in 1980 and 1981 for Helene (20 obs.), Telesto (2 obs.), and Calypso (4 obs.), the observations from Oberti et al. (1989) obtained in 1981, 1982, 1984, and 1985 for Helene (133 obs.), Telesto (46 obs.), and Calypso (93 obs.), the observations from Veiga & Vieira Martins (2000) obtained in 1985 and 1987 for Helene (22 obs.), and the observations from Veiga et al. (2002) obtained in 1995 and 1996 for Helene (46 obs.), Telesto (13 obs.), and Calypso (19 obs.), the latter being deduced from the positions of the major satellites as exposed in Vienne et al. (2001).

For each satellite, three different fits were performed using first the observations from 1980 to 1987, then the observations from 1995 to 1996, and finally all the observations. The precisions of the resulting solutions are summarized in Table 2, in arcsec. The residuals from Veiga et al. (2002) were 0.15 arcsec for Helene, 1.43 arcsec for Telesto, and 0.42 arcsec for Calypso.

The selected solutions are the ones computed from the whole set of observations. Deduced from these solutions, the respective libration angles for Helene, Telesto, and Calypso are of order 28.5, 4.5, and 10.5 degrees, and their respective long-period libration periods are of order 783, 662.5, and 663.5 days.

   

Table 2: Results of the fits to the observations (arcsec).

Data	Helene	Telesto	Calypso
1980-1987	$0.27\pm0.27$	$0.28\pm0.26$	$0.25\pm0.22$
1995-1996	$0.13\pm0.19$	$0.40\pm0.45$	$0.25\pm0.29$
1980-1996	$0.65\pm0.79$	$0.36\pm0.35$	$0.83\pm1.03$

  

Table 3: Helene's frequencies (rd/day) and phases (rd).

  

Table 4: Helene's positions (a.u.) and velocities (a.u./day).

7 Back to physical variables

For each satellite, the fitted solution gives the Cartesian coordinates and their conjugate variables in the theoretical space. They are developed in trigonometric polynomials involving three time-depending angles $\varphi_1$, $\varphi_2$, $\varphi_3$:

$$\varphi_1=\omega_1\tau+\varphi^0_1,\quad \varphi_2=\omega_2\tau+\varphi^0_2,\quad \varphi_3=\omega_3\tau+\varphi^0_3,$$

where $\tau$ is the dimensionless time computed from the Saturnian Julian day t (Terrestrial Time minus light time) by:

$$\tau={2\pi\over k}\left({t-t_0\over P}\right),$$

with t₀=2 451 545.0, P being the companion's observed period (in days): 1.887803 for Tethys and 2.736916 for Dione (Vienne & Duriez 1995).

Obtaining the saturnicentric positions X, Y, Z and velocities $\dot X$, $\dot Y$, $\dot Z$ (where dots stand for derivation with respect to the real time) referred to the mean ecliptic and equinox of J2000 requires several steps. With the notations of the first canonical transformation:

$$\left(\matrix{X\cr\noalign{\smallskip } Y\cr\noalign{\smallsk... ...gn{\smallskip } y\cr\noalign{\smallskip } z\cr}\right)\right],$$

$$\left(\matrix{\dot X\cr\noalign{\smallskip } \dot Y\cr\noalig... ...(1-\gamma\right)\cr\noalign{\smallskip } P_z\cr}\right)\right],$$

where r₀ is the scaling factor, R₁ is a time-depending rotation of angle $\theta$ in the ring plane devoted to returning to the fixed frame (and adding a new trigonometric variable in the series), where:

$$\theta=k\tau+\theta_0,$$

and R₂ is a composition of numerical rotations linking the ring plane to the chosen frame.

The signs of y and P_x are changed for Calypso to take into account its position near L₅.

Then, the series are shortened giving the three amplitudes I⁰₁, I⁰₂, I⁰₃ their numerical values and removing non-significant terms, and the frequencies of the four angles are multiplied by $2\pi/kP$ (and called $\nu_1$, $\nu_2$, $\nu_3$, and $\lambda$) to deal directly with a more realistic time.

The signs of $\nu_1$, $\nu_2$, and $\nu_3$ are changed for Calypso to take into account its position near L₅.

8 User's manual

The frequencies (rd/day) and phases (rd) for Helene, Telesto, and Calypso, together with their saturnicentric positions (a.u.) and velocities (a.u./day), referred to the mean ecliptic and equinox of J2000, are summarized in Tables 3-8.

Each coordinate is the sum of the corresponding numerical coefficients multiplied by the cosine (f=0) or sine (f=1) of the angle $\varphi$ computed by:

$$\begin{eqnarray*}\varphi=n_1(\nu_1T+\varphi^0_1)+n_2(\nu_2T+\varphi^0_2)+n_3(\nu_3T+\varphi^0_3)+\ell(\lambda T+\theta_0), \end{eqnarray*}$$

where T is a Saturnian time (Terrestrial Time minus light time) in days computed from 2 451 545.0.

A downloadable FORTRAN code for the theory is available at: ftp://www.imcce.fr/pub/ephem/satel/htc20/

  

Table 5: Telesto's frequencies (rd/day) and phases (rd).

  

Table 6: Telesto's positions (a.u.) and velocities (a.u./day).

  

Table 7: Calypso's frequencies (rd/day) and phases (rd).

  

Table 8: Calypso's positions (a.u.) and velocities (a.u./day).

9 Conclusion

The comparison with numerical simulations shows the model's efficiency to deal with Lagrangian satellites similar to Tethys and Dione's companions (with small and moderate libration amplitudes). With almost twenty years of observations, the precision obtained when fitting each satellite's theory gives a rather good level of confidence for predictions.

References

 

• Campbell, J. K., & Anderson, J. D. 1989, AJ, 97, 1485 In the text NASA ADS  
• Deprit, A. 1969, Celest. Mech., 1, 12 In the text NASA ADS  
• Deprit, A., & Rom, A. 1970, A&A, 5, 416 In the text NASA ADS  
• Duriez, L., & Vienne, A. 1997, A&A, 324, 366 In the text NASA ADS  
• Moons, M. 1991, Dept. of Math., Facultés Universitaires Notre Dame de la Paix, Namur, Belgium In the text  
• Oberti, P., Veillet, C., & Catullo, V. 1989, A&AS, 80, 289 In the text NASA ADS  
• Oberti, P. 1990, A&A, 228, 275 In the text NASA ADS  
• Oberti, P. 1994, A&A, 284, 281 In the text NASA ADS  
• Reitsema, H. J. 1981a, Icarus, 48, 23 In the text NASA ADS  
• Reitsema, H. J. 1981b, Icarus, 48, 140 In the text NASA ADS  
• Sharma, R. K., & Subba Rao, P. V. 1976, Celest. Mech., 13, 137 In the text NASA ADS  
• Veiga, C. H., & Vieira Martins, R. 2000, A&AS, 143, 405 In the text NASA ADS  
• Veiga, C. H., Vieira Martins, R., Vienne, A., Thuillot, W., & Arlot, J.-E. 2002, preprint In the text  
• Vienne, A., & Duriez, L. 1995, A&A, 297, 588 In the text NASA ADS  
• Vienne, A., Thuillot, W., Veiga, C. H., Arlot, J.-E., & Vieira Martins, R. 2001, A&A, 380, 727 In the text NASA ADS  
• Szebehely, V. 1967, Theory of orbits (Academic Press) In the text