A&A 390, 733-749 (2002)
DOI: 10.1051/0004-6361:20020625

Solar system models with a selected set of frequencies

G. Gómez1,3 - J. J. Masdemont2,3 - J. M. Mondelo2


1 - Departament de Matemàtica Aplicada i Anàlisi, Universitat de Barcelona, Gran Via 585, 08007 Barcelona, Spain
2 - Departament de Matemàtica Aplicada I, Universitat Politècnica de Catalunya, E.T.S.E.I.B., Diagonal 647, 08028 Barcelona, Spain
3 - Institut d'Estudis Espacials de Catalunya, Edifici Nexus, Gran Capità 2-4, 08034 Barcelona, Spain

Received 4 March 2002 / Accepted 19 April 2002

Abstract
The purpose of this paper is to develop a methodology to generate simplified models suitable for the analysis of the motion of a small particle, such as a spacecraft or an asteroid, in the Solar System. The procedure is based on applying refined Fourier analysis methods to the time-dependent functions that appear in the differential equations of the problem. The equations of the models obtained are quasi-periodic perturbations of the Restricted Three Body Problem that depend explicitly on natural frequencies of the Solar System. Some examples of these new models are given and compared with other ones found in the literature. For one of these new models, close to the Earth-Moon system, we have computed the dynamical substitutes of the collinear libration points.

The methodology developed in this paper can also be used for the analytical construction of simplified models of systems governed by differential equations which have a quasi-periodic (in time) external excitation and such that the form of the equations is rather cumbersome.

Key words: celestial mechanics - Solar System - reference systems - ephemerides


1 Introduction

The main goal of this paper is the construction of quasi-periodic analytic models suitable for the study of the motion of a small particle in the Solar System. Without any simplification, the equations of the general problem form a set of 60 first order differential equations difficult to analyze. It is well known that very simple models, such as the Two Body Problem or the Restricted Three Body Problem (RTBP), are suitable for many purposes, since they give a good insight of dynamics in large regions of the phase space of the problem. Some of these models are restricted, which means that the small particle does not have any influence in the motion of the remaining bodies. The models introduced in this paper will be also restricted but not so simple as the ones already mentioned. We will try to keep within them the behavior of the dynamics related to the resonances between natural and excitation frequencies.

Most of the well known restricted problems take as starting point the RTBP. We recall that it models the motion of a particle under the gravitational attraction of two primaries which are assumed to be point masses revolving in circular orbits around their center of mass. The Hamilton function of this system is, in a coordinate system that revolves with the primaries (such a system is called synodic, see Szebehely 1967),

\begin{eqnarray*}H(x,y,z,p_x,p_y,p_z)& = & \frac{1}{2}(p_x^2+p_y^2+p_z^2)+yp_x-x... ...mu)^2+y^2+z^2)^{1/2}}- \frac{\mu}{((x-\mu+1)^2+y^2+z^2)^{1/2}}, \end{eqnarray*}


being $\mu=m_2/(m_1+m_2)$, where m1>m2 are the masses of the primaries. In order to get closer to more realistic situations, this model is modified in different ways. For instance
1.
Assuming that the primaries move in an elliptic orbit instead of a circular one. This is what is known as the elliptic RTBP. Its main difference with the RTBP is that it is non-autonomous and in fact is a time-periodic perturbation of it.
2.
Adding a third primary and setting four body problem restricted models. The simplest one is the Bicircular Restricted Problem. It can be also considered a periodic perturbation of the the RTBP in which one primary has been splitted in two that move around their common center of mass. This model can be suitable to take into account the gravitational effect of the Sun in the Earth-Moon RTBP. In coordinate system revolving with Earth and Moon, the Hamiltonian, H(x,y,z,px,py,pz,t), of this problem is (see Simó et al. 1995)

\begin{eqnarray*}H& = & \frac{1}{2}(p_x^2+p_y^2+p_z^2)+yp_x-xp_y - \frac{1-\mu... ...}}- \frac{m_{\rm S}}{a_{\rm S}^2}(y \sin \theta-x\cos \theta), \end{eqnarray*}


with $\theta=w_{\rm S}t+\theta_0$, where $w_{\rm S}$ is the mean angular velocity of the Sun, $m_{\rm S}$ its mass and $a_{\rm S}$ the distance from the Earth-Moon barycenter to the Sun.
3.
Asking the three primaries to move along a true solution of the three body problem. This provides the so called coherent models. These models have been introduced for the study of the motions around the geometrically defined collinear and triangular equilibrium points of the Earth-Moon system and the Sun-Jupiter system, perturbed by Saturn (see Andreu 1999; Gabern & Jorba 2001). The Hamiltonian of these problems can be written as

\begin{eqnarray*}H& = & \frac{1}{2}\alpha_1(p_x^2+p_y^2+p_z^2)+\alpha_2(xp_x+yp_... ...m_{\rm S}}{((x-\alpha_7)^2+ (y-\alpha_8)^2+z^2)^{1/2}} \right), \end{eqnarray*}


where the $\alpha_i$ are time periodic functions, with the same basic frequency as the Bicircular Problem.

In our approach, instead of taking as starting equations those of the RTBP, we will consider Newton's equation for the motion of an infinitesimal body in the force field created by the bodies of the Solar System

 \begin{displaymath} {\bf R}''=G\sum_{i}m_i\frac{({\bf R}_i-{\bf R})}{\Vert{\bf R}-{\bf R}_i\Vert^3}, \end{displaymath} (1)

where G is the gravitational constant, ${\bf R}$ is the position of the infinitesimal body, ${\bf R}_i$ is the position of the Solar System body i in some suitable inertial system and mi its mass.

Following the ideas of Gómez et al. (2001b), in Sect. 2 we introduce suitable reference systems and units such that, after selecting two bodies of the Solar System as primaries, the above equations are set as a perturbation of the RTBP. These new equations will be a generalization of the ones that we have already introduced for the intermediate models just discussed. In Sect. 3, and for two particular choices of primaries, we will perform the Fourier analysis of all the time periodic functions that appear in the new equations. In this way we are able to introduce a graded set of models with an increasing number of frequencies, that can be considered between the RTBP and the true equations. This is done in Sect. 4. Finally, in the last section we compute the dynamical substitutes of the collinear equilibrium points for one of the intermediate models introduced, close to the Earth-Moon system.

   
2 The equations of motion in the Solar System as a perturbation of the RTBP

Through the full paper, the set of bodies of the Solar System will be denoted by

 \begin{displaymath} {\cal S} = \{ P_1, \dots, P_{11} \} \end{displaymath} (2)

where P1,..., P11 are the nine planets, the Moon and the Sun, respectively. The mass of $P\in{\cal S}$ will be denoted by $m_{\rm P}$.

In an inertial reference system, the Lagrangian related to Newton's equations of motion (Eq. (1)) of an infinitesimal body Q under the gravitational action of the bodies in ${\cal S}$, is

\begin{eqnarray*}L({\bf R},{\bf R}',t^*)&=&\frac12\langle{\bf R}',{\bf R}'\rangle +\sum_{i\in{\cal S}}\frac{Gm_i}{\Vert{\bf R}-{\bf R}_i\Vert}, \end{eqnarray*}


where ${\bf R}=(X,Y,Z)^T$ is the position of Q, the prime denotes the derivative with respect to time, t*, $\langle{\bf R}', {\bf R}'\rangle$ is the dot product between ${\bf R}'$ and ${\bf R}'$, G is the gravitational constant, ${\bf R}_i$ is the position of the body $i\in{\cal S}$ and $\Vert\cdot\Vert$ denotes the Euclidean norm. In practice, it is convenient that the reference frame and units, both in space and time, are consistent with the ephemerides data files used for the determination of ${\bf R}_i$.

Since we are interested in writing the equations of motion for Q as a perturbation of the RTBP equations, we must select two bodies $I,J\in{\cal S}$ with mI>mJ, which will play the role of primaries. In this way, the mass parameter, $\mu$, is defined as $\mu=m_J/(m_I+m_J)$, and so $1-\mu=m_I/(m_I+m_J)$. Next, we must introduce the synodic reference frame. Recall that the origin of this system is set at the barycenter of I, Jand that the positions of the primaries are fixed at $(\mu,0,0)$ and $(\mu-1,0,0)$ (see Szebehely 1967).

Following Gómez et al. (2001b), the transformation from synodical coordinates, ${\bf r}=(x,y,z)^T$, to sidereal ones, ${\bf R}$, is defined by

 \begin{displaymath}{\bf R}={\bf B}+kC{\bf r}, \end{displaymath} (3)

where
-
The translation, ${\bf B}$, is given by

\begin{displaymath}{\bf B}=\frac{m_I{\bf R}_I+m_J{\bf R}_J}{m_I+m_J}, \end{displaymath}

that clearly puts the barycenter of the primaries at the origin,
-
The orthogonal matrix $C=({\bf e}_1,{\bf e}_2,{\bf e}_3)$, sets the primaries on the x-axis and turns the instantaneous plane of motion of the primaries into the xy plane (by requiring that the relative velocity of one primary with respect to the other has its third component equal to zero). The columns of C are

\begin{displaymath}{\bf e}_1=\frac{{\bf R}_{JI}}{\Vert{\bf R}_{JI}\Vert},\quad {... ...s{\bf R}_{JI}'\Vert},\quad {\bf e}_2={\bf e}_3\times{\bf e}_1, \end{displaymath}

being ${\bf R}_{JI}={\bf R}_J-{\bf R}_I$.
-
$k=\Vert{\bf R}_{JI}\Vert$ is a scaling factor which makes the distance between the primaries to be constant and equal to 1.
It is important to remark that this change of variables is non-autonomous, since ${\bf B}$, k and C depend on time through the components of ${\bf R}_I$and ${\bf R}_J$.

The change of coordinates given by Eq. (3) is checked to preserve the Lagrangian form of the equations and the new Lagrangian becomes

\begin{eqnarray*}L({\bf r},{\bf r}',t^*)&=& \frac12\langle{\bf B}',{\bf B}'\rang... ...+\sum_{i\in{\cal S}^*}\frac{Gm_i}{k\Vert{\bf r}-{\bf r}_i\Vert}, \end{eqnarray*}


where ${\bf s}=C{\bf r}$, ${\bf r}_i$ is the position of the body i in dimensionless coordinates and ${\cal S}^*$ represents the set of Solar System bodies without the two primaries I, J. To get the above expression of L, we use that C defines an orthogonal transformation and, hence, it preserves the scalar product and the Euclidean norm.

Finally, we want to use the same time units as those usual for the RTBP, where $2\pi$ time units correspond to one revolution of the primaries. If t* is some dynamical time and n is the mean motion of J with respect to I, then we perform the change of independent variable through

 
t=n(t*-t*0), (4)

where t*0 is a fixed epoch. From now on, t will be called dimensionless time. In Table 1 we give the values of n for the Earth-Moon and the Sun-(Earth+Moon) systems. In this second case, Earth+Moon means the Earth-Moon barycenter and, for this system, the Earth and the Moon are substituted in ${\cal S}$ by a fictitious body of mass mE+mM behaving as their barycenter. The values in Table 1 are averaged values through the 6000 years covered by the JPL ephemerides file DE406 (see Standish 1998), and have been computed from this file. Using Kepler's third law, G(mI+mJ)=n2a3, we can also define the mean semi-major axis of the orbit of one primary around the other; these values are also given in Table 1.


 

 
Table 1: Values for the mean motion and mean semi-major axis used in the Earth-Moon and Sun-(Earth+Moon) cases. The units are $(\mbox{Julian Days})^{-1}$ and km.
  Earth-Moon Sun-(Earth+Moon)
n 0.22997154619514 0.01720209883844
a 384 601.25606767 149 598 058.09228115


If we denote with a dot the derivative with respect to t, then the new Lagrangian can be written as

\begin{eqnarray*}L({\bf r},\dot{\bf r},t) &=& n^2\left(\frac12\langle\dot{\bf B}... ...m_{i\in{\cal S}^*}\frac{Gm_i}{k\Vert{\bf r}-{\bf r}_i\Vert}\cdot \end{eqnarray*}


From this Lagrangian we can remove the term $\langle \dot{\bf B},\dot{\bf B} \rangle$, since it is independent of ${\bf r}$ and $\dot{\bf r}$, and multiply by the scaling factor a/(G(mI+mJ))=1/(n2 a2) without affecting the equations of motion. In this way we get

\begin{eqnarray*}L({\bf r},\dot{\bf r},t) &=& \frac1{a^2}\Bigl(\dot k\langle\dot... ...i\in{\cal S}^*}\frac{\mu_i}{\Vert{\bf r}-{\bf r}_i\Vert}\biggr), \end{eqnarray*}


where $\mu_i=m_i/(m_I+m_J)$.

Since ${\bf e}_1,{\bf e}_2,{\bf e}_3$ form an orthogonal basis, we have that $\langle{\bf e}_i,{\bf e}_j\rangle=\delta_{ij}$, $\langle\dot{\bf e}_i,{\bf e}_j\rangle= -\langle{\bf e}_i,\dot{\bf e}_j\rangle$ and $\langle\dot{\bf e}_i,{\bf e}_i\rangle=0$ for i,j=1,2,3. It can be further shown that $\langle\dot{\bf e}_1,\dot{\bf e}_2\rangle= 0$, $\langle\dot{\bf e}_2,\dot{\bf e}_3\rangle= 0$ and $\langle\dot{\bf e}_1,{\bf e}_3\rangle= 0$. Recalling that ${\bf r}=(x,y,z)^T$, writing ${\bf s}=C{\bf r}={\bf e}_1x+{\bf e}_2y+{\bf e}_3z$ and using the previous relations, we get

\begin{eqnarray*}L({\bf r},\dot{\bf r},t) &=& a_1(\dot x^2+\dot y^2+\dot z^2) +... ...S}^*}\frac{\mu_i}{[(x-x_i)^2+(y-y_i)^2+(z-z_i)^2]^{1/2}}\biggr), \end{eqnarray*}


where the ai are time dependent functions that can be computed in terms of the positions, velocities, accelerations and over-accelerations of the two primaries. Their explicit expressions are

\begin{displaymath}\begin{array}{@{}l@{~}c@{~}l@{\quad}r@{~}c@{~}l@{\quad}r@{~}c... ...}_3\rangle, & &a_{15} =& \displaystyle\frac ak\cdot \end{array}\end{displaymath}

With the above expressions we can write the second-order differential equations of motion as

 \begin{displaymath}\left\{ \begin{array}{rcl} \ddot x &=& b_1+b_4\dot x+b_5\dot ... ...playstyle\frac{\partial\Omega}{\partial z}, \end{array}\right. \end{displaymath} (5)

being
 
$\displaystyle \Omega$ = $\displaystyle \frac{1-\mu}{\sqrt{(x-\mu)^2+y^2+z^2}}+ \frac{\mu}{\sqrt{(x-\mu+1... ...2}} +\sum_{i\in{\cal S}^*} \frac{\mu_{i}}{\sqrt{(x-x_i)^2+(y-y_i)^2+(z-z_i)^2}}$ (6)

where $\mu_i=m_i/(m_I+m_J)$ and the bi time-dependent functions are defined as

\begin{displaymath}\begin{array}{@{}l@{~}c@{~}l@{\quad}r@{~}c@{~}l@{\quad}r@{~}c... ... & &b_{13} =& \displaystyle\frac{a^3}{k^3}\cdot \\ \end{array}\end{displaymath}

We note that setting bi=0 for $i\neq5,7,10,13$, b5=2, b7=b10=b13=1 and skipping the sum over ${\cal S}^*$ in Eq. (6), the Eqs. (5) become the usual RTBP equations with mass parameter $\mu$. Therefore, we can see Eqs. (5) as a perturbation of the RTBP equations. Once the primaries have been fixed, we will get an idea of the order of magnitude of the perturbation, by looking at the first coefficient of the Fourier expansions of the bi functions. The Fourier analysis of this functions will be done in the next sections for two different systems.

   
3 Fourier analysis of the time-dependent functions of the equations of motion

This section is devoted to the results of the Fourier analysis applied to all the time-dependent functions appearing in Eqs. (5), this is: the bj functions and the coordinates, xi, yi, zi, of the bodies of the Solar System included in ${\cal S}^*$. The Fourier analysis follows the methodology developed in Gómez et al. (2001c), which is a refined procedure that allows a very accurate determination of frequencies and amplitudes for analytic quasi-periodic functions. Here we will discuss the selection of the main parameters used in the method as well as the results obtained. Although the analysis can be done for any set of primaries, we have selected two different couples - the Earth and the Moon and the Sun and the Earth-Moon barycenter - because of their relevance in many spacecraft mission analysis simulations.

3.1 Fourier analysis of the b $_\mathsfsl{i}$ functions

Using the algorithm described in Gómez et al. (2001c), we have performed Fourier analysis of the $\{b_i\}_{i=1,\dots,13}$ functions, both for the Earth-Moon case and the Sun-(Earth+Moon) case. This means that for each bifunction we have obtained a set of frequencies and amplitudes that define its quasi-periodic approximation as a trigonometric polynomial, Qbi(t). As for any Fourier procedure, the most relevant parameters to be specified are the size, T, of the time (sampling) interval and the number, N, of equally spaced sampling points chosen in the interval. These parameters define, for instance, the Nyquist critical frequency, $\omega_{\rm c} = N/(2T)$, that fixes the window within we will find all the frequencies (true or aliased) of our time series. So, the first thing that we need is some criteria to choose properly T and N.

Due to our implementation of the Fourier analysis procedure, the parameter N must range over powers of two. For consistency, the length, T, of the time-interval has also been chosen to range over a geometric progression, and the time-interval has always started at January 1st, 2001. The smallest time-interval length, $T_{\rm min}$, has been taken of 95 years (34698.75 Julian days) and the greatest time-interval length, $T_{\rm max}$, has been chosen as the maximum time-interval covered by the JPL DE406 ephemerides after Jan 1st 2001, which is 364938 Julian days (999.15 years). Therefore, we have let T range over the set $\{\delta^n T_{\rm min}\}_{n=0}^{10}$ where $\delta=(T_{\rm max}/T_{\rm min})^{1/10}$. The time units used are revolutions of the secondary (J) around the primary (I) or, equivalently, dimensionless time divided by $2\pi$. The reason for this is that, in this way, the frequency 1.0 corresponds to one revolution of J around I, which has a more intuitive meaning (one lunar month in the Earth-Moon case, one sidereal year in the Sun-(Earth+Moon) case) that will help in the elaboration of the intermediate models of motion. Moreover, in order to evaluate the trigonometric approximations of the bi functions, we only have to multiply the frequencies found by the dimensionless time, without the need of an additional $2\pi$ factor.

The maximum number of samples $N_{\rm max}$ has been chosen to be 220, in order to allow for "comfortable'' runs on machines with 64MB of memory (or, equivalently, bi-processor machines with 128MB). For each value of T, the minimum number of samples has been chosen such that $\frac N{2T}\geq1.5$, in order to make the maximum detectable frequency to be at least 1.5.

Assume that, for certain fixed values of T and N, we have performed Fourier analysis of a given function bi(t) obtaining the trigonometric polynomial Qbi(t). Then, we can easily compute the maximum difference between the analyzed function and its quasi-periodic approximation at the sampling points, that is,

 \begin{displaymath}{d}_{\rm max}=\max_{l=0,\dots,N-1} \bigl\vert b_i(t_l)-Q_{b_i}(t_l)\bigr\vert, \end{displaymath} (7)

where tl= l  (T/N), $l=0,\dots,N-1$ are the sampling epochs. In Figs. 1 (Earth-Moon case) and 2 (Sun-(Earth+Moon) case) we have represented, for all the bi functions, the minimum of ${d}_{\rm max}$ with respect to the different values of N explored, when varying T according to the preceding discussion.


  \begin{figure} \par\includegraphics[width=4.3cm,clip]{H3524F01.eps}\hspace*{5mm}... ...m} \includegraphics[width=4.3cm,clip]{H3524F13.eps}\hspace*{5.8mm}\end{figure} Figure 1: Error results of the Fourier analysis of the bi functions in the Earth-Moon case. For each value of T explored, we have represented the minimum value of ${d}_{\rm max}$with respect to N. The values of T are given in revolutions of the Moon around the Earth.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=4.3cm,clip]{H3524F14.eps}\hspace*{5mm}... ...3mm] \hspace*{5.8mm} \includegraphics{H3524F26.eps}\hspace*{5.8mm}\end{figure} Figure 2: Same as Fig. 1 but for the Sun-(Earth+Moon) case. The values of T are given in revolutions of the Earth+Moon barycenter around the Sun.
Open with DEXTER

To reduce the leakage effect, in all the computations we have multiplied our data by a Hanning function of order two

\begin{eqnarray*}H_T^2(t) &=& \frac{2}{3} \left(1-\cos \left( 2 \pi \frac{t}{T}\right) \right)^2\cdot \end{eqnarray*}


The advantages of the Hanning function with respect to other well-known window functions (see Brigham 1974) are its simplicity and its degree of differentiability. For instance, HTn(t) has degree 2n, whereas a general "triangle window function'' TTn(t) has degree just n. This higher degree of regularity implies a faster decay of the Fourier coefficients (see Gómez et al. 2001c for more details).

In order to control aliasing, two different strategies have been followed. The first one is based on time-domain, and consists in computing the maximum difference between the initial function and its quasi-periodic approximation, over a refinement of the initial grid of data used for the Fourier analysis. This difference will be denoted as $\alpha _1$. If it increases significantly when increasing the number of points of the grid, then aliasing is very likely to occur. For this test, we have used a refinement with 16N equally spaced points in [0,T].

The second anti-aliasing strategy is based on frequency-domain. It consists in computing the number of rightmost consecutive harmonics of the residual Discrete Fourier Transform (DFT) that have modulus less than a fraction of the maximum modulus of the residual DFT. Then, we divide this number by N/2, the total number of harmonics, and this defines the parameter $\alpha _2$. That is, if Ei(t) = bi(t)-Qbi(t) is the error of the trigonometric approximation of bi(t) and cEi,T,N(j), sEi,T,N(j), $j=0,\dots,N/2$, are the cosine and sine coefficients of its DFT, we compute

\begin{eqnarray*}\lefteqn{{p}_{E_i,T,N}(j) = ((c_{E_i,T,N}(j))^2+(s_{E_i,T,N}(j)... ...T,N}(l)\leq{p}_{\max}/25\mbox{ for }l=j,\dots,N/2\} }{N/2}\cdot} \end{eqnarray*}


Then, for instance, a value of 0.2 for $\alpha _2$ means that there are no frequencies greater than $0.8\cdot\omega_{\rm max} = 0.8\cdot (N/2T)$, with amplitude greater than 1/25 times the modulus of the residual DFT, so we do not expect aliasing in the corresponding Fourier analysis. We are assuming here that amplitudes decrease as frequencies increase, which is ensured by the Cauchy estimates of the Fourier coefficients for an analytic quasi-periodic function.

As an example of aliasing and how the two previously-described strategies detect it, we have represented in Fig. 3 the residual DFT of some of the Fourier analysis of the b1 function in the Earth-Moon case. Some numerical values of these analysis are given in Table 2. In the left plot, we see that for N=16384 there are frequencies of high amplitude near $\omega_{\rm max}=4.02903$. As we increase N, the amplitude of the frequencies near $\omega _{\rm max}$ decrease and the values of ${d}_{\rm max}$ as well as the parameter $\alpha _1$ of the first anti-aliasing strategy become closer.

According to this, for the results displayed in Figs. 1 and  % latex2html id marker 2015 $\ref{fig:tscs}$ only those analyses with $\alpha_1<1.2{d}_{\rm max}$ and $\alpha_{2}\geq 0.2$ have been taken into account.

For the generation of simplified models for the Solar System, among all the analysis performed we have selected the best ones in terms of minimum ${d}_{\rm max}$. The corresponding parameters of these "best'' analysis are given in Tables 3 (Earth-Moon) and 4 (Sun-(Earth+Moon)).

3.2 Fourier analysis of the positions of the planets

In order to complete the quasi-periodic approximation of all the time-dependent part in the vector-field given by Eqs. (5), we give in this section the results of the Fourier analysis of the positions of the Solar System bodies in dimensionless coordinates. For each coordinate xi, yi, zi, we have performed Fourier analysis using the same procedure as for the analysis of the bi functions. The minimum value of ${d}_{\rm max}$ with respect to N, for fixed values of T, is plotted in Fig. 4 for the Earth-Moon system. The values of the parameters for the best analysis are given in Table 5. In the electronic version of the paper the reader can find the results for the Sun-(Earth+Moon)system, which can be also provided by the authors.


  \begin{figure} \par\includegraphics[width=4cm,clip]{H3524F27.eps}\hspace*{4mm}\i... ...F28.eps}\hspace*{4mm}\includegraphics[width=4cm,clip]{H3524F29.eps} \end{figure} Figure 3: Modulus of the residual DFT some of the Fourier analysis of the b1 function in the Earth-Moon case. The values of the parameters of these analysis are given in Table 2.
Open with DEXTER


 

 
Table 2: Parameters associated to the Fourier analysis of Fig. 3. From left to right: day0, initial Julian day since Jan 1st, 2001; day$_{\rm f}$, final Julian day (same units); T, time interval in J-revolutions; N, number of points used; $\omega _{\rm max}$, maximum detectable frequency; ${p}_{\rm max}$, maximum modulus of the residual DFT; ${d}_{\rm max}$, maximum difference between b1and its quasi-periodic approximation over the samples; $\alpha _1$, $\alpha _2$, values of the two anti-aliasing parameters.
day0 day$_{\rm f}$ T N $\omega _{\rm max}$ ${p}_{\rm max}$ ${d}_{\rm max}$ $\alpha _1$ $\alpha _2$
366 55 917.4 2033.24 16 384 4.02903 2.66E-05 4.90E-04 2.29E-03 0.0007
366 55 917.4 2033.24 32 768 8.05806 2.66E-05 5.30E-04 5.67E-04 0.1633
366 55 917.4 2033.24 65 536 16.1161 2.66E-05 5.63E-04 5.67E-04 0.5816



 

 
Table 3: Values of the parameters for the best Fourier analyses of the bi functions for the Earth-Moon case.
function T (days) T (J-rev.) N ${p}_{\rm max}$ ${d}_{\rm max}$
b1 55 551.4 2033.24 65 536 2.66E-05 5.63E-04
b2 55 551.4 2033.24 65 536 2.67E-05 5.49E-04
b3 55 551.4 2033.24 32 768 3.30E-06 5.58E-05
b4 55 551.4 2033.24 65 536 2.31E-06 5.01E-05
b5 43 904.0 1606.94 32 768 4.85E-06 9.16E-05
b6 70 288.7 2572.64 32 768 3.92E-08 1.13E-06
b7 55 551.4 2033.24 65 536 3.51E-06 7.81E-05
b8 55 551.4 2033.24 524 288 1.96E-07 5.94E-06
b9 70 288.7 2572.64 65 536 1.97E-08 5.69E-07
b10 55 551.4 2033.24 65 536 3.51E-06 7.83E-05
b11 70 288.7 2572.64 65 536 1.67E-08 5.05E-07
b12 43 904.0 1606.94 32 768 1.58E-06 3.29E-05
b13 55 551.4 2033.24 65 536 3.51E-06 7.99E-05



 

 
Table 4: Values of the parameters for the best Fourier analyses of the bi functions for the Sun-(Earth+Moon) case. Note that, in this case, J-revolutions are sidereal years.
function T (days) T (J-rev) N ${p}_{\rm max}$ ${d}_{\rm max}$
b1 142 382.6 389.815 65 536 4.95E-08 4.40E-07
b2 142 382.6 389.815 65 536 4.95E-08 4.33E-07
b3 112 529.5 308.083 131 072 2.28E-09 2.68E-08
b4 34 698.8 94.998 4096 8.34E-06 6.74E-05
b5 34 698.8 94.998 4096 1.75E-05 1.26E-04
b6 88 935.7 243.488 262 144 1.76E-08 5.71E-07
b7 34 698.8 94.998 4096 1.36E-05 9.17E-05
b8 288 422.1 789.642 524 288 9.65E-08 1.67E-06
b9 88 935.7 243.488 131 072 9.71E-09 3.19E-07
b10 34 698.8 94.998 4096 1.36E-05 9.17E-05
b11 70 288.7 192.436 524 288 2.35E-08 2.38E-06
b12 34 698.8 94.998 4096 3.92E-06 4.06E-05
b13 34 698.8 94.998 4096 1.34E-05 9.47E-05



 

 
Table 5: Values of the parameters for the best Fourier analysis of the coordinates of the Solar System bodies in dimensionless units. Earth-Moon system.
body coord. T (days) T (years) T (J-rev) N ${p}_{\rm max}$ ${d}_{\rm max}$
Mercury x 70 288.7 192.440 2572.64 65 536 1.37E-02 3.41E-01
Mercury y 70 288.7 192.440 2572.64 65 536 1.08E-02 2.89E-01
Mercury z 70 288.7 192.440 2572.64 32 768 3.18E-03 9.99E-02
Venus x 55 551.4 152.091 2033.24 65 536 5.13E-03 1.53E-01
Venus y 55 551.4 152.091 2033.24 65 536 5.60E-03 1.65E-01
Venus z 88 935.7 243.493 3255.14 65 536 1.25E-03 4.10E-02
Mars x 55 551.4 152.091 2033.24 65 536 3.61E-02 8.43E-01
Mars y 180 155.5 493.239 6593.89 131 072 3.21E-02 7.53E-01
Mars z 180 155.5 493.239 6593.89 131 072 3.26E-03 1.38E-01
Jupiter x 55 551.4 152.091 2033.24 32 768 1.40E+00 1.53E+01
Jupiter y 112 529.5 308.089 4118.71 65 536 5.39E-01 1.31E+01
Jupiter z 70 288.7 192.440 2572.64 32 768 1.37E-01 1.31E+00
Saturn x 70 288.7 192.440 2572.64 32 768 6.07E+00 6.19E+01
Saturn y 142 382.6 389.822 5211.36 65 536 2.53E+00 6.46E+01
Saturn z 180 155.5 493.239 6593.89 65 536 3.87E-01 1.04E+01
Uranus x 142 382.6 389.822 5211.36 131 072 2.33E+00 3.75E+01
Uranus y 142 382.6 389.822 5211.36 131 072 2.33E+00 3.76E+01
Uranus z 364 938.0 999.146 13 357.14 131 072 2.42E-01 4.14E+00
Neptune x 288 422.1 789.657 10 556.57 262 144 3.12E+00 4.52E+01
Neptune y 364 938.0 999.146 13 357.14 262 144 2.37E+00 4.51E+01
Neptune z 364 938.0 999.146 13 357.14 131 072 1.80E+00 2.72E+01
Pluto x 364 938.0 999.146 13 357.14 262 144 4.15E+00 1.69E+02
Pluto y 364 938.0 999.146 13 357.14 262 144 2.08E+01 2.93E+02
Pluto z 364 938.0 999.146 13 357.14 131 072 2.42E+00 5.16E+01
Sun x 55 551.4 152.091 2033.24 65 536 4.41E-03 9.73E-02
Sun y 55 551.4 152.091 2033.24 65 536 4.41E-03 9.21E-02
Sun z 34 698.8 95.000 1270.01 16 384 8.49E-04 8.65E-03



  \begin{figure} \par\includegraphics[width=4.3cm,clip]{H3524F30.eps}\hspace*{5mm}... ....eps}\hspace*{5mm} \includegraphics[width=4.3cm,clip]{H3524F44.eps} \end{figure} Figure 4: Error results of the Fourier analysis of the coordinates of the Solar System bodies (in dimensionless coordinates) for the Earth-Moon case. For each value of T explored, we have represented the minimum value of ${d}_{\rm max}$ (in RTBP units) with respect to N. They are continued in Fig. 5.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=4.3cm,clip]{H3524F45.eps}\hspace*{5mm}... ....eps}\hspace*{5mm} \includegraphics[width=4.3cm,clip]{H3524F56.eps} \end{figure} Figure 5: Continuation of Fig. 4.
Open with DEXTER

4 Generation of simplified Solar System models

In this section we will generate several simplified Solar System models using the Fourier approximations computed according to the previous section. The models obtained will be compared with other ones through the computation of residual accelerations.

   
4.1 Adjustment using linear combinations of basic frequencies

In order to give a more physical meaning to the results obtained from the Fourier analysis, we will write the computed frequencies as linear combinations, with integer coefficients, of basic ones. These basic frequencies can be identified as "natural'' frequencies of the planetary and lunar theories. The introduction in the Fourier expansions of the basic frequencies will be the key point for the construction of models of motion with increasing dynamical complexity.

In principle, the basic frequencies will be extracted from the list of frequencies computed in the Fourier analysis and using the procedure explained below. Nevertheless, in some cases it can be convenient to introduce a fixed set of basic frequencies obtained by other means, for instance from an analytical lunar theory, and then write all the computed frequencies as linear combinations of the ones in this fixed set. Both approaches will be considered in what follows.

To set up the algorithms we need two definitions. Assume that $\omega_1,\dots,\omega_n$ is a set of basic frequencies and that a frequency f can be written as $f=k_1 \omega_1+\ldots +k_n \omega_n$ with $k_1,\dots,k_n$integer numbers, then we say that f is a linear combination of $\omega_1,\dots,\omega_n$ of order $k=\vert k_1\vert+\dots+\vert k_n\vert$. We say that fis a linear combination of $\omega_1,\dots,\omega_n$ of order kwithin tolerance $\epsilon>0$ if, for some $k_1,\dots,k_n$ such that $k=\vert k_1\vert+\dots+\vert k_n\vert$, we have

\begin{displaymath}\left\vert f - (k_1 \omega_1+k_2 \omega_2+\ldots +k_n \omega_n) \right\vert < \epsilon. \end{displaymath}

A simple approach for the determination of the basic frequencies is:

1.
choose a maximum order of the linear combinations to be found;
2.
choose a tolerance for the adjustment of frequencies as linear combination of basic ones;
3.
for each frequency, try out all the linear combinations of the current set of basic frequencies up the chosen maximum order;
4.
if any of the linear combinations fulfills the tolerance requirements, add the current frequency to the set of basic ones.
This procedure may add extra basic frequencies (and thus end up with a rationally dependent set) in some cases, for instance, if the current frequency is an integer divisor of one of the basic frequencies. To avoid this, instead of trying to adjust the current frequency as linear combination of the basic ones, we will try to adjust zero as linear combination of the current frequency and the basic ones. If we succeed to do this and the current frequency gets a coefficient different from $\pm$1, it may be necessary to divide some basic frequencies by this coefficient.

These considerations lead to the following procedure for the determination of a basic set of frequencies:

Algorithm 4.1   Given $\{f_1,\dots,f_{N_{\rm f}}\}$ the set of frequencies to be adjusted as linear combination of basic ones to be selected in the set, a tolerance tol for the adjustments and a maximum order maxor for the linear combinations to be found, compute the basis $\{\omega_1,\dots,\omega_{n_{\rm b}}\}$ and the linear combinations $\{(k^i_1,\dots,k^i_{n_{\rm b}})\}_{i=1,\dots,{N_{\rm f}}}$ as 


  $\omega_1\leftarrow f_1$,
$k^1_1\leftarrow1$,
${n_{\rm b}}\leftarrow1$

for $i=2,\dots,{N_{\rm f}}$

  if $0\in{\rm lc}(\{f_i,\omega_1,\dots,\omega_{n_{\rm b}}\},{\rm tol},{\rm maxor})$

     $(k_0,k^i_1,\dots,k^i_{n_{\rm b}})= {\rm adjust}(0,\{f_i,\omega_1,\dots,\omega_{n_{\rm b}}\},{\rm tol},{\rm maxor})$

      for $j=1,\dots,{n_{\rm b}}$

       if $k_0\mathop{\rm divides}\nolimits k^i_j$

    $k^i_j\leftarrow k^i_j/k_0$

  else

     $\omega_j\leftarrow\omega_j/k_0$

     for $l=1,\dots,i-1$

      $k^l_j\leftarrow k^l_j k^{\phantom{l}}_0$

   else 
      ${n_{\rm b}}\leftarrow{n_{\rm b}}+1$

      $(k^i_1,\dots,k^i_{{n_{\rm b}}-1},k^i_{n_{\rm b}})\leftarrow(0,\dots,0,1)$

      for $l=1,\dots,i-1$

         $k^l_{n_{\rm b}}\leftarrow 0.$
In this formulation, we have introduced two functions ${\rm lc}$ and ${\rm adjust}$, defined as follows: Of course, in the actual implementation, the role of these functions is accomplished by the same code.  

In the second case, in which the basic frequencies $\{\omega_1,\dots,\omega_{n_{\rm b}}\}$ are known, we can just take the best linear combination for each frequency. This can be stated as

Algorithm 4.2   Given $\{f_1,\dots,f_{N_{\rm f}}\}$ the set of frequencies to be adjusted as linear combination of the frequency basis $\{\omega_1,\dots,\omega_{n_{\rm b}}\}$, a tolerance tol for the adjustments and a maximum order maxor for the linear combinations to be found, compute the linear combinations $\{(k^i_1,\dots,k^i_{n_{\rm b}})\}_{i=1,\dots,{N_{\rm f}}}$ as 

  for $i=1,\dots,N_{\rm f}$

  $(k^i_1,\dots,k^i_{n_{\rm b}})={\rm adjust}(f_i,\{\omega_l\}_l,{\rm tol},{\rm maxor}).$

   
4.2 Simplified models for the Earth-Moon case

In a rather accurate theory for the lunar motion, as the simplified Brown theory given in Escobal (1968), the fundamental parameters can be expressed in terms of five basic frequencies. In terms of cycles per lunar revolution, their numerical values are

the value of $\omega_5$ is close to the lower frequencies computed in our Fourier expansions and, at the same time, is close to the precision we can expect in the determination of frequencies with the data used (see Gómez et al. 2001c). By these reasons and in order to have also a set of basic frequencies with astronomical meaning, we have adopted for the Earth-Moon models these frequencies as the basic set, instead of the ones provided by Algorithm 2 of the preceding section.

For the Earth-Moon models to be developed in this section, and leaving aside the two primaries - Earth and Moon - the Sun will be the only perturbing body in ${\cal S}^*$. As it will be shown, this provides rather accurate models and, at the same time, avoids the introduction of additional basic frequencies. In this way, in the equations of motion (5) we will only use the Fourier expansions of $b_1,\dots,b_{13}$ and $x_{\rm S},y_{\rm S},z_{\rm S}$ and its general expression for the equations of motion will be

\begin{eqnarray*}\left\{ \begin{array}{rcl} \ddot x &=& b^{(i)}_1+b^{(i)}_4\dot ... ...ystyle\frac{\partial\Omega^{(i)}}{\partial z} \end{array}\right. \end{eqnarray*}


being

\begin{eqnarray*}\Omega^{(i)}&=&\frac{1-\mu_{\rm E,M}}{\sqrt{(x-\mu_{\rm E,M})^2... ...)}_{\rm S})^2+(y-y^{(i)}_{\rm S})^2+(z-z^{(i)}_{\rm S})^2}}\cdot \end{eqnarray*}


The super-index (i) that we have used for the b(i)j, $j=1,\dots,13$, $\Omega^{(i)}$ and $x^{(i)}_{\rm S}$, $y^{(i)}_{\rm S}$, $z^{(i)}_{\rm S}$ functions will be used as a label for the different intermediate models, according to the number of basic frequencies retained in the Fourier expansions.

In Tables 6 and 7 we give partial results relative to the Fourier analyses of $x_{\rm S}$, $y_{\rm S}$, $z_{\rm S}$ and the bifunctions. The full trigonometric expansions can be found in Mondelo (2001). For these tables we have used Algorithm 2 of the preceding section to adjust the frequencies, found by the Fourier method, as linear combinations of the $\{\omega _i\}_{i=1, \dots ,5}$. The values of the parameters used for the adjustment are: ${\it tol}=10^{-6}$ and ${\it maxor}=20$. Only the frequencies associated to the five largest amplitudes are given in the tables.


  

 
Table 6: Fourier analysis results of the $x_{\rm S}$, $y_{\rm S}$ and $z_{\rm S}$ functions. Only five frequencies, associated to the five largest amplitudes, are given. The frequencies have been adjusted as linear combinations, $\sum k_i \omega _i$, of the five basic frequencies. The order of the linear combination, k, and the corresponding error are also displayed.
Func Frequency Amplitude Error k1 k2 k3 k4 k5 k
$x_{\rm S}$ 0.00000000000 -6.27023E-02 0.00000E+00 0 0 0 0 0 0
  0.92519578630 3.86480E+02 -2.11130E-07 0 1 0 0 0 1
  1.91674083000 3.17140E+01 -3.88890E-07 1 1 -1 0 0 3
  0.06634926280 1.32180E+01 3.87440E-08 1 -1 -1 0 0 3
  0.99999608230 1.03360E+01 -3.43580E-07 1 0 0 0 -1 2
$y_{\rm S}$ 0.00000000000 1.60785E-05 0.00000E+00 0 0 0 0 0 0
  0.92519578630 3.87760E+02 -2.11130E-07 0 1 0 0 0 1
  1.91674083000 3.17800E+01 -3.88890E-07 1 1 -1 0 0 3
  0.99999608230 1.03360E+01 -3.43590E-07 1 0 0 0 -1 2
  0.06634926280 8.30700E+00 3.87360E-08 1 -1 -1 0 0 3
$z_{\rm S}$ 0.00000000000 4.24394E-04 0.00000E+00 0 0 0 0 0 0
  0.07882283000 3.40520E+01 -1.09480E-08 1 -1 0 1 0 3
  0.91272219540 9.30940E-01 -1.85070E-07 0 1 -1 -1 0 3
  0.00402231670 9.11850E-01 -9.57650E-08 0 0 0 1 1 2
  1.07036785680 9.31450E-01 -2.05600E-07 2 -1 -1 1 0 5

 
Table 7: First five frequencies of the Fourier analysis of bi functions. The frequencies have been adjusted as linear combinations of $\{\omega _i\}_{i=1, \dots ,5}$. From left to right the columns are: function; frequency; in cycles per lunar revolution, amplitude; error ( ${\rm freq}-k_1\omega _1-\dots -k_5\omega _5$); coefficients of the linear combination that approximates freq, and order of the linear combination ( $\vert k_1\vert+\dots \vert k_5\vert$).
Func Frequency Amplitude Error k1 k2 k3 k4 k5 k
b1 0.00000000000 3.49728E-04 0.00000E+00 0 0 0 0 0 0
  0.92519578630 2.16240E+00 -2.11120E-07 0 1 0 0 0 1
  1.91674083000 1.77450E-01 -3.88880E-07 1 1 -1 0 0 3
  0.85039537680 7.53250E-02 -1.92240E-07 -1 2 0 0 1 4
  0.06634926290 7.39730E-02 3.88600E-08 1 -1 -1 0 0 3
b2 0.00000000000 -6.70000E-09 0.00000E+00 0 0 0 0 0 0
  0.92519578630 2.16960E+00 -2.11120E-07 0 1 0 0 0 1
  1.91674083000 1.77820E-01 -3.88890E-07 1 1 -1 0 0 3
  0.85039537680 7.58320E-02 -1.92220E-07 -1 2 0 0 1 4
  0.06634926260 4.64680E-02 3.85950E-08 1 -1 -1 0 0 3
b3 0.00000000000 -1.41400E-07 0.00000E+00 0 0 0 0 0 0
  0.07882283210 1.90520E-01 -8.87040E-09 1 -1 0 1 0 3
  0.15362345870 6.56920E-03 1.89270E-07 2 -2 0 1 -1 6
  0.91272221270 5.20890E-03 -1.67780E-07 0 1 -1 -1 0 3
  1.07036787670 5.21170E-03 -1.85760E-07 2 -1 -1 1 0 5
b4 0.00000000000 0.00000E+00 0.00000E+00 0 0 0 0 0 0
  0.99154505160 1.07920E-01 -1.69890E-07 1 0 -1 0 0 2
  1.85039157300 2.94710E-02 -4.21940E-07 0 2 0 0 0 2
  0.85884652970 1.68610E-02 -2.43690E-07 -1 2 1 0 0 4
  1.98309009370 8.82140E-03 -3.49210E-07 2 0 -2 0 0 4
b5 0.00000000000 2.00003E+00 0.00000E+00 0 0 0 0 0 0
  0.99154503470 2.17650E-01 -1.86770E-07 1 0 -1 0 0 2
  1.85039156830 4.29420E-02 -4.26650E-07 0 2 0 0 0 2
  0.85884653190 3.81670E-02 -2.41550E-07 -1 2 1 0 0 4
  1.98309007300 1.48070E-02 -3.69960E-07 2 0 -2 0 0 4
b6 0.00000000000 0.00000E+00 0.00000E+00 0 0 0 0 0 0
  0.84637295300 1.44550E-03 -2.03520E-07 -1 2 0 -1 0 4
  1.00401861550 1.44530E-03 -2.22890E-07 1 0 0 1 0 2
  0.01247357960 1.89340E-04 -3.72940E-08 0 0 1 1 0 2
  0.14517208260 1.88980E-04 1.76520E-08 2 -2 -1 1 0 6
b7 0.00000000000 1.00478E+00 0.00000E+00 0 0 0 0 0 0
  0.99154504270 1.65040E-01 -1.78730E-07 1 0 -1 0 0 2
  0.85884652970 3.24780E-02 -2.43700E-07 -1 2 1 0 0 4
  1.85039157280 1.84070E-02 -4.22070E-07 0 2 0 0 0 2
  1.98309009370 1.35090E-02 -3.49200E-07 2 0 -2 0 0 4



 
Table 7: continued.
Func Frequency Amplitude Error k1 k2 k3 k4 k5 k
b8 0.00000000000 -7.00000E-10 0.00000E+00 0 0 0 0 0 0
  1.85039159880 8.24730E-03 -3.96070E-07 0 2 0 0 0 2
  2.84193667480 9.04550E-04 -5.41620E-07 1 2 -1 0 0 4
  0.85884652020 9.17510E-04 -2.53210E-07 -1 2 1 0 0 4
  1.77559103310 5.07100E-04 -5.33340E-07 -1 3 0 0 1 5
b9 0.00000000000 -0.00000E+00 0.00000E+00 0 0 0 0 0 0
  0.84637295300 7.24530E-04 -2.03520E-07 -1 2 0 -1 0 4
  1.00401861550 7.24450E-04 -2.22890E-07 1 0 0 1 0 2
  0.01247357980 4.82170E-05 -3.71330E-08 0 0 1 1 0 2
  0.14517208260 4.80940E-05 1.76380E-08 2 -2 -1 1 0 6
b10 0.00000000000 1.00478E+00 0.00000E+00 0 0 0 0 0 0
  0.99154504270 1.65030E-01 -1.78730E-07 1 0 -1 0 0 2
  0.85884652970 3.24780E-02 -2.43700E-07 -1 2 1 0 0 4
  1.85039157280 1.84070E-02 -4.22070E-07 0 2 0 0 0 2
  1.98309009370 1.35090E-02 -3.49200E-07 2 0 -2 0 0 4
b11 0.00000000000 -0.00000E+00 0.00000E+00 0 0 0 0 0 0
  1.00401861560 7.20820E-04 -2.22850E-07 1 0 0 1 0 2
  0.84637295300 6.06950E-04 -2.03500E-07 -1 2 0 -1 0 4
  0.14517208280 4.66020E-05 1.78760E-08 2 -2 -1 1 0 6
  1.99556364910 3.64300E-05 -4.10800E-07 2 0 -1 1 0 4
b12 0.00000000000 -1.61183E-03 0.00000E+00 0 0 0 0 0 0
  0.99154502640 5.38970E-02 -1.95110E-07 1 0 -1 0 0 2
  1.85039157030 2.69200E-02 -4.24600E-07 0 2 0 0 0 2
  0.85884654110 8.04870E-03 -2.32340E-07 -1 2 1 0 0 4
  1.98309004860 7.32970E-03 -3.94350E-07 2 0 -2 0 0 4
b13 0.00000000000 1.00747E+00 0.00000E+00 0 0 0 0 0 0
  0.99154504270 1.64840E-01 -1.78730E-07 1 0 -1 0 0 2
  0.85884652970 3.15620E-02 -2.43700E-07 -1 2 1 0 0 4
  1.85039157290 2.66550E-02 -4.22010E-07 0 2 0 0 0 2
  1.98309009370 1.34800E-02 -3.49210E-07 2 0 -2 0 0 4



 

 
Table 8: Mean residual accelerations between several models and the real Solar System over selected halo orbits of the RTBP around L2 in the Earth-Moon case. The first column displays the z-amplitude of the halo orbit used as test orbit. The remaining columns show the mean residual acceleration between the corresponding model and the real Solar System over the test orbit.
z-a. RTBP BCP QBCP $\sssm_1$ $\sssm_2$ $\sssm_3$ $\sssm_4$ $\sssm_5$
0.020 0.140126 0.146459 0.138580 0.365299 0.095769 0.010674 0.001374 0.000727
0.025 0.136603 0.142856 0.135174 0.353302 0.093293 0.010388 0.001346 0.000720
0.031 0.132882 0.139025 0.131578 0.340305 0.090590 0.010076 0.001315 0.000711
0.038 0.129087 0.135080 0.127914 0.326550 0.087699 0.009744 0.001282 0.000702
0.048 0.125352 0.131141 0.124312 0.312235 0.084643 0.009393 0.001247 0.000691
0.059 0.121813 0.127324 0.120905 0.297505 0.081429 0.009024 0.001210 0.000678
0.073 0.118614 0.123757 0.117835 0.282462 0.078045 0.008637 0.001171 0.000664
0.091 0.115905 0.120571 0.115249 0.267173 0.074461 0.008229 0.001128 0.000646
0.113 0.113823 0.117895 0.113283 0.251690 0.070634 0.007796 0.001081 0.000625
0.141 0.112471 0.115836 0.112037 0.236056 0.066510 0.007331 0.001030 0.000598
0.175 0.111872 0.114443 0.111533 0.220325 0.062042 0.006831 0.000973 0.000566
0.217 0.111928 0.113663 0.111672 0.204551 0.057196 0.006292 0.000910 0.000526
0.269 0.112400 0.113311 0.112201 0.188782 0.051978 0.005716 0.000840 0.000481
0.300 0.112678 0.113200 0.112492 0.180899 0.049240 0.005417 0.000802 0.000456


From the preceding tables it is clear that the largest bi functions are: b2, whose Fourier expansion has an independent term close to 2, as could be expected from the RTBP equations and b1, b2, which are the "independent terms'' of the planar (z=0) equations of motion (5). We will use these in what follows.

Starting from the frequency basis $\{\omega _i\}_{i=1, \dots ,5}$, we will look for a new basis $\{\nu_i\}_{i=1,\dots,5}$, in terms of which we will generate 5 models, that will be labeled as $\sssm_i$, $i=1,\dots,5$. The index i will indicate the number of frequencies, $\nu_j$, retained in the model. We take $\nu_1=\omega_2$ as the first frequency of the new basis. The reason for that is that it is the main frequency of b1, b2, $x_{\rm S}$ and $y_{\rm S}$, and, in this way, it can be considered the main "planar frequency''. According to this, for the first model, $\sssm_1$, we only keep in the Fourier expansions, the independent terms and those harmonics with the $\nu_1$ frequency. This is coherent with the fact that $\omega_2$ is also the frequency of the BCP and QBCP models.

The second model, $\sssm_2$, will depend on two frequencies: $\nu_1$ and $\nu_2$. We observe that, except for b3, b6, b9, b11 and $z_{\rm S}$, the main frequencies of the remaining functions can be expressed as linear combinations of $\omega_2$ and $\omega_1-\omega_3$. Thus, we will take $\nu_2=\omega_1-\omega_3$. Note that, in this way, bi for i=3,6,9,11 and $z_{\rm S}$ will be poorly approximated in $\sssm_2$, but this will not give a too bad global approximation, because bi for i=3,6,9,11are smaller than the remaining bi, and $z_{\rm S}$ is also smaller than $x_{\rm S},y_{\rm S}$.

The remaining $\nu_i$ have been taken in order to make the sequence of models $\sssm_3$, $\sssm_4$, $\sssm_5$ decreasing in error in the residual accelerations test that will be discussed below. After some trials, we have set

In this way, we have

\begin{displaymath}\left(\begin{array}{c}\nu_1\\ \nu_2\\ \nu_3\\ \nu_4\\ \nu_5\e... ...ega_2\\ \omega_3\\ \omega_4\\ \omega_5 \end{array}\right)\cdot \end{displaymath}

Since the matrix in the above transformation is unimodular, $\{\nu_i\}_{i=1,\dots,5}$ is a valid basic set of frequencies.

Once the different models have been produced, it is desirable to see if they are close or not the the "real'' one, that is: the full equations of motion in which the time periodic functions, bi and xi,yi,zi, are computed using the JPL ephemeris files. We will also compare these intermediate models with more standard simplifications, such as the RTBP, the Bicircular Problem and the Quasi-Bicircular Problem, with the Earth-Moon mass ratio. For these purposes we first select a set of trajectories,

\begin{displaymath}\begin{tabular}{cccc} $\gamma_z$ : & ${\mathbb R}$\space & $\... ...arrow$\space & $({\bf r} (t), \dot{\bf r} (t))$ , \end{tabular}\end{displaymath}

along which the position, ${\bf r} (t)$, and velocity, $\dot{\bf r} (t)$, are known. We have done two kinds of selections. In the first one we have chosen for $\gamma_z$ a family of periodic halo orbits with different z-amplitudes; these orbits are true solutions of the RTBP (see Gómez et al. 2001b for their computation) and cover a large set of solutions with very different sizes. Then, given two models to be compared, with differential equations $\ddot{\bf r}=f({\bf r}, \dot{\bf r}, t)$ and $\ddot{\bf r}=g({\bf r}, \dot{\bf r}, t)$, respectively, and given a trajectory, $\gamma_z$, which does not need to be a true solution of any of the two models, we compute the "mean relative residual acceleration over $\gamma''$ as

\begin{displaymath}\frac1{L}\int_0^T \frac{\Vert f\bigl(\gamma_z(s),t\bigr)-g\bi... ...igl(\gamma_z(s),t\bigr)\Vert}\Vert\dot{\bf r}(s)\Vert{\rm d}s, \end{displaymath} (8)

where t is a fixed epoch (in dimensionless units) and Lis the length of the trajectory $\gamma_z$ (in configuration space).

For the second test the computations are similar except that we have taken instead of $\gamma_z(t)$ a set of points uniformly distributed around a large neighborhood of the equilibrium points. We have also required to their energy (Jacobi constant) to be in a certain interval around the value associated to the equilibrium points. The results obtained are analogous to the ones obtained for the halo orbits, and will not be given here.

It must be noted that the BCP and the QBCP, as stated in the Introduction, assume that for t=0 the Moon and the Sun are in opposition with respect to the Earth. Therefore, we must set the origin of dimensionless time, both in the $\sssm_i$ models and the real Solar System, such that Earth, Moon and Sun are in a configuration close to the one of the BCP and the QBCP for t=0. We have chosen as t=0 the first epoch after Jan. 1st, 2001 at which this happens, this is the Julian day 2451919.3489 (Jan. 9th, 2001).

The results for the residual accelerations are given in Table 8, using as test paths several halo orbits around the collinear equilibrium point L2. As it has already been mentioned, the results with other trajectories, or other equilibrium points, give the same qualitative information. From this table, it becomes clear that the best one-frequency models that we can be used, according to the residual acceleration criteria, are the BCP and the QBCP. But, when we allow two or more frequencies, the models we get fit the JPL one much better. As it has been said, only the Sun has been taken into account in all the intermediate models. By adding additional Solar System bodies, the residual accelerations are of the same order of magnitude than the ones obtained just using the Sun. It is also clear that, from this point of view, there is not a significant improvement between the RTBP and the non-autonomous Bicircular and Quasi Bicircular models.

4.3 Simplified models for the Sun-(Earth+Moon) case

In this case, we proceed as in the Earth-Moon system, except that the basic frequencies will be obtained using the results of the Fourier analysis of Sect. 3 and Algorithm 1 for their determination.

From the numerical data obtained (see the electronic version of the paper), we first observe that the maximum modulus of the highest Fourier coefficient of b1, b2, b3, b6, b8, b9, b11 is 3.521E-05, whereas the minimum modulus of the highest Fourier coefficient of the remaining bi is 1.669E-02. Therefore, in order to detect basic frequencies, we will only take into consideration the b4, b5, b7, b10, b12 and b13 functions. In addition to this simplification, we will not consider any Solar System body in Eq. (6), since, just using the bi, we are already taking the Sun into account.

Applying Algorithm 1 to the b13function, setting ${\it tol}=10^{-5}$, ${\it maxor}=20$, we get the following four basic frequencies:

\begin{displaymath}\nu_1=0.9999926164,\ \nu_2=0.6255242728,\ \nu_3=0.9147445983,\ \nu_4=1.8313395538. \end{displaymath}

These four frequencies allow to adjust the frequencies of the Fourier analysis of the b4, b5, b7, b10 and b12 functions. For that, we have applied the second algorithm of Sect. 4.1 with ${\it tol}={\rm 1E--5}$ and ${\it maxor}=20$. The numerical results can be found in the electronic version of the paper. With these frequencies, we construct the $\sssm_1$, ..., $\sssm_4$ as we did in the Earth-Moon case.


 

 
Table 9: Mean relative residual accelerations between several models and the real Solar System over selected halo orbits of the RTBP around L2 in the Sun-(Earth+Moon) case.
z-a. RTBP $\sssm_1$ $\sssm_4$
0.020000 3.446497E-02 9.901526E-05 8.905454E-04
0.024838 3.411184E-02 9.779360E-05 8.768670E-04
0.030846 3.366579E-02 9.616913E-05 8.589500E-04
0.038308 3.313580E-02 9.416327E-05 8.364166E-04
0.047575 3.254789E-02 9.175134E-05 8.092527E-04
0.059084 3.194355E-02 8.895610E-05 7.776813E-04
0.073376 3.137381E-02 8.582841E-05 7.420444E-04
0.091126 3.089082E-02 8.236183E-05 7.026421E-04
0.113169 3.053770E-02 7.859979E-05 6.597243E-04
0.140545 3.033772E-02 7.450252E-05 6.135638E-04
0.174543 3.028516E-02 7.020714E-05 5.643885E-04
0.216766 3.034115E-02 6.579492E-05 5.127031E-04
0.300000 3.047577E-02 5.898080E-05 4.323859E-04


In Table 9, we compare the models RTBP, $\sssm_1$ and $\sssm_4$ with the real Solar System using the same residual acceleration tests that we used in the Earth-Moon case. We note that the $\sssm_4$ model gives worse results than $\sssm_1$. This is not a contradiction. In the electronic version of the paper we can see that the maximum amplitude of the frequencies of b4, b5, b7, b10 and b12 that are not multiple of $\nu_1$is 6.695E-05. Because of that, adding frequencies does not improve significantly the approximation of the bi functions, and in this way the structure of Eqs. (5) takes over the fact that the bi terms of $\sssm_4$ are closer to the ones of the real Solar System than the corresponding terms of $\sssm_1$.

Therefore, for the Sun-(Earth+Moon) case, we will give $\sssm_1$ as simplified Solar System model. Note that this is a model with very few frequencies that significantly improves the RTBP.

5 Dynamical substitutes of the collinear equilibrium points

As it is well known, the RTBP has five equilibrium points: three of them (L1, L2, L3) are collinear with the primaries and the other two (L4 and L5) form an equilateral triangle with them. Although the intermediate models introduced in the preceding section are close to the RTBP, they are non autonomous, so they do not have any critical point. If we consider the $\sssm_1$ model, since it depends on only one frequency, it can be seen as a periodic perturbation of the RTBP so, under very general non-resonance conditions between the natural modes around the equilibrium points and the perturbing frequency, the libration points can be continued to periodic orbits of the model. These periodic orbits, which have the same period as the perturbation, are the dynamical substitutes of the equilibrium points. In this section we will show these substitutes for the three collinear equilibrium points for $\sssm_1$, in the Earth-Moon system. For the other models, $\sssm_2$,..., $\sssm_5$, as the perturbation is quasi-periodic, the corresponding substitutes will be also quasi-periodic solutions. The methodology for their efficient computation, as well as the results obtained, will appear elsewhere. The dynamical substitutes of the triangular points in the Earth-Moon system, for models close to the ones of this paper, have been studied in Gómez et al. (2001a), Simó et al. (1995) and Jorba (2000) and will not be considered here.


  \begin{figure} \par\includegraphics[width=6.6cm,clip]{H3524F57.eps}\hspace*{7mm}... ...24F58.eps}\\ [5mm] \includegraphics[width=6.3cm,clip]{H3524F59.eps} \end{figure} Figure 6: Dynamical substitutes for the $\sssm_1$model of the three collinear equilibrium points.
Open with DEXTER

The numerical computation of the periodic orbits of $\sssm_1$ that substitute L1 and L3 has no problem and the results obtained are shown in Fig. 6. We can see that L1 is replaced by a very small size periodic orbit and that the substitute of L3 is also almost planar but rather large in the (x,y)-plane. The computation of the substitute of L2, also displayed in Fig. 6, requires more care. Mainly, we need to introduce a continuation parameter between the RTBP and $\sssm_1$, so we consider the 1-parameter family of vector-fields which can be formally written as

\begin{displaymath}(1-\varepsilon) \times ({\rm RTBP}) + \varepsilon \times (\sssm_1). \end{displaymath}

If $\varepsilon=0$ we get the RTBP and when $\varepsilon=1$ we get the desired final model $\sssm_1$.


  \begin{figure} \par\includegraphics[width=6.9cm,clip]{H3524F60.eps}\hspace*{7mm}\includegraphics[width=7.2cm,clip]{H3524F61.eps} \end{figure} Figure 7: a) Evolution with the continuation parameter $\varepsilon $ of the family of periodic orbits starting at the libration point L2. See explanations in the text. b) Planar Lyapunov orbit labeled as A and C in a).
Open with DEXTER

The evolution with the continuation parameter $\varepsilon $ of the family of isoperiodic orbits that start at the libration point L2 is shown in Fig. 7a. In this figure, the point O represents the equilibrium point for $\varepsilon=0$. When the continuation is done starting at O with increasing positive values of $\varepsilon $, we never reach $\varepsilon=1$. The continuation of this branch has been stopped, for the representation, at the point A, where the value $\varepsilon=0$has been reached again. This point corresponds to a planar Lyapunov orbit of the RTBP with frequency equal to the double of the frequency of the model and displayed in Fig. 7b. The other branch, starting for decreasing values of $\varepsilon $, reaches the value $\varepsilon=1$ at the point C. In fact, this curve goes through $\varepsilon=0$ at the point B, which again corresponds to the previously mentioned planar Lyapunov orbit. In the diagram, both points, A and B have y=0 and the different x coordinates correspond to the two orthogonal crossings with y=0 of the planar periodic orbit.

Acknowledgements
The research has been supported by DGICYT grant PB94-0215. Partial support of the Catalan grants CIRIT 1998S0GR-00042 and 2000 SGR-27 is also acknowledged. The numerical explorations have been carried out on HIDRA. We are grateful to CIRIT, DGICYT and the University of Barcelona who are providing funds for this array. J. M. M. wishes also to acknowledge the support of the doctoral research grant 1997FI00136PG from the Generalitat de Catalunya.

References

 

Copyright ESO 2002