- 2.1 Definition of units
- 2.2 Initial conditions and input data
- 2.3 Treatment of the N-body problem
- 2.4 The instantaneous behavior of the spin axis of a rigid body - rotation of the Earth
- 2.5 External torques acting on the planet
- 2.6 The summer solstice for the northern hemisphere
- 2.7 Rotation of Mars
- 2.8 The length of the summer half year
- 2.9 Spin axis precession and orbital precession
- 2.10 Torques due to Phobos and Deimos

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
=
m). Moreover the gravitational constant is chosen
as *G*=1. In this system of units, 1 Solar mass is =
.

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.

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) | |

Moon/Earth | 0.012300034 |

Phobos/Mars | |

Deimos/Mars |

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*^{2} =0.3307007,
*R*=6378.136 km,
(JPL). This gives that
and
.
For Mars, the moments of inertia are:
corresponding to
*J*_{z}/*mR*^{2} =0.3662,
(JPL) and *R*=3397.2 km (Bouquillon & Souchay 1999);
and
from
*J*_{z}/*J*_{x} =1.005741and
*J*_{z}/*J*_{y} =1.005044 (Bouquillon & Souchay 1999).

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

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
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. 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 ,
we check that for two distinct simulations the numerical truncation error
is negligible:
.
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
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.

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.
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
is always
diagonal. Furthermore, the body is considered to be a rigid (not necessarily
homogeneous) oblate ellipsoid with principal moments of inertia
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 (
)
acting
on the body w.r.t. the center of mass *G* is related to the angular momentum
(
)
according to:

The similarity transformation between the

where remains diagonal at all times and is the Eulerian rotation matrix in the

If we define
,
the matrix elements of
then depend on time according
to:

We thus obtain the following first order nonhomogeneous differential equation:

By replacing d*t* with
,
we end up with
the numerical expression:

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: , where radians/d. The second factor is the ratio between a Julian day and a sidereal day (JPL). The obliquity at MJD 51600.5 is degrees (JPL). will certainly at a later stage precess away from this direction but never away from the

This means that given in

The exact expression for external torques from particles *j* acting on
an ellipsoid at time *t* may in the body fixed system
be written:

where is pointing at an arbitrary mass element d

where is the moment of inertia about the direction of .

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 is computed in

Suppose is the radius vector from the center of the Sun pointing at the Earth-Moon system's center of mass and is its velocity. The normal vector to the ecliptic plane is then given by . Let be the spin angular velocity of Earth. If is fulfilled we have either summer or winter solstice. It may then be checked if the angle in which case we have summer solstice. However, for simulations backward in time the condition for summer solstice instead becomes since is then replaced by . An analogous calculation can also be performed to find the northern summer solstices in the case 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
radians/d. The obliquity (the angle between the
Mars equator and its instantaneous orbital plane) of Mars is given by
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
and
are known,
the equations

completely determine the components of . Here, ( ) is the normal vector to the orbital plane of Mars and is the heliocentric position of Mars at its summer solstice. The obtained results in radians/d are:

The initial conditions for the Mars fixed system ( ) is . The Euler angle is no longer directly related to the obliquity because the Mars orbital plane has an inclination of about relative to the ecliptic. This gives that the angle between (

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:
.
The general relation between *E* and the true anomaly *v* is:

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

The angle from perihelion to the vernal equinox is given by and for the autumnal equinox, . The change in time d

where

where and are the eccentric anomalies at autumnal and vernal equinoxes respectively. According to Eq. (17) we have:

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

The precession angle of the spin axis (
)
is calculated by

where and are normal vectors to the instantaneous and present orbital plane, respectively, and analogously for and .

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
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:

where and (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 , the argument of perihelion and the inclination

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.,
.
The expression for *I* in Eq. (9) may be rewritten as:

where the angles to a moon-position are measured relative to the body fixed system . Let us now use a coordinate system

We treat a mean virtual moon orbit as a simple circular orbit which in polar coordinates becomes: . We thus obtain:

As usual, the moon revolution () and the rotation of Mars () are extremely rapid compared to the change in allowing a time average over an integer number of periods to be performed (the symbol refers to a time average, e.g. for a function

In the same way one obtains that and . The average moment of inertia can now be written:

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

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

The torque is in the direction of , where is the normal vector to the instantaneous Mars orbit and 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

Copyright ESO 2002