Chapter 6, Section 6.10.1

Here we are going to derive explicit expressions for the secular and resonant terms associated with the motion of an asteroid at the 3:1 resonance with Jupiter.

We start by loading the DisturbingFunction package.  In this example it is assumed that the file DisturbingFunction.m is located in the current directory/folder.

<<DisturbingFunction.m;

Show how the function ValidArguments is used.

?ValidArguments

Use ValidArguments to find the secular arguments (associated with terms up to second order in the eccentricities and inclinations), i.e. those independent of both mean longitudes, or, of a form 0λ'+0λ+ ... .  There are 3 such arguments to second order.

argsec = ValidArguments[0,0,2]

Show how the function DisturbingFunctionE is used.

?DisturbingFunctionE

Map DisturbingFunctionE on to each of the secular arguments in order to find the terms (to second order) associated with each.

secular = Map[DisturbingFunctionE[#,2]&,argsec]

Now use ValidArguments to find the terms (up to second order) associated with the arguments of the form 3λ'–λ+... .  There are 6 such arguments.

args31 = ValidArguments[3,-1,2]

Map DisturbingFunctionE on to each of the resonant arguments in order to find the terms (to second order) associated with each.

resonant = Map[DisturbingFunctionE[#,2]&,args31]

Define the value of α to be used in the calculation of the Laplace coefficients (see Eq.(6.206)).

α = N[(1/3)^(2/3)*(1/(1+1/1047.355))^(1/3)]

Show how the function LaplaceCoefficient is used.

?LaplaceCoefficient

Calculate numerical values for the terms associated with the secular arguments ...

secterms = 
secular /. {𝒢->1, m'->1, a'->1,
 b[s_, j_, nd_] :> LaplaceCoefficient[s,j,α,nd]}

... and the resonant arguments.

resterms = 
resonant /. {𝒢->1, m'->1, a'->1,
 b[s_, j_, nd_] :> LaplaceCoefficient[s,j,α,nd]}

The numerical values of secterms and resterms are identical to those given in Table 6.1.  However, we can also make use of the explicit expressions given in Eqs.(6.195)–(6.205) to check these values.  Here we use the defined function LaplaceSeries.

?LaplaceSeries

Now we can find the values of a0, a1, ..., a10.

a0 = LaplaceSeries[1/2,0,α,0,30]/2

a1 = (1/8)*(
2*α*LaplaceSeries[1/2,0,aa,1,30]+
α^2*LaplaceSeries[1/2,0,aa,2,30]) /. aa->α

a2 = (1/2)*(
-α*LaplaceSeries[3/2,1,aa,0,30]) /. aa->α

a3 = 
(1/4)*(
    2*LaplaceSeries[1/2,1,aa,0,30]
-2*aa*LaplaceSeries[1/2,1,aa,1,30]
-aa^2*LaplaceSeries[1/2,1,aa,2,30]
) /. aa->α

a4 = 
aa*LaplaceSeries[3/2,1,aa,0,30] /. aa->α

a5 = 
(1/8)*(
    21*LaplaceSeries[1/2,3,aa,0,30]
+10*aa*LaplaceSeries[1/2,3,aa,1,30]
+aa^2*LaplaceSeries[1/2,3,aa,2,30]
) /. aa->α

a6 = 
(1/4)*(
    -20*LaplaceSeries[1/2,2,aa,0,30]
-10*aa*LaplaceSeries[1/2,2,aa,1,30]
-aa^2*LaplaceSeries[1/2,2,aa,2,30]
) /. aa->α

a7 = 
(1/8)*(
    17*LaplaceSeries[1/2,1,aa,0,30]
+10*aa*LaplaceSeries[1/2,1,aa,1,30]
 +aa^2*LaplaceSeries[1/2,1,aa,2,30])-
 (27/8)*aa /. aa->α

a8 = 
(1/2)*aa*LaplaceSeries[3/2,2,aa,0,30] /. aa->α

a9 = 
-aa*LaplaceSeries[3/2,2,aa,0,30] /. aa->α

a10 = 
(1/2)*aa*LaplaceSeries[3/2,2,aa,0,30] /. aa->α

Again we can compare these values with those given in Table 6.1.

TableForm[
      Transpose[
{Range[0, 10], 
{a0, a1, 
            a2, a3, a4, a5, a6, a7, a8, a9, a10}}]]