2 Theoretical considerations

2.1 Definition of units

In all calculations below we use the following units. The unit of time is one Julian day (1 d), where 1 d is exactly 86400 s. The sidereal day is given by 86164.09054 s (JPL). The unit of length is 1 astronomical unit (1 AU = $1.49597870691 \times 10^{11}$ m). Moreover the gravitational constant is chosen as G=1. In this system of units, 1 Solar mass is = $2.9591221 \times 10^{-4}$.

2.2 Initial conditions and input data

Although the number of figures in numerical input data below is sometimes greater than the estimated accuracy, those are the input parameters used in our program. We present them here simply for reproducibility reasons.

All planetary masses are derived from the mass ratios listed at the JPL database (JPL). However, the Earth and the Moon are integrated separately. Their individual masses are therefore computed using the mass ratios of Sun/(Earth+Moon) and Moon/Earth. The mass ratios are listed in Table 1.

 

Table 1: IAU estimates of the mass ratios used in present work, see (JPL).

Sun/Mercury	6023600
Sun/Venus	408523.71
Sun/(Earth+Moon)	328900.56
Sun/(Mars system)	3098708
Sun/(Jupiter system)	1047.3486
Sun/(Saturn system)	3497.898
Sun/(Uranus system)	22902.98
Sun/(Neptune system)	19412.24
Sun/(Pluto system)	$1.35 \times 10^{8}$
Moon/Earth	0.012300034
Phobos/Mars	$1.683 \times 10^{-8}$
Deimos/Mars	$2.804 \times 10^{-9}$

Relative masses for Phobos and Deimos as well as the positions and velocities for all the bodies are obtained from "Horizons On-Line Ephemeris System'' (DE-406) (JPL). The reference frame is ICRF/J2000.0. All our simulations start at MJD 51600.5 (2000-Feb.-26 00:00).

The moments of inertia for the Earth are based on the JPL data (J_z-J_y)/J_z =0.0032737634, J_z/mR² =0.3307007, R=6378.136 km, $m=8.8877456\times 10^{-10}$(JPL). This gives that $J_{z} =5.3427328\times 10^{-19}$ and $J_{x} = J_{y} =5.3252420 \times 10^{-19}$. For Mars, the moments of inertia are: $J_{z} =1.8034092 \times 10^{-20}$ corresponding to J_z/mR² =0.3662, $m=9.5495922 \times 10^{-11}$ (JPL) and R=3397.2 km (Bouquillon & Souchay 1999); $J_{x} =1.7931150 \times 10^{-20}$ and $J_{y} =1.7943585 \times 10^{-20}$ from J_z/J_x =1.005741and J_z/J_y =1.005044 (Bouquillon & Souchay 1999).

The obliquity of the Earth is given by $\epsilon =23.4393$ degrees at MJD 51600.5 (JPL). The present obliquity of Mars is $\epsilon =25.189417$ degrees (Folkner et al. 1997). Finally, the spin angular velocities (denoted $\omega$) for Earth or Mars are listed in Sects. 2.4 and 2.7.

2.3 Treatment of the N-body problem

The N-body problem is here treated in a purely classical way. In order to keep the model simple we will not include relativistic effects, solar mass-loss or tides. These effects are studied and discussed in Appendix. Two independent engines to integrate the equations of motion are applied. First, we use the SWIFT library (Levison & Duncan 2000) which contains both the Mixed Variable Symplectic method (MVS) and the Bulirsch-Stoer method. Secondly, we have implemented a simple, but general algorithm for repeated Richardson extrapolation of Euler's method. The time step in this method is cut in half an arbitrary number of times (usually five) for each time interval $\triangle t$ and then extrapolations are repeatedly carried out in order to obtain very good positions and velocities (analogous to Romberg's method). This method is applied for the tests in the Appendix.

We have written an ANSI C program where all the above numerical methods can be used. The stability of the numerical results are tested using all methods. For appropriate time steps they all yield identical results (e.g. $\sim$15 figures in the positions are the same after thousands of years of simulation). For the final simulations, we decided to apply the MVS method because of its high speed and high stability properties (e.g. stability of the total angular momentum).

The program tests that the total mechanical energy and angular momentum are conserved at every time step. If conservation is violated the particle system is moved back to the previous time step. The integration is then attempted using half the previous time step. This procedure may then be repeated until conservation is fulfilled. The simulation then continues a certain number of steps (user defined) using the smaller time step before restoring the original time step. Another check of the accuracy is to study rapid oscillations (such as moon orbits) for peculiarities, since those are most sensitive to numerical errors. Standard modifications in order to gain numerical stability are also applied in the code, e.g. summation order. Trivial checks such as that compiler optimizations do not alter the results are also done. For every computed output quantity $\chi$, we check that for two distinct simulations the numerical truncation error $\varepsilon$is negligible: $\varepsilon =\left\vert \chi (\triangle t)-\chi (\frac{\triangle t}{2})\right\vert$. We have also investigated the sensitivity in obliquity (output) to variation in some input parameters by computing the condition number C = (relative change in output)/(relative change in input). Variations in solar mass, Earth mass, Moon mass, moments of inertia and spin angular velocity typically result in C= 10^-4-10^-2. This indicates that the problem is well conditioned. In present computations, we use the time step $\triangle t =0.182625$ d or lower (as low as 0.0036525 d in extreme cases). Finally, the consistency of our results is also checked by comparisons with measurements and findings of other authors, whenever possible.

2.4 The instantaneous behavior of the spin axis of a rigid body - rotation of the Earth

Instantaneous studies of spin axis precession, nutation and obliquity require a method for treating the momentary dynamics of a rotating rigid body. The approach below is general for rigid bodies and in principle exact; it can be used to study oscillations at any time scale and is still convenient for numerical simulations. We use a center of mass reference system S(x,y,z), i.e. a system that always follows the center of mass G of the body and has its axes fixed in inertial space. These axes are parallel with the external heliocentric system. $S^{\prime }(x^{\prime },y^{\prime },z^{\prime })$ is also a reference system with the origin at the center of mass, but its axes are instead fixed in the precessing body. These axes coincide with the body's principal axes so the inertia tensor $\overline{\overline{J}}_{_{\rm G}}(S^{\prime })$ is always diagonal. Furthermore, the body is considered to be a rigid (not necessarily homogeneous) oblate ellipsoid with principal moments of inertia $J_{x}\neq J_{y}\neq J_{z}$in the general case. The actual numerical values are taken from satellite measurements (see Sect. 2.2). We shall thus ignore the phenomenon of variation of latitude and assume that the spin axis is rigidly fixed along the z'-axis of S'. In the S-system, the total external torque ( $\overline{M}_{_{\rm G}}$) acting on the body w.r.t. the center of mass G is related to the angular momentum ( $\overline{L}_{_{\rm G}}$) according to:

$$\overline{M}_{_{\rm G}}=\frac{\rm d}{{\rm d}t}\overline{L}_{... ...J}}_{_{\rm G}}(S)\frac{\rm d\overline{\omega }}{{\rm d}t}\cdot$$ (1)

The similarity transformation between the S- and S'-system is given by:

$$\overline{\overline{J}}_{_{\rm G}}(S)=\overline{\overline{R}... ...\prime })\overline{\overline{R}}\equiv \overline{\overline{A}}$$ (2)

where $\overline{\overline{J}}_{_{\rm G}}(S^{\prime })$ remains diagonal at all times and $\overline{\overline{R}}$ is the Eulerian rotation matrix in the x-convention (Goldstein 1980):

$$\overline{\overline{R}}= \left[\! \begin{array}{ccc} \cos \p... ... \theta \cos \phi & \cos \theta \end{array}\!\right]\!\!\cdot$$ (3)

If we define ${\rm d}\overline{\overline{J}}_{_{\rm G}}(S)/{\rm d}t\equiv \overline{\overline{B}}$, the matrix elements of $\overline{\overline{B}}$ then depend on time according to:

$$B_{ij}=\frac{\rm d}{{\rm d}t}A_{ij}(\phi ,\theta ,\psi )= \le... ...partial }{\partial \psi }\right] A_{ij} (\phi ,\theta ,\psi ).$$ (4)

We thus obtain the following first order nonhomogeneous differential equation:

$$\overline{M}_{_{\rm G}}=\overline{\overline{B}}\, \overline{... ...ine{\omega }(t+{\rm d}t)-\overline{\omega }(t)}{{\rm d}t}\cdot$$ (5)

By replacing dt with $\triangle t$, we end up with the numerical expression:

$$\overline{\omega }(t+\triangle t)= \left[ \overline{\overline... ...\overline{\omega }(t)\right] \triangle t+\overline{\omega }(t)$$ (6)

which subsequently can be improved using a few Richardson extrapolations. The explicit expression within the brackets is too long and complicated to be given here; it is fully implemented in our C-program. Initial conditions for Earth are given for simulations forward in time as follows: $\overline{\omega }(0)=\omega \sin \epsilon \overline{e}_{y}+\omega \cos \epsilon \overline{e}_{z}$, where $\omega =2\pi \times 1.00273791$ radians/d. The second factor is the ratio between a Julian day and a sidereal day (JPL). The obliquity at MJD 51600.5 is $\epsilon =23.4393$ degrees (JPL). $\overline{\omega }(t)$will certainly at a later stage precess away from this direction but never away from the z'-axis. $\overline{\overline{A}}\, ^{^{-1}}(0)$ is determined by $\phi (0)=\pi ,\theta (0)=\epsilon ,\psi (0)=0$. In the case of $\overline{\overline{B}}(0)$, we set $\frac{\rm d\phi }{{\rm d}t}(0)=\frac{\rm d\theta }{{\rm d}t}(0)=0$ and $\frac{\rm d\psi }{{\rm d}t}(0)=\omega$. For simulations backward in time, the initial conditions that change are: $\overline{\omega }(0)\rightarrow -\overline{\omega }(0)$, $\frac{\rm d\psi }{{\rm d}t}(0)\rightarrow -\omega$ and $\phi (0)=0,\theta (0)=\pi -\epsilon$. $\overline{M}_{_{\rm G}}(0)$ also needs to be calculated at the start of the simulation (and also of course at every other time step). Note that an arbitrary vector must obey the transformation $\overline{u}=\overline{\overline{R}}\, ^{^{-1}}\overline{u}'$between the S and $S^{\prime }$-system. In  $S^{\prime }$, the basis vector $\overline{e}_{\omega }$ is always given by (0,0,1)while in S, $\overline{e}_{\omega }=(x,y,z)$. We therefore have the specific relation:

$$\overline{e}_{\omega }(t)=\left[ \begin{array}{c} x(t)\\ y(... ...a (t)\cos \phi (t)\\ \cos \theta (t) \end{array}\right] \cdot$$ (7)

This means that given $\overline{e}_{\omega }(t)$ in S, we also obtain $\theta (t)$ and $\phi (t)$. The last angle is simply given by $\psi (t)=\omega t$ in the case J_x=J_y, otherwise $\psi (t+\triangle t)=\psi (t)+\omega \triangle t$. Finally, the angular velocities in Eq. (4) needed to determine $\overline{\overline{B}}(t)$ are approximated with $\triangle \phi /\triangle t$and $\triangle \theta /\triangle t$, while $\frac{{\rm d}\psi }{{\rm d}t}=\omega (t)$.

2.5 External torques acting on the planet

The exact expression for external torques from particles j acting on an ellipsoid at time t may in the body fixed system $S^{\prime }$be written:

$$\overline{M}^{\prime }_{_{\rm G}}(t)=-G\sum _{j}m_{j}\int _{... ...^{\prime }+\overline{r}^{\prime }_{j}\right\vert ^{3}}{\rm d}m$$ (8)

where $\overline{r}^{\prime }$ is pointing at an arbitrary mass element dm in the ellipsoid and $\overline{r}^{\prime }_{j}=(x^{\prime }_{j},y^{\prime }_{j},z^{\prime }_{j})$is the radius vector to the particle j. However, this expression becomes far too complicated to compute, even for a spheroid. An extremely good approximation is due to MacCullagh (see Goldstein 1980, p. 228). In this approximation the potential energy of an oblate ellipsoid of mass m and a particle j is:

$$W(i,j)=-\frac{Gmm_{j}}{r_{j}}-\frac{Gm_{j}}{2r_{j}^{3}}\left( J_{x}+J_{y}+J_{z}-3I\right)$$ (9)

where $I=(J_{x}x^{\prime 2}_{j}+J_{y}y^{\prime 2}_{j}+J_{z}z^{\prime 2}_{j})/r_{j}^{2}$is the moment of inertia about the direction of $\overline{r}^{\prime }_{j}$. J_x, J_y and J_z are of course the moments of inertia w.r.t. the principal axes. The torque about the center of mass G of body i exerted by all particles j may now be written as:

$$\overline{M}^{\prime }_{_{\rm G}}(t)= 3G\sum _{j}\frac{m_{j}}... ...}z^{\prime }_{j},(J_{y}-J_{x})x^{\prime }_{j}y^{\prime }_{j}).$$ (10)

The instantaneous torques from all the other planets, the Moon and the Sun are included in the summation at every time step. In the case of Mars we are also including the moons Phobos and Deimos (see Sect. 2.10). Since $\overline{M}^{\prime }_{_{\rm G}}(t)$is computed in $S^{\prime }$, it is subsequently transformed back to S ( $\overline{M}_{_{\rm G}}(t)=\overline{\overline{R}}\, ^{^{-1}}\overline{M}'_{_{\rm G}}(t)$) so it can be used in Eq. (6).

2.6 The summer solstice for the northern hemisphere

Suppose $\overline{r}$ is the radius vector from the center of the Sun pointing at the Earth-Moon system's center of mass and $\overline{v}$ is its velocity. The normal vector to the ecliptic plane is then given by $\overline{n}=\overline{r}\times \overline{v}/\left\vert \overline{r}\times \overline{v}\right\vert$. Let $\overline{\omega }$ be the spin angular velocity of Earth. If $\left( \overline{r}\times \overline{\omega }\right) \cdot \overline{n}=0$is fulfilled we have either summer or winter solstice. It may then be checked if the angle $\varphi =\arccos \left( \overline{\omega }\cdot \overline{r}/\left\vert \overline{\omega }\right\vert \left\vert \overline{r}\right\vert \right) >90^{\circ }$in which case we have summer solstice. However, for simulations backward in time the condition for summer solstice instead becomes $\varphi <90^{\circ }$since $\overline{\omega }$ is then replaced by $-\overline{\omega }$. An analogous calculation can also be performed to find the northern summer solstices in the case of Mars.

2.7 Rotation of Mars

One day at Mars is slightly longer than the day for Earth, namely 24.622962 h compared to 23.93419 h (JPL). The angular velocity for Mars is therefore given by $\omega (0) =6.1241471$ radians/d. The obliquity (the angle between the Mars equator and its instantaneous orbital plane) of Mars is given by $\epsilon =25.189417$degrees (Mars pathfinder mission, Folkner et al. 1997). Currently, the summer solstice of the northern hemisphere of Mars occurs at Dec. 16, 2000, 04h UT. Given that $\overline{n},\, \overline{r}$ and $\epsilon$ are known, the equations

 

$$\overline{n}\cdot \overline{\omega } n\omega \cos (\epsilon )$$ = (11)

$$\overline{r}\cdot \overline{\omega } r\omega \cos (\epsilon +\frac{\pi }{2})$$ = (12)

$$\omega ^{2} \omega ^{2}_{x}+\omega ^{2}_{y}+\omega ^{2}_{z}$$ = (13)

completely determine the components of $\overline{\omega }$. Here, $\overline{n}$( $\overline{n}=\overline{r}\times \overline{v}$) is the normal vector to the orbital plane of Mars and $\overline{r}$ is the heliocentric position of Mars at its summer solstice. The obtained results in radians/d are:

 

$$\omega _{x}(0)=$$ (14)

$$\omega _{y}(0)=$$ (15)

$$\omega _{z}(0)=$$ (16)

The initial conditions for the Mars fixed system ( $S^{\prime }$) is $\phi (0)= 82.826460 ^{\circ }$. The Euler angle $\theta (0)$ is no longer directly related to the obliquity because the Mars orbital plane has an inclination of about $i=1.8504^{\circ }$relative to the ecliptic. This gives that the angle between $\overline{\omega }$(z'-axis) and the ecliptic normal (z-axis) becomes $\theta (0) =26.717957 ^{\circ }$. The last Eulerian angle is arbitrarily given by $\psi (0)=0$. The time derivatives are set analogously to what we did in Sect. 2.4. For simulations backwards in time the changes become as follows: $\overline{\omega }(0)\rightarrow -\overline{\omega }(0)$, $\frac{\rm d\psi }{{\rm d}t}(0)\rightarrow -\omega$ and $\phi (0)\rightarrow (\phi (0)-\pi ),\theta (0)\rightarrow (\pi -\theta (0))$.

2.8 The length of the summer half year

Given the semi major axis a, eccentricity e and the distance r at summer solstice for a specific year, we get from the definition of eccentric anomaly (E) that at summer solstice: $E_{\rm s}=\arccos ((a-r)/ae)$. The general relation between E and the true anomaly v is:

$$\tan \frac{E}{2}=\sqrt{\frac{1-e}{1+e}}\tan \frac{v}{2}\cdot$$ (17)

The specific true anomaly $\varphi$ from perihelion to the planet at summer solstice then becomes:

$$\varphi =2\arctan \left( \sqrt{\frac{1+e}{1-e}}\tan \left( \frac{E_{\rm s}}{2}\right) \right) \cdot$$ (18)

The angle from perihelion to the vernal equinox is given by $v_{\rm v}=\varphi -\pi /2$and for the autumnal equinox, $v_{\rm a}=\varphi +\pi /2$. The change in time dt is directly related to the change in eccentric anomaly dE through (Danby 1964):

$$n{\rm d}t={\rm d}E(1-e\cos E)$$ (19)

where n is the mean motion. Integrating between the equinoxes then gives that

$$nt=E_{\rm a}-e\sin E_{\rm a}-E_{\rm v}+e\sin E_{\rm v}$$ (20)

where $E_{\rm a}$ and $E_{\rm v}$ are the eccentric anomalies at autumnal and vernal equinoxes respectively. According to Eq. (17) we have:

 

$$E_{\rm a} 2\arctan \left( \sqrt{\frac{1-e}{1+e}}\tan \left( \frac{\varphi }{2}+\frac{\pi }{4}\right) \right)$$ = (21)

$$E_{\rm v} 2\arctan \left( \sqrt{\frac{1-e}{1+e}}\tan \left( \frac{\varphi }{2}-\frac{\pi }{4}\right) \right)\cdot$$ = (22)

The length of the summer half year relative to a full year is finally given by $nt/2\pi$. For more discussions regarding variations of the astronomical seasons, see e.g. Berger & Loutre (1994).

2.9 Spin axis precession and orbital precession

The precession angle of the spin axis ( $\phi _{\rm p}$) is calculated by

$$\phi _{\rm p}={\rm arccos}\left( \frac{(\overline{n}\times \... ...rline{n}_{0}\times \overline{\omega }_{0}\right\vert }\right)$$ (23)

where $\overline{n}$ and $\overline{n}_{0}$ are normal vectors to the instantaneous and present orbital plane, respectively, and analogously for $\overline{\omega }$ and $\overline{\omega }_{0}$.

The orbital precession is studied by calculating the Laplace-Runge-Lenz vector which is conserved for the two-particle problem. It is directed along the radius vector $\overline{r}$ to the perihelion point. In the solar system, many particle interactions (as well as the smaller general relativistic effect) cause the orbit of a planet of mass m to precess. The orbital precession can easily be studied by observing the rotating Laplace-Runge-Lenz vector which is given by:

$$\overline{A}=\overline{p}\times \overline{L}-\frac{GM_{\rm s}m^{2}}{r}\overline{r}$$ (24)

where $\overline{p}$ and $\overline{L}$ (of the planet) are the linear and angular momentum, respectively. An alternative method would be to compute the combined changes in the orbital parameters: the longitude of the ascending node $\Omega$, the argument of perihelion $\omega$and the inclination i: d $\omega /{\rm d}t+\cos (i){\rm d}\Omega /{\rm d}t$.

2.10 Torques due to Phobos and Deimos

Since Phobos and Deimos are so close to Mars (for Phobos the distance is less than 3 Mars radii), their motions are incredibly rapid compared to all other bodies in our simulation. For the time periods of interest, it is impossible to integrate these moons in the same manner as our own Moon. The objective here is therefore instead to obtain an average formula for the torque on Mars due to its moons. First a full integration will be carried out including Phobos and Deimos for a time span that is long enough to contain all relevant oscillations. The torques acting on Mars are calculated instantaneously, i.e. at every time step. We then choose circular moon orbits in the orbital plane of Mars (virtual moons) that yield exactly the same average torques as those obtained from the full integration. This treatment is similar to that given by Goldstein (1980, p. 229). Account will also be taken to the triaxiality of Mars, i.e., $J_{x}\neq J_{y}\neq J_{z}$. The expression for I in Eq. (9) may be rewritten as:

$$I=J_{x}\cos ^{2}\alpha +J_{y}\cos ^{2}\beta +J_{z}\cos ^{2}\gamma$$ (25)

where the angles $(\alpha ,\beta ,\gamma )$ to a moon-position $\overline{r}'$are measured relative to the body fixed system $S^{\prime }$. Let us now use a coordinate system xyz for which the effective mean orbit of a virtual moon is in the xy-plane and the spin axis of Mars ( $\overline{\omega }$) is in the xz-plane. The xy-plane is chosen as the instantaneous Mars orbit plane with the origin at the Mars center. The transformation between $S^{\prime }$ and this new xyz-system obeys: $\overline{r}'=\overline{\overline{R}}\, \overline{r}$, i.e. with $(\phi ,\theta ,\psi )=(\frac{\pi }{2},\theta ,\psi )$ one obtains:

$$\left[ \begin{array}{c} x^{\prime }\\ y^{\prime }\\ z^{\pr... ...ight] \left[ \begin{array}{c} x\\ y\\ 0 \end{array}\right] .$$ (26)

We treat a mean virtual moon orbit as a simple circular orbit which in polar coordinates becomes: $(x,y)=r(\cos \eta ,\sin \eta )$. We thus obtain:

$$\cos^{2}\alpha \equiv \left( \frac{x^{\prime }}{r}\right) ^{2... ...si \cos^{2}\eta\! -\!\frac{1}{2}\cos\theta \sin2\psi\sin2\eta.$$ (27)

As usual, the moon revolution ($\eta$) and the rotation of Mars ($\psi$) are extremely rapid compared to the change in $\theta$ allowing a time average over an integer number of periods to be performed (the symbol $\left\langle \, \, \right\rangle$refers to a time average, e.g. for a function f(t) this is given by $\left\langle f\right\rangle =\frac{1}{T}\int _{0}^{T}f(t){\rm d}t$):

$$\left\langle \cos^{2}\alpha \right\rangle =\left\langle \cos^... ...{2}\theta \left\langle \sin^{2}\psi \cos^{2}\eta \right\rangle$$

$$-\frac{1}{2}\cos\theta \left\langle \sin2\psi \sin2\eta \right\rangle =(\cos^{2}\theta +1)/4.$$ (28)

In the same way one obtains that $\left\langle \cos^{2}\beta \right\rangle =(\cos^{2}\theta +1)/4$and $\left\langle \cos^{2}\gamma \right\rangle =\frac{1}{2}\sin^{2}\theta$. The average moment of inertia can now be written:

$$\left\langle I\right\rangle \!=\!J_{x} (\cos^{2}\theta \!+\!... ...y} (\cos^{2}\theta +1)/4\!+\!J_{z}\frac{1}{2}\sin^{2}\theta .$$ (29)

The part of the averaged potential energy that has an orientational dependence is (see Eq. (9)):

 

$$\left\langle W\right\rangle \frac{Gm}{2r^{3}}\left( 3\left\langle I\right\rangle -(J_{x}+J_{y}+J_{z})\right)$$ =

$$\frac{Gm}{2r^{3}}(J_{z}-\frac{J_{x}+J_{y}}{2})\left( \frac{1}{2}-\frac{3}{2}\cos^{2}\theta \right) \cdot$$ = (30)

Finally the average torque acting on Mars due to a moon becomes:

$$\left\langle \overline{M}_{_{\rm G}}\right\rangle =\frac{3Gm... ...-\frac{J_{x}+J_{y}}{2}\right)\sin2\theta \overline{e}_{\xi } .$$ (31)

The torque is in the direction of $\overline{e}_{\xi }=\overline{n}\times \overline{\omega }/\left\vert \overline{n}\times \overline{\omega }\right\vert$, where $\overline{n}$ is the normal vector to the instantaneous Mars orbit and $\overline{\omega }$ is the spin axis of Mars. In order to apply this average formula at every time step we need to know the best possible value of r. Fortunately, we already have access to the accurate instantaneous Eq. (10). We therefore first start a simulation performing an exact integration of the Mars moons ( $\triangle t=1.82625 \times 10^{-3}$ d). We then divide the resulting instantaneous torques by the instantaneous $\sin2\theta$. By forming the average (k)

$$k=\sqrt{\left\langle \frac{M_{x}}{\sin2\theta }\right\rangle... ...{2}+\left\langle \frac{M_{z}}{\sin2\theta }\right\rangle ^{2}}$$ (32)

r can be computed from $k=\frac{3Gm}{4r^{3}}(J_{z}-\frac{J_{x}+J_{y}}{2})$. In order to get an accurate value of r, one has to be very careful in selecting the period over which the average is taken. During the simulation we found that the angle between the Mars equator and the orbital planes of Phobos and Deimos have periods of 6050 and 6050/4 years, respectively (see Figs. 1 and 2). Averages must be taken over the full periods. The resulting averages of r for Phobos and Deimos are $7.222214 \times 10^{-5}$ and $1.804138 \times 10^{-4}$AU, respectively. Note that this results in much greater average torques acting on Mars than do averages taken over shorter periods. A 10 year mean treatment would instead result in $r=1.7251 \times 10^{-4}$ and $4.485 \times 10^{-4}$AU for Phobos and Deimos, respectively. For long simulations, it is important to use the first set of values for r.