A&A 449, 821-826 (2006)
DOI: 10.1051/0004-6361:20054451

Orbital resonance widths in a uniformly rotating second degree and order gravity field

Ø. Olsen

Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway

Received 31 October 2005 / Accepted 6 December 2005

Abstract
Aims. In this work I seek to derive analytical expressions for the widths of mean motion resonances in uniformly rotating second degree and order gravity fields.
Methods. A Deprit perturbation series is used to derive Hamiltonians valid close to mean motion resonances. These expressions are reduced to the Hamiltonian of a pendulum through canonical transformations.
Results. Analytical expressions for the orbital resonance widths in uniformly rotating second degree and order gravity field are derived and tested against numerical simulations. The widths of these mean motion resonances are independent of the rotation rate and the mass of the central body. The general stability of retrograde orbits are explained by the overlap criteria.

Key words: celestial mechanics - minor planets, asteroids

1 Introduction

It has been shown that resonances between a rotating body and mean orbital motion may play a strong role (Hu & Scheeres 2004) in determining the stability of the orbits. The Hamiltonian of an orbiter close to a resonance was derived by Olsen (2001). The object of this paper is to study the widths of these resonances both numerically and analytically. Section 2 gives a short recap of the perturbation theory used to derive the resonant Hamiltonian, Sect. 3 is used to derive analytically the width of the resonances while Sect. 4 compares the analytical results with numerical experiments.

   
2 Motion close to a resonance

I first state the potential of the problem. Then I show how to express the Hamiltonian in the uniformly rotating body-fixed frame and give a series expansion for this potential in the Delaunay elements. Following this, I use a perturbation theory to study orbits without any assumptions about resonance, and I explain why many orbits are shape stable. Finally, I derive a simple analytical model for the trajectories close to resonances.

2.1 The potential

The coordinates are body-fixed such that the principal moments of inertia are ordered as $I_{xx} \leq I_{yy} \leq I_{zz}$ with axes x, yand z. For simplicity I assume the body to be a tri-axial ellipsoid. Then the second degree potential (Heiskanen & Moritz 1967) is

$$V = -\frac{\mu}{r} + V_1,$$ (1)

where

$$V_1 = -\frac{\mu}{r}\left( \frac{R^2 C_{20}}{2r^2}\left(3\sin... ...ght) - \frac{3 R^2 C_{22}}{r^2}\cos^2\phi \cos 2\theta\right),$$ (2)

where $\phi$ is the latitude, $\theta$ the longitude and R is the equatorial radius of the boy. The orbits will be analyzed as two-body orbits perturbed by the potential V₁.

2.2 The Hamiltonian

The Hamiltonian of the two-body problem expressed in terms of the Delaunay variables L, G, H, l, g and h(Murray & Dermott 1999) is

$$F_0 = -\frac{\mu^2}{2L^2}\cdot$$ (3)

I assume the body to rotate uniformly around the z-axis. The sign of the angular velocity, $\omega_{\rm T}$, is given according to the right-hand rule. The only difference between the orbital elements of a two-body problem in body-fixed orbital elements and inertial orbital elements is the argument of the ascending node. The longitude of the ascending node, $h = -\partial F_0/\partial H$, is decreasing with a rate equal to $\omega_{\rm T}$. Therefore, the Hamiltonian can be expressed in the body fixed elements as

$$F_0 = -\frac{\mu^2}{2L^2} -\omega_{\rm T} H.$$ (4)

The transformation to body-fixed elements has been achieved by subtracting $\omega_{\rm T}$ times the z-component of the angular momentum from the original Hamiltonian (Hand & Finch 1998). F₀ denotes he zeroth order Hamiltonian, and its exact solution is known.

The orbits are perturbed by a first order Hamiltonian, F₁ = V₁, which also must be expressed in terms of the Delaunay variables. The relations

   

$$\sin \phi = \sin i \sin(g+f),$$ (5)

$$x = r(\cos(g+f)\cos h - \sin(g+f)\sin h \cos i)$$

$$y = r(\cos(g+f)\sin h + \sin(g+f)\cos h\cos i)$$ (6)

$$z = r\sin(g+f)\sin i,$$

$$\frac{a^3}{r^3} = \mathcal A_0+\sum_{n=1}^\infty \mathcal A_n\cos nl$$

$$\frac{a^3}{r^3}\cos 2f = \sum_{n=1}^\infty \mathcal B_n\cos nl$$ (7)

$$\frac{a^3}{r^3}\sin 2f = \sum_{n=1}^\infty \mathcal C_n\sin nl$$

and the definitions of the Delaunay variables yield

$$-\frac{\mu^4 R^2C_{20}}{4L^3G^3} \left(1-\frac{3H^2}{G^2}\right)$$                            F_C20 =

$$- \frac{\mu^4 R^2C_{20}}{4G^2L^6}\left(\sum_{n=1}^{\infty} (G^2-3H) \mathcal A_n\cos nl\right.$$ (8)

$$- \sum_{n=1}^{\infty} 3(G^2-H^2)\cos(2g-nl)(\mathcal B_n-\mathcal C_n)$$

$$- \left.\sum_{n=1}^{\infty} 3(G^2-H^2)\cos(2g+n~l)(\mathcal B_n+\mathcal C_n)\right)$$

and

 

$$- \frac{3\mu^4 R^2C_{22}}{8G^2L^6} \biggl(4(G^2-H^2)\mathcal A_0\cos2h$$                                   F_C22 =

$$+ \sum_{n=1}^{\infty} 2(G^2-H^2) (\cos(2h-nl)+ \cos(2h-nl))\mathcal A_n$$

$$+ \sum_{n=1}^{\infty} \left((G-H)^2\cos(2g-2h-nl)\right.$$

$$+ \left.(G+H)^2 \cos(2g+2h-nl)\right) (\mathcal B_n-\mathcal C_n)$$ (9)

$$+ \left.\sum_{n=1}^{\infty} \left( (G-H)^2 \cos(2g-2h+nl)\right.\right.$$

$$+ \left. (G+H)^2\cos(2g+2h+nl)\right) (\mathcal B_n+\mathcal C_n)\bigg),$$

where F₁ = F_C20 + F_C22. A theory to compute these Fourier coefficients have been developed by Hagihara (1970). Olsen (2001) has given coefficients up 8th order in the eccentricity.

   
2.3 Canonical perturbation theory

In this paper I'll use a Deprit perturbation series (Lichtenberg & Lieberman 1992) to first order. I first introduce a time $\bar t$ associated with the reference orbit through the equations

$$\frac{\partial (L, G, H)}{\partial \bar t} = -\frac{\partial F_0}{\partial (l, g, h)}$$ (10)

and

$$\frac{\partial (l, g, h)}{\partial \bar t} = \frac{\partial F_0}{\partial (L, G, H)}\cdot$$ (11)

One possible algorithm is

 

$$\bar F_0$$ = F₀

$$\bar F_1 \lim_{t\to\infty} \frac{1}{T}\int_0^T F_1 {\rm d}\bar t$$ = (12)

$$-\int F_1 {\rm d}t.$$ w₁ = (13)

The new Hamiltonian is

$$\bar F = \bar F_0 + \bar F_1 = -\frac{\mu^2}{2\bar L^2} -\ome... ...L^3 \bar G^3} \left(1-\frac{3 \bar H^2}{\bar G^2}\right)\cdot$$ (14)

The transformations between the mean elements and the instantaneous elements are given by

$$x = \bar x + \{w_1, \bar x\}$$ (15)

and

$$\bar x = x - \{w_1, x\},$$ (16)

where the brackets denotes the Poisson brackets. These transformations contain only periodic terms. Furthermore, singularities exist whenever there is an integer n such that

$$n\frac{\mu^2}{L^3} = 2\omega_{\rm T},$$ (17)

because of Eqs. (4), (10) and (12). I define these frequencies as n'th degree primary resonances, since they emerge in a first order perturbation theory from the nth Fourier component.

Sufficiently far from resonances, the orbits are on average similar to orbits in gravity fields without a second degree field. In other words the orbits' size and shape are constant on average. This theory confirms the results of Hu & Scheeres (2001). The orientation parameters $\Omega$ and $\omega$ changes linearly with time, while the inclination is constant on average. Before I can procede to define sufficiently far from resonances, I must derive a Hamiltonian valid at a nth degree primary resonance.

2.4 Adiabatic ordering

At a nth degree primary resonance, the argument n l + 2 h is expected to change slowly compared to other combinations of l and h. The canonical transformation given by

$$S_2 = (n~ l + 2 h)\tilde L + g~\tilde G + h~\tilde H,$$ (18)

introduces a new set of variables through

$$\tilde l\, \frac{\partial S_2}{\partial{\tilde L}} = n l + 2h$$ =

$$\tilde g \frac{\partial S_2}{\partial{\tilde G}} = g$$ =

$$\tilde h \frac{\partial S_2}{\partial{\tilde H}} = h$$ = (19)

$$\frac{\partial S_2}{\partial l} = n\tilde L$$ L =

$$\frac{\partial S_2}{\partial g} = \tilde G$$ G =

$$\frac{\partial S_2}{\partial h} = 2\tilde L + \tilde H.$$ H =

Both $\tilde l\,$ and $\tilde g\,$ are slowly changing at the resonance. I now rewrite the function w₁ as w₁₀+w₁₁ and assume w₁₁ to be much smaller than w₁₀. The time derivative of $\tilde h\,$ is by definition much larger than the time derivative of $\tilde l\,$. Hence, the averaging of the Hamiltonian can be performed over the angular variable $\tilde h$ instead of a reference time (Lichtenberg & Lieberman 1992). The results are

$$\bar F_0 = \tilde F_0$$

$$\bar F_{10} = \tilde F_{1s} = \frac{1}{2\pi}\int_0^{2\pi}\tilde F_1 {\rm d}\tilde h$$ (20)

$$w_{10} = \frac{1}{\omega_{\rm T}}\int (\tilde F_1 -\tilde F_{1s}) {\rm d}\tilde h$$ (21)

$$\bar F_{11} = 0$$ (22)

$$w_{11} = \frac{1}{\omega_{\rm T}}\int \frac{\partial \tilde F_0}{\partial \tilde L} \frac{\partial w_{10}}{\partial \tilde l} {\rm d}\tilde h$$ (23)

and

 

$$-\frac{\mu^2}{2 L^2} - \omega_{\rm T} H - \frac{(G^2 - 3H^2)} {4G^5 L^3} \mu^4 R^2 C_{20}$$                                              F =

$$- \frac{3(G + H)^2}{8G^2 L^6}\mu^4 R^2 C_{22} (\mathcal B_n + \mathcal C_n) \cos(n l + 2 g + 2 h)$$ (24)

$$- \frac{3(G^2 - H^2)}{4G^2 L^6} \mu^4 R^2 C_{22} \mathcal A_n \cos(n l + 2 h)$$

$$- \frac{3(-G + H)^2}{8 G^2 L^6} \mu^4 R^2 C_{22} (\mathcal B_n - \mathcal C_n) \cos(n l - 2 g + 2 h) .$$

In the last expression I have reversed the canonical transformation given by Eq. (18) and dropped any superscripts. Implicitly I have also assumed $\omega_{\rm T} \gg \vert\partial F/\partial G\vert$.

   
3 The widths of the resonances

The magnitude of the different terms in Eq. (24) can be estimated from the relations

$$\mathcal B_n - \mathcal C_n \propto e^2 \mathcal A_n \propto e^4( \mathcal B_n + \mathcal C_n).$$ (25)

A reasonable approximation is to ignore the terms with the factors $\mathcal B_n - \mathcal C_n$ and $\mathcal A_n$. The canonical transformation

$$S_2 = (n~ l + 2 g + 2 h)\hat L + g~\hat G + h~\hat H,$$ (26)

then reduces the number of angular variables to one. The new variables are given by the following set of equations

 

$$\hat l\, \frac{\partial S_2}{\partial{\hat L}} = n l+2 g+2 h$$ =

$$\hat g \frac{\partial S_2}{\partial{\hat G}} = g$$ =

$$\hat h \frac{\partial S_2}{\partial{\hat H}} = h$$ = (27)

$$\frac{\partial S_2}{\partial{ l}} = n\hat L$$ L =

$$\frac{\partial S_2}{\partial{ g}} = 2\hat L + \hat G$$ G =

$$\frac{\partial S_2}{\partial{ h}} = 2\hat L + \hat H,$$ H =

and it follows that the new momenta $\hat G$ and $\hat H$ are constants and that the new coordinates $\hat g$ and $\hat h$ vary linearly with time. To study the evolution of $\hat l\,$ and $\hat L$, the Hamiltonian can be rewritten as

$$\hat{H} = -\frac{\mu^2}{2 n^2\hat{L}^2} - \omega_{\rm T}(\hat H + 2\hat L) - f_{C20}(\hat L) - f_{C22}(\hat L) \cos(\hat l),$$ (28)

where the functions $f_{C20}(\hat L)$ and $f_{C22}(\hat L)$ are introduced to make the notation more concise. This Hamiltonian was derived by assuming the orbit is close to a resonance. Therefore, for a given eccentricity and inclination, the point L₀ given by

$$\frac{\mu^2}{n^2\hat L_0^3} -2\omega_{\rm T} - f_{C20}'(\hat L_0) = 0,$$ (29)

may be used to expand the Hamiltonian in a Taylor series:

 

$$\hat{H} \approx -2\omega_{\rm T}(\hat H + \hat L_0 ) - \frac{\mu^2}{2n^2\hat L_0^2} - f_{C20}(\hat L_0) - f_{C22}(\hat L_0) \cos(\hat l)$$

$$+ \Delta \hat L \left(\frac{\mu^2}{n^2\hat L_0^3} -2 \omega_{\rm T} - f_{C20}'(\hat L_0) - f_{C22}'(\hat L_0) \cos(\hat l)\right)$$ (30)

$$+ \frac{\Delta \hat L^2}{2} \left( -\frac{3\mu^2}{n^2\hat L_0^4} - f_{C20}''(\hat L_0) - f_{C22}''(\hat L_0) \cos(\hat l)\right) + ...$$

The first three terms are constants and can be removed form the Hamiltonian without changing the equations of motion. Assuming any combination of $\Delta\hat L$ and C₂₀ or C₂₂ to be small and using the definition of L₀, the Hamiltonian has been reduced to the Hamiltonian of a pendulum:

$$\hat{H} = -\frac{1}{2}\frac{3\mu^2}{n^2\hat L_0^4}\Delta \hat L^2 - f_{C22}(\hat L_0) \cos(\hat l).$$ (31)

The width of the separatrix is

$$\Delta\hat L_{\max} = 2\left(\frac{n^2 L_0^4~ f_{C22}(\hat L_0)}{3\mu^2}\right)^{1/2},$$ (32)

or in average Keplerian elements

$$\Delta {a} = \cos^2\left(\frac{ i}{2}\right) \sqrt{8 R^2~C_{22}~(\mathcal B_n + \mathcal C_n)}.$$ (33)

The separatrix of a pendulum gives the boundary between rotational and oscillatory motion. A natural definition for the width of resonance is therefore given by Eq. (33).

Stochastic regions exist around the separatrix (Lichtenberg & Lieberman 1992) and global stochasticity may be expected when separatrices from different resonances start to overlap. The stability of retrograde orbits are therefore explained by the oscillation width's dependence on the inclination. Also, orbits around slow rotators are also expected to be more stable than orbits around similarly shaped fast rotators, since the resonances are further apart and not because they are weaker. A resonance's strength depends neither on the semi-major axis of the orbit nor on the angular rotational velocity of the central object.

From Eq. (33) it is evident that one of the ignored terms in Eq. (24) may dominate the Hamiltonian at certain inclinations. Using the same procedure for these two terms, the two different expressions

$$\Delta {a} = \vert\sin i\vert\sqrt{4 R^2~C_{22}~\mathcal A_n}$$ (34)

and

$$\Delta {a} = \sin^2\left(\frac{ i}{2}\right) \sqrt{8 R^2~C_{22}~(\mathcal B_n - \mathcal C_n)}$$ (35)

are derived for the oscillation width. The assumption that parts of the resonant Hamiltonian can be ignored, doesn't hold for all combinations of eccentricities and inclinations, but these formulae show the oscillation widths to be strongly dependent on both the inclination and eccentricity of the orbit. Therefore, they can still be used to explain the general effects of the resonances. Equation (25) also shows that the oscillations widths for retrograde orbits are smaller than those for prograde orbits.

   
4 Numerical tests of the resonance widths

I performed numerical experiments to test the resonance widths given by Eqs. (33)-(35). For numerical analyses in this paper I used the Burlisch Stoer (Press et al. 1996a) method to integrate the equations of motions in cartesian coordinates. The equations of motion were integrated in non-rotating coordinates, and they were derived from Eqs. (1) and (2) by replacing the angle $\theta$ with $\theta - \omega_{\rm T} t$.

Figure 1: Resonance widths for a set of prograde orbits in the gravity field of asteroid Castalia. The color coding depicts the relative change in the mean semi-major axis as a function of mean semi-major axis and eccentricity. The horizontal coordinate represents the semi-major axis in kilometers. Dark gray tones represent a decrease in the semi-major axis, whilst light tones represents an increase. The angular coordinates where chosen to select orbits close the stable fixed point of Eq. (31).

Open with DEXTER

	Figure 2: A closeup of the 6:2 resonance.
Open with DEXTER

The numerical experiments were performed by integrating a number of orbits with different initial semi-major axes and eccentricities. The initial values of the inclination, mean anomaly, argument of periapse and ascending node were constant within each experiment. Each orbit was integrated over a given time, and initial and final average values were computed for the semi-major axes and eccentricities. The averaging was done over a time corresponding to five times the period of the initial orbit.

The asteroid Castalia was used for the first three experiments. Initial circular equatorial orbits around Castalia were studied by Hu & Scheeres (2004). The gravitational constant, reference radius, rotational period and the coefficients C₂₀ and C₂₂ were taken from Scheeres et al. (1996), and the integration time was set to 2000 times the rotational period of the asteroid. Two of these experiments were performed on prograde orbits with an inclination of 60 degrees. The initial semi-major axis was varied from 1.0 km to 2.25 km using a step size of $2.5 \times 10^{-3}$, and the eccentricity was varied from 0 to 0.5 with a step size of $5.0 \times 10^{-3}$. The selected initial mean anomaly, argument of periapse and the ascending node were based on the pendulum Eq. (31). The stable fixed point is at $\hat l = \pi$ for an inclination of 60 degrees. The transformation from mean to instantaneous elements is periodic and with an amplitude proportional to the perturbing Hamiltonian. Any combination of l, g and h such that $n~ l + 2 g + 2 h = \pi$, will give a value $\hat l \approx \pi$. In the first two experiments I used $n~ l ~+~ 2 g ~+~ 2 h =\pi$ and n  l  +  2 g  +  2 h = 0. The results from these two experiments are displayed in Figs. 1 and 3. I repeated the first experiment with a smaller range and stepsize to produce a closeup of the 6:2 resonance, which is shown in Fig. 2.

The gray level indicate the relative change between the mean of the final and the initial semi-major axes. Light gray tones indicate an increase in the semi-major axis while dark tones indicate a decrease. The coordinates of each point are the numerically computed mean initial semi-major axis in kilometers and eccentricity for each orbit. The resonances' locations and widths were computed using Eqs. (29) and (33) and are drawn as black lines. These figures span the resonances 4:2 through 9:2 at respectively 1.26 km and 2.18 km. The white area at small semi-major axes and eccentricities are caused by the averaging process. Several orbits ended either up on an escape trajectory or passed below the reference radius. These orbits are not plotted.

A set of retrograde orbits was tested in the third experiment. Since retrograde orbits are more stable than prograde orbits, I chose to study orbits closer to the asteroid. The initial semi-major axis was varied from 0.55 km to 1.8, while the interval of eccentricities remained the same. I used an inclination of 150 degrees and $\hat l \approx \pi$. The results are shown in Fig. 4 together with the computed widths of the resonances given by Eqs. (33) and (34). Trajectories below the reference radius is the cause of the white area in the upper left part of Fig. 4.

	Figure 3: Resonance widths for orbits close to the unstable fixed point of Eq. (31) in the the gravity field of Castalia.
Open with DEXTER

There are two families of stable orbits. The first family consists of orbits which are not inside any resonance. These orbits are described by Eqs. (12) and (14), and the condition sufficiently far away from resonances is given by the formulas for resonances' widths and loactions. The second family of stable orbits are orbits close to the center of the resonance. These orbits require a combination of angular coordinates which put them close to the stable fixed point of the pendulum equation. They are the spikes within each resonance shown in Fig. 1. Orbits with $\hat l \approx \pi$ and which start at the left within a resonance, migrate to the right and vice versa. When $\hat l \approx 0$ the effect is reversed and resonances pushes the orbits away from the center of the resonances. Orbits in regions where the separatrices from different resonances overlap are unstable as expected from Lichtenberg & Lieberman (1992).

	Figure 4: Resonance widths for a set of retrograde orbits around Castalia. This figure covers the resonances 2:2 through 6:2.
Open with DEXTER

Reasonable strong resonances are present between every n:2 and n+1:2resonance. These are not explained by a first order perturbation theory, but a second order theory would involve an integral (Lichtenberg & Lieberman 1992) of the form

$$w_2 = \int 2\bar F_2 -\{w_1, \bar F_1 - F_1\}~{\rm d}t.$$ (36)

$\bar F_2$ would have been selected to remove any secular terms from the above integral. The poisson bracket introduces many new singularites at locations given by

$$(\pm m \pm n)\frac{\mu^2}{L^3} - k_i\omega_{\rm T},$$ (37)

where m and n are integers such that $(\pm m \pm n) >0$ and k_itakes the value 2 or 4. A new adiabatic ordering will give Hamiltonians valid close to these secondary resonances.

Figure 2 and the overlap criteria explain why most orbits within any resonance are unstable. As the secondary resonances emerge within a primary resonace, they soon start to overlap and destroy the possibility of stable orbits. The secondary resonances are proportional to R⁴/r⁶ as opposed to R²/r³ for the primary resonances. They are therefore less important for orbits further away from the asteroid, and regions covered by a single primary resonances will become increasingly stable.

The disagreement between theory and experiment close to the asteroid is caused mainly by the Taylor expansion of f_C22 in Eq. (30). It introduces the partial derivative of the eccentricity with respect to the L, which leads to negative powers of the eccentricity. f_C22' and f_C22'' may therefore be non-negligible even when combined with small excursions of $\Delta\hat L$. Retaining f_C22 will complicate the analysis and most likely introduce separatrices which are not symmetric around the stable fixed points. On the other hand, the computed locations are reasonable good as is demonstrated by the 2:2 resonance. Its nominal position given by Eq. (17), is 0.80 km. This compares to the computed location of 0.86 km at zero eccentricity.

I chose the asteroid Itokawa in the final numerical experiment. This asteroid is a slow rotator as opposed to Castalia, and it can be used to double check the resonance widths' independence of the rotational speed of the central body. Scheeres et al. (2004) provide up to date information for this asteroid. The initial inclination was 60 degrees and I used the stable fixed point to select initial angular coordinates. Figure 5 confirms the previous results and shows resonances from 2:2, synchronous orbits, and up. No stable orbits were found at the 1:2resonance using this selection algorithm.

	Figure 5: Resonance widths for the asteroid Itokawa.
Open with DEXTER

5 Conclusions

A study of mean motion resonances in uniformly rotating second degree and order gravity field is presented. The mean motion resonances are an important factor when studying the stability of orbits in such a gravity field. Surprisingly, the widths of these mean motion resonances are independent of the rotation rate and the mass of the central body. Using the overlap criteria, stability in slowly rotating gravity fields is explained by the distance between the resonances and not by the strength of a particular resonance. On the other hand, the widths are strongly dependent on both eccentricity an inclination with resonance widths for retrograde orbits being much smaller than those for prograde orbits. Finally, the evolution of orbits within a resonance is not fully understood and warrants a closer study. A second order perturbation theory is probably necessary.

References

• Hand, L. N., & Finch, J. D. 1998, Analytical Mechanics (Cambridge: Cambridge University Press) (In the text)
• Hagihara, Y. 1970, Celestial mechanics (The Massachusetts Institute of Technology Press), Vol. 1, Chap. 5.11 (In the text)
• Hu, W., & Scheeres, D. J. 2001, Periodic orbits in rotating 2nd degree and order gravity fields. Ninth International Space Conference of Pacific-basin Societies (ISCOPS), Pasadena, CA (In the text)
• Hu, W., & Scheeres, D. J. 2004, Planet. Space Sci., 52, 685 [NASA ADS] [CrossRef] (In the text)
• Lichtenberg, A. J., & Lieberman, M. A. 1992, Regular and Chaotic Dynamics, 2d ed. (New York: Springer-Verlag) (In the text)
• Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) (In the text)
• Olsen, Ø. 2001, An Analytical Study of Orbits Around Rotating Nonsymmetric Objects (Cand. Scient./Msc. thesis, University of Oslo, Norway) (In the text)
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1996, Numerical Recipes in Fortran 90, 2d ed. (Cambridge: Cambridge University Press) (In the text)
• Scheeres, D. J., Ostro, S. J., Hudson R. S., & Werner R. A. 1996, Icarus, 121, 67 [NASA ADS] [CrossRef] (In the text)
• Scheeres, D. J., Broschart, S. J., Ostro, S. J., & Beener, L. A. 2004, The Dynamical Environment About Asteroid 25143 Itokawa: Target of the Hayabusa Mission, AIAA/AAS Astrodynamics Specialist Conference and Exhibit (In the text)
• Heiskanen, W. A., & Moritz, M. 1967, Physical Geodesy (San Francicso: W.H. Freeman) (In the text)