Perturbations`DisturbingFunction`

S and S' being two planets or satellites, internal and external respectively, the disturbing functions of the two bodies are expressed by means of the following formulae, where the symbols having a prime refer to the external object:

\(R = \(\(G\ m'\)\/a'\) \((RD - \(r\/r'\^2\) 
              cos \((ψ)\))\)
\)
\(R' = \(\(G\ m\)\/a'\) \((RD - \(r'\/\(\(r\)\(\ \)\)\^2\) 
              cos \((ψ)\))\)
\)
\(RD = 
      a'\/\@\(r\^2 + r'\^2 - 2\ r\ r' cos \((ψ)\)\)
\)
\({m, m'} = mass
\)
\({r, r'} = radius\ vector
\)
\({a, a'} = semi - major\ axis
\)
\(ψ = angle\ \((r\ , r')\)
\)
\(α = a\/a'
\)
G = gravitational constant

RD is the so-called direct part; the other two are denominated indirect parts. For further clarity, here is a list of the symbols adopted for the notation of elements:

a = semi - major axis

e = eccentricity

\(ℑ = inclination
\)
\(Ω = 
      longitude\ of\ the\ ascending\ node
\)
\(ϖ = 
      longitude\ of\ the\ pericentre
\)
\(L = \(M + ϖ = mean\ longitude\)
\)
σ = sin ((ℑ/2))

In expanding RD, Laplace coefficients are used:

((Subscript[b, s])^((j))) (α)

s = (1/2), (3/2), (5/2), …

and in the package, these are represented in this way:

\(b[s, j, 0] = \(b\_s\^\((j)\)\) \((α)\)
\)
b[s, j, 
        p] = ((d ((Subscript[b, s])^((j))) (α))/(d (α^p)))

The expansion of RD has the following form:

RD = (∑_(j = -∞)^(+∞))RDj

j is left unevaluated here, that is, an expansion of RDj is given in which the coefficients of the various cos( ) are expressed as a function of parameter j, of Laplace coefficients and of the elements e, e', σ, σ'. Such coefficients contain a limited number of terms, according to the condition:

\(\(σ\^w\) \(σ'\^x\) \(e\^y\) e'\^z\)=0if w+x+y+z ≤ ord

ord is the order (degree) of the expansion.
I have found that:

\(RDj\ \((ord)\) = ∑\+\(h = H1\)\%H2\(∑\+\(k = K1\)\%K2\(∑\+\(s = S1\)\%S2\(∑\+\(t = T1\)\%T2 coef \((h, k, s, t)\) 
                cos \((\((2\ h - j + s + t)\)\ L + j\ L' - s\ ϖ - 
                    t\ ϖ' - \((h + k)\)\ Ω - \((h - 
                          k)\)\ Ω')\)\)\)\)
\)
\(coef \((h, k, s, 
          t)\) = ∑\+\(z1 = Z1\)\%ZZ1\(∑\+\(z2 = Z2\)\%ZZ2\(∑\+\(z3 = Z3\)\%ZZ3\(∑\+\(z4 = Z4\)\%ZZ4 
                coef1 \((h, k, s, t, z1, z2, z3, 
                  z4)\) σ\^\(\(\(2\ z1\)\(+\)\) | h + k | \)\ σ'\^\(\(\(2\ z2\)\(+\)\) | h - k | \)\ e\^\(\(\(2\ z3\)\(+\)\) | s | \)\ e'\^\(\(\(2\ z4\)\(+\)\) | t | \)\)\)\)
\)
coef1 (h, k, s, t, z1, z2, z3, 
          z4) = (∑_(i = I1)^I2)(∑_(p = P1)^P2)(∑_(u = U1)^U2) coef2 (h, k, s, t, z1, z2, z3, z4, i, p, 
                u) (α^(i + u)) b[i + (1/2), -h + j + p - t, 
                u]

coef2 is a polynomial in j, and has rational coefficients. I elaborated all the summations in order to supply different terms, preventing the emergence of opposite arguments or of equal monomials. coef2 is a combination of numerical functions, that I determined by studying the general formula for the series expansion of the Kaula functions and of Hansen coefficients, making no use of any recursive relations.
The package defines the following functions:

({{"\!\(\*StyleBox[\"\\\"DFExpand[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"ord\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\)", "gives the expansion of RDj up to order \!\(\*StyleBox[\"\\\"ord\\\"\", \"TI\"]\)"}, {"\!\(\*StyleBox[\"\\\"DFArguments[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"ord\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\)", "gives the list of arguments of RDj up to order \!\(\*StyleBox[\"\\\"ord\\\"\", \"TI\"]\)"}, {"\!\(\*StyleBox[\"\\\"DFCoefficient[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"arg\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\",\\\"\", \"MR\"]\) \!\(\*StyleBox[\"\\\"ord\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\) ", "gives the coefficient of argument \!\(\*StyleBox[\"\\\"arg\\\"\", \"TI\"]\) of RD or of RDj up to order \!\(\*StyleBox[\"\\\"ord\\\"\", \"TI\"]\)"}, {"\!\(\*StyleBox[\"\\\"AddIndirectTerms[\\\"\", \"MR\"]\)\!\(\*StyleBox[\"\\\"expr\\\"\", \"TI\"]\)\!\(\*StyleBox[\"\\\"]\\\"\", \"MR\"]\)", "introduces the contribution of indirect terms in the expression \!\(\*StyleBox[\"\\\"expr\\\"\", \"TI\"]\)"}})

Structure of the disturbing function.

This loads the package.

<<Perturbations`DisturbingFunction`

DFExpand

The Appendix B of the textbook "Solar System Dynamics" di C.D.Murray e S.F.Dermott (Cambridge University Press, 1999) [here referred to as SSD] contains the expansion of RD up to the fourth order. This is the result obtained with DFExpand[4]:

DFExpand[4]

The book shows the grouping of monomials within some coefficients. I have not yet been able to identify the general law of this phenomenon or to introduce it in the general architecture of the formula. Therefore, I attributed to the DFExpand function the option Simplification, which may assume either True or False values (the latter being the default choice). With Simplification->True, the program applied makes use of Collect to identify the groupings; I notice that, in the case of more complicated expressions, simplification may take a long time.

I could easily have introduced groupings of Laplace coefficients in pairs, which we may see in some coefficients in SSD, but I decided against it, as this may have interfered with the right functioning of the function AddIndirectTerms.

DFExpand[4,Simplification->True]

The term 4D0.9 found in SSD does not appear in the previous result, but it is replaced by a slightly different one. Bearing in mind that in each term of RD, parameter j must assume the values ranging from -infinity to +infinity, it is deduced that it is possible to use the replacement j->-j, without altering the value of the result. By doing this, it can easily be seen that the term in SSD coincides with the one generated by the function.

DFArguments

DFArguments[ord] supplies an ordered list of all the arguments present in the DFExpand[ord]. Therefore, it proves useful when only some terms of the expansion are to be studied with DFCoefficient.

By looking through the arguments in SSD, we will find j L1-j L+ϖ-ϖ1-Ω+Ω1 in place of j L1-j L-ϖ+ϖ1+Ω-Ω1, but we accounted for this apparent divergence in the previous section.

DFArguments[4]

DFCoefficient

DFCoefficient[arg, ord] is the coefficient of cos(arg) in RDj (or RD).

Let us check the term 4D1.5 on page 543 in SSD:

DFCoefficient[(1 - j) L + j L1 + ϖ - 2 Ω, 
      4]

If we effect the replacement j->-j in the argument, we can see it can be found correctly in the coefficient as well (functions b are normalised so that j is always preceded by the sign +, changing the sign in the second argument when necessary).

DFCoefficient[(1 + j) L - j L1 + ϖ - 2 Ω, 
      4]

The argument may also have all its coefficients in numerical form (without j); for instance, let us consider j=-5 in the previous result:

DFCoefficient[(1 + j) L - j L1 + ϖ - 2 Ω /. 
        j -> -5, 4]

DFExpand[ord] generates certain terms which, as the values for j change, produce infinite other arguments, though not their opposite values (this is because, as I stated at the beginning, the algorithm prevents repetitions). However, cos(arg)=cos(-arg), so using an argument like arg or like its opposite is perfectly equivalent. Since the function DFCoefficient recognizes this, it may be used for any valid forms of the arguments:

DFCoefficient[-(5 L - 2 L1 - ϖ - 
              2 Ω1), 6] === 
      DFCoefficient[5 L - 2 L1 - ϖ - 2 Ω1, 
        6]

AddIndirectTerms

AddIndirectTerms[expr] introduces the contribution of the indirect terms  in one given expression, on condition that this is in an entirely numerical form (without parameter j). The option Perturber may assume the values: External (default choice) and Internal:

df = 
        DFExpand[2, Simplification -> True] /. 
          j -> 2;
AddIndirectTerms[df]

AddIndirectTerms[df, Perturber -> Internal]

co = 
      DFCoefficient[(1 - j) L + j L1 + ϖ - 2 Ω, 
          4] /. j -> 1

AddIndirectTerms[co]

The internal perturber does not supply any indirect terms in this case:

AddIndirectTerms[co, Perturber -> Internal]