For those following this thread, you will undoubtedly have seen some Mathematica, Python and C code for calculating the perturbation terms of a body's gravitational force. And it may well have seemed that the force calculations in these notes seemed to have been manufactured 'from thin air' and bearing little relation to the mathematics of earlier posts in this thread. This post aims to remedy that by setting out the underlying maths used. It will also serves as a 'for the record' document so that I don't entirely forget what I did and how I did it.
I'm going to break this 'for the record' statement into two parts. In the first part, I'll just focus on the algorithm for calculating the perturbations due to the 'J_n' terms - i.e., the kind of gravitational perturbations that Orbiter permits. In the second part, I'll work through the derivation of the force calculations used in the full perturbation model (e.g. the EGM96 model - as per the above posts) and show how that derivation relates to the C code that I posted.
So, let's begin:
The Orbiter gravitational perturbation model - i.e., J_n terms only
At the start of this thread, 'cristiapi' wanted analytical expressions for the J_4 and J_5 perturbation terms of the gravitational perturbation force expressed in cartesian (x-y-z) coordinates. Explicit analytics expressions were given for terms relating to J_2, J_3, J_4 and J_5 - but if you look at those expressions, they quickly become rather cumbersome as the perturbation order is increased. Going beyond n=5 yields analytical expressions that can only be described as hideous. Surely there must be a better way of calculating the perturbation terms rather than relying upon these cumbersome expansion terms?
And, yes, there is. As it turns out, we don't need the analytical expressions at all. For a specific cartesian point (x,y,z), we can calculate the perturbation force at that point using a recursion scheme - rather than relying upon an explicit analytical expression that describes the perturbation force at all spatial points.
And here is how it works...
If our gravitational model contains only the [imath]J_n[/imath] terms, we can write the gravitational potential, [imath]U[/imath] as:
[math]U(x,\,y,\,z) = -\frac{\mu}{r} + \sum_{n=2}^{N}J_n\,\frac{1}{r^{n+1}}\,P_n(\xi)[/math]
where:
[imath]\mu[/imath] is the gravitational parameter ('GM') for the gravitating body in question;
[imath]r = \sqrt(x^2 + y^2 + z^2)[/imath] is the radius of the cartesian point [imath](x,\,y,\,z)[/imath] from the centre of mass of the gravitating body;
[imath]N \ge 2[/imath] is the maximum order of the perturbation expansion;
[imath]J_n[/imath] is the [imath]n[/imath]-th order perturbation expansion coefficient. It is a given tabulated constant. (Note, though, that this not the normalised expansion coefficients [imath]\tilde{J}_n[/imath] quoted in Orbiter's configuration files - but we'll come back to this a bit later.);
[imath]P_n(\xi)[/imath] is the Legendre polynomial of order [imath]n[/imath]; and
[imath]\xi[/imath] is [imath]z / r[/imath] where the [imath]z[/imath]-axis of the cartesian [imath]x[/imath]-[imath]y[/imath]-[imath]z[/imath] reference frame is chose such that it points along the instantaneous rotation axis of the gravitating body.
To calculate the gravitational force vector from this expression, the general procedure is to calculate:
[math](g_x,\,g_y,\,g_z) = (-\frac{\partial\,U}{\partial x}, -\frac{\partial\,U}{\partial y}, -\frac{\partial\,U}{\partial z})[/math]
If we apply this rule to the first term in the expression for [imath]U(x,\,y,\,z)[/imath] - i.e., [imath]-\mu/r[/imath], we get the usual Keplerian gravitational force vector. If we apply the rule to the perturbation terms, then we have to calculate the partial derivatives of [imath]\Delta U_n = J_n\,\frac{1}{r^{n+1}}\,P_n(\xi)[/imath] with respect to [imath]x[/imath], [imath]y[/imath] and [imath]z[/imath]. If we remember that both [imath]r[/imath] and [imath]\xi[/imath] are themselves functions of [imath]x[/imath], [imath]y[/imath] and [imath]z[/imath], we get:
[math]\begin{aligned}
-\frac{\partial\,\Delta U_n}{\partial x} &= J_n\,r^{-n-4} \left((n+1)\,x\,r \,P_n(\xi )+x\,z\,P_n'(\xi )\right) \\
-\frac{\partial\,\Delta U_n}{\partial y} &= J_n\,r^{-n-4} \left((n+1)\,y\,r \,P_n(\xi )+y\,z\,P_n'(\xi )\right) \\
-\frac{\partial\,\Delta U_n}{\partial z} &= J_n\,r^{-n-4} \left((n+1)\,r\,z \,P_n(\xi )-\left(x^2+y^2\right)\,P_n'(\xi )\right)
\end{aligned}[/math]
And since we know [imath]x[/imath], [imath]y[/imath], [imath]z[/imath] and [imath]J_n[/imath] - and, hence, [imath]r[/imath] and [imath]\xi[/imath] - we know immediately evaluate everything in this expression except the Legendre Polynomial, [imath]P_n(\xi)[/imath], and its derivative, [imath]P_n'(\xi)[/imath]. Now, we could use explicit expressions for these functions, but for large [imath]n[/imath], this becomes very unwieldy. A better approach is to use the known recursion relations for the Legendre Polynomials - i.e.,
[math]P_n(\xi) = ((2\,n - 1)\,\xi\,P_{n-1}(\xi) + (1-n)\,P_{n-2}(\xi)) /n[/math]
So, if we know [imath]P_0(\xi)[/imath] and [imath]P_1(\xi)[/imath] we can then calculate [imath]P_2(\xi)[/imath]. And having calculated this, we can then calculate [imath]P_3(\xi)[/imath] and so on. And since we know that:
[math]P_0(\xi) = 1[/math]
and
[math]P_1(\xi) = \xi[/math]
we have all that we need to calculate the Legendre Polynomials to arbitrary order. Equally, we can calculate a similar recursion scheme for the derivatives of the Legendre polynomials. Simply by taking the derivative with respect to [MATH]\xi[/MATH], we get:
[math]P_n'(\xi) = ((2\,n - 1)\,P_{n-1}(\xi) + (2\,n - 1)\,\xi\,P_{n-1}'(\xi) + (1-n)\,P_{n-2}'(\xi)) /n[/math]
and again, we know the starting values for this recursion scheme, namely:
[math]P_0'(\xi) = 0[/math]
and
[math]P_1'(\xi) = 1[/math]
We can easily translate this recursion scheme into C code:
Code:
r = sqrt(x*x + y*y + z*z);
xi = z / r;
P[1] = xi;
P[0] = 1.0;
DP[1] = 1.0;
DP[0] = 0.0;
for (int n = 2; n <= nmax; n++) {
P[n] = (2*n - 1) * xi * P[n-1] + (1-n) * P[n-2]; P[n] /= n;
DP[n] = (2*n - 1) * P[n-1] + (2*n - 1) * xi * DP[n-1] + (1-n) * DP[n-2]; DP[n] /= n;
}
where the array 'P' is used to store the Legendre polynomial values up to and including order 'nmax'; and the array 'DP' is used to store the derivatives of the Legendre polynomials up to and including order 'nmax'. Of course, there are some array declarations and other 'boilerplate' code missing from this - but this contains the core of the recursion scheme.
So, now we are in a position to calculate the total gravitational perturbation force. Before we do so, we need to talk about the normalised perturbation coefficients, [imath]\tilde{J}_n[/imath]. In the above, we have used the unnormalised perturbation coefficients, [imath]J_n[/imath]. This is how Orbiter's config files present the coefficients and, indeed, this is how most people represent these coefficients. What is the connection between the two? Well, it's simply this:
[math]J_n = \tilde{J}_n\,\mu\,R^n[/math]
where [imath]\mu[/imath] is a gravitational parameter ('GM') used in the normalisation - usually just the gravitating body's GM value; and [imath]R[/imath] is a radius used to normalise the coefficients - usually the gravitating body's physical radius. In terms of the normalised perturbation coefficients, then, we get:
[math]\begin{aligned}
-\frac{\partial\,\Delta U_n}{\partial x} &= \tilde{J}_n\,\frac{\mu}{r^4}\,(\frac{R}{r})^n \left((n+1)\,x\,r \,P_n(\xi )+x\,z\,P_n'(\xi )\right) \\
-\frac{\partial\,\Delta U_n}{\partial y} &= \tilde{J}_n\,\frac{\mu}{r^4}\,(\frac{R}{r})^n \left((n+1)\,y\,r \,P_n(\xi )+y\,z\,P_n'(\xi )\right) \\
-\frac{\partial\,\Delta U_n}{\partial z} &= \tilde{J}_n\,\frac{\mu}{r^4}\,(\frac{R}{r})^n \left((n+1)\,r\,z \,P_n(\xi )-\left(x^2+y^2\right)\,P_n'(\xi )\right)
\end{aligned}[/math]
Ans since we have calculated the Legendre polynomials and their derivatives using our recursion scheme, all terms and factors in these expressions are known - so we can immediately write C code to calculate the perturbation contributions up to order max:
Code:
phi = R / r;
k = r * r; k *=k; k = mu / k;
q = phi;
for (int n = 2; n <= nmax; n++) {
q *= phi;
F[n][0] = k * J[n] * q * ((n+1) * P[n] * x * r + DP[n] * x * z) ;
F[n][1] = k * J[n] * q * ((n+1) * P[n] * y * r + DP[n] * y * z) ;
F[n][2] = k * J[n] * q * ((n+1) * P[n] * z * r - DP[n] * (x * x + y * y)) ;
}
Then, having, calculated these, all that we need to do is to sum the contributions:
Code:
f[0] = 0.0;
f[1] = 0.0;
f[2] = 0.0;
for (int n = 0; n <= nmax; n++) {
f[0] += F[n][0];
f[1] += F[n][1];
f[2] += F[n][2];
And aside from boilerplate code, that really is all there is to it. Simple, really.
The more general expansion case used to generate the code for EGM96 (as described earlier in this thread) proceeds in much the same way - although life is made a little more complicated because one must work with the fully-normalised associated Legendre polynomials rather than their simpler Legendre polynomial brethren - but the basic scheme for constructing the perturbation forces of various orders and degrees remains the same. In a later post in this thread, and for completeness, I'll outline the logic behind the construction of the C code.
