The purpose of this post is to explain the derivation of the C code used to calculate the gravity perturbation in the 'EGM96' Earth Gravity Model. As such, it is mathematical 'backfill' and a 'for the record' contribution to this thread.
Part of this thread's earlier dialogue included:
Even if I don’t need that accuracy, is there any way to write the code to include at least the sectorial coefficients (degree = order)?
Yes, give me a few days to get the recursion scheme right and I can upgrade the code to a full gravity model
In this exchange, 'cristiapi' asked for an algorithm to calculate a more complete gravity model of the Earth that went beyond the relatively straightforward one used by Orbiter involving only the assorted 'J_n' contributions to the perturbed gravity field. I said yes, I thought that could be done and produced some C code which (after a little debugging) did as requested.
That code was based on the EGM96 model (i.e., a year 1996 version) which was the latest full dataset for a full, Earth gravity model that I could readily identify. There is a newer version of the same, EGM2008 (i.e., a year 2008 version) but at the time of writing this note I was unable to identify the model coefficients for that model. So, the ensuing C code that I developed was based on the earlier EGM96 model instead. The table of coefficients for the model can be found here:
Data file for the EGM96 gravity model
This first few lines of this data file look like:
Code:
2 0 -0.484165371736E-03 0.000000000000E+00 0.35610635E-10 0.00000000E+00
2 1 -0.186987635955E-09 0.119528012031E-08 0.10000000E-29 0.10000000E-29
2 2 0.243914352398E-05 -0.140016683654E-05 0.53739154E-10 0.54353269E-10
3 0 0.957254173792E-06 0.000000000000E+00 0.18094237E-10 0.00000000E+00
3 1 0.202998882184E-05 0.248513158716E-06 0.13965165E-09 0.13645882E-09
3 2 0.904627768605E-06 -0.619025944205E-06 0.10962329E-09 0.11182866E-09
3 3 0.721072657057E-06 0.141435626958E-05 0.95156281E-10 0.93285090E-10
In total, the data file contains:
"The file contains fully-normalized, unitless spherical harmonic coefficients and their standard deviations, for a gravitational model complete from degree 2 order 0, to degree and order 360."
In total, there 65338 records records in the file. In each row, there are two values of interest (in the third and fourth columns), so the total number of model parameters is around 130,000. Compare this with the gravity model for the Earth used by Orbiter which uses just four perturbation terms ([imath]J_2[/imath], [imath]J_3[/imath], [imath]J_4[/imath] and [imath]J_5[/imath]).
The underlying perturbation series used in the EGM96 model is:
[math]\Delta U_(x,y,z) = -\sum_{n=2}^{360}\frac{\mu\,R^n}{r^{n+1}}\,\sum_{m=0}^{n}\left(C_{n,m}\,\cos(m\,\phi) + S_{n,m}\,\sin(m\,\phi) \right)\,\tilde{P}_{n,m}(\xi)[/math]
where:
[imath](x,\,y,\,z)[/imath] is the spatial point (in the ECEF reference frame) at which one wishes to calculate the gravitational perturbation such that the z-axis points along the Earth's instantaneous rotation axis; and the x-axis points towards the prime meridian in that frame.
[imath]r = \sqrt{x^2+y^2+z^2}[/imath]
[imath]\xi = z / r[/imath]
[imath]\phi = \tan^{-1}(y/x)[/imath] is the longitudinal angle in the equatorial reference frame
[imath]\mu[/imath] is the gravitational parameter (the product 'GM') used to normalise the perturbation coefficients
[imath]R[/imath] is the reference radius used to normalise the perturbation coefficients
[imath]C_{n,m}[/imath] and [imath]S_{n,m}[/imath] are the coefficients picked up from the EGM96 model data file such that [imath]n[/imath] is the value in the first column of the data file; [imath]m[/imath] is the value in the second column; [imath]C_{n,m}[/imath] is the value from the third column; and [imath]S_{n,m}[/imath] is the value from the fourth column.
[imath]\tilde{P}_{n,m}(\xi)[/imath] are the 'fully-normalised associated Legendre functions'. Details of these functions can be found at
fully-normalised associated Legendre functions
As always, the perturbation force due to all of these non-spherical gravity terms is given by:
[imath](g_x,\, g_y,\, g_z) = (-\frac{\partial\,U}{\partial x}, -\frac{\partial\,U}{\partial y}, -\frac{\partial\,U}{\partial z})[/imath]
So, to calculate the perturbation contributions, we have to take the partial derivatives of:
[imath]r^{-n-1}\,\cos(m\,\phi)\,\tilde{P}_{n,m}(\xi)[/imath]
and
[imath]r^{-n-1}\,\sin(m\,\phi)\,\tilde{P}_{n,m}(\xi)[/imath]
with respect to [imath]x[/imath], [imath]y[/imath] and [imath]z[/imath].
Bearing in mind that [imath]r=\sqrt{x^2 + y^2 + z^2}[/imath], [imath]\xi=z/r[/imath]; and [imath]\phi = \tan^{-1}(y/x)[/imath], we find that the gravitational perturbation contributions from the potential term:
[math]\Delta UC_{n,m} = -C_{n,m}\,\frac{\mu\,R^n}{r^{n+1}}\,\cos(m\,\phi)\,\tilde{P}_{n,m}(\xi)[/math]
are:
[math]\begin{aligned}
-\frac{\partial}{\partial\,x}\Delta UC_{n,m} &= -C_{n,m}\frac{\mu}{r^4}(\frac{R}{r})^n\,\left(\left((1+n) \,x \,r \,\cos (m \,\phi )-m \,\left(\frac{r}{\rho }\right)^2 \,y \,r \,\sin (m \,\phi )\right) \tilde{P}[I]{n,m}(\xi ) + x \,z \,\cos (m \,\phi ) \, \tilde{P}[/I]{n,m}'(\xi)\right) \\
-\frac{\partial}{\partial\,y}\Delta UC_{n,m} &= -C_{n,m}\frac{\mu}{r^4}(\frac{R}{r})^n\,\left(\left((1+n) \,y \,r \,\cos (m \,\phi )+m \,\left(\frac{r}{\rho }\right)^2 \,x \,r \,\sin (m \,\phi )\right) \tilde{P}[I]{n,m}(\xi ) + y \,z \,\cos (m \,\phi ) \,\tilde{P}[/I]{n,m}'(\xi )\right) \\
-\frac{\partial}{\partial\,z}\Delta UC_{n,m} &= -C_{n,m}\frac{\mu}{r^4}(\frac{R}{r})^n\,\left(\left((1+n)\,z\,r\,\cos (m \,\phi ) \,\tilde{P}[I]{n,\,m}(\xi ) - \left(x^2+y^2\right)\, \cos (m\, \phi ) \,\tilde{P}[/I]{n,m}'(\xi )\right)\right)
\end{aligned}[/math]
and the gravitational perturbation contributions from the potential term:
[math]\Delta US_{n,m} = -S_{n,m}\,\frac{\mu\,R^n}{r^{n+1}}\,\sin(m\,\phi)\,\tilde{P}_{n,m}(\xi)[/math]
are
[math]\begin{aligned}
-\frac{\partial}{\partial\,x}\Delta US_{n,m} &= -S_{n,m}\frac{\mu}{r^4}(\frac{R}{r})^n\left(\left((1+n) \,x \,r \,\sin (m \,\phi )+m \,\left(\frac{r}{\rho }\right)^2 \,y \,r \,\cos (m \,\phi )\right) \tilde{P}[I]{n,m}(\xi )+x \,z \,\sin (m \,\phi ) \,\tilde{P}[/I]{n,m}'(\xi )\right) \\
-\frac{\partial}{\partial\,y}\Delta US_{n,m} &= -S_{n,m}\frac{\mu}{r^4}(\frac{R}{r})^n\left(\left((1+n) \,y \,r \,\sin (m \,\phi )-m \,\left(\frac{r}{\rho }\right)^2 \,x \,r \,\cos (m \,\phi )\right) \tilde{P}[I]{n,m}(\xi )+y \,z \,\sin (m \,\phi ) \,\tilde{P}[/I]{n,m}'(\xi )\right) \\
-\frac{\partial}{\partial\,z}\Delta US_{n,m} &= -S_{n,m}\frac{\mu}{r^4}(\frac{R}{r})^n\left(\left((1+n)\,z\,r\,\sin (m \,\phi ) \,\tilde{P}[I]{n,m}(\xi )-\left(x^2+y^2\right)\, \sin (m\, \phi ) \,\tilde{P}[/I]{n,m}'(\xi )\right)\right)
\end{aligned}[/math]
where [imath]\rho^2 = x^2 + y^2[/imath].
Notwithstanding that these expressions are a bit of a mess, we can readily calculate most of the terms and factors on the right-hand side of each. However, as before fro the simpler ‘[imath]J_n[/imath] only' case, we have to calculate [imath]\tilde{P}[I]{n,\,m}(\xi)[/imath] and [imath]\tilde{P}'[/I]{n,\,m}(\xi)[/imath] for all [imath]2\le n \le N[/imath] and for all [imath]0\le m \le N[/imath] for some cutoff value of [imath]n[/imath], [imath]N[/imath]. And, as before, we rely upon recursion to allow us to do this. Since we are dealing with the 'fully-normalised associated Legendre functions', the recursion relation is different. But following
fully-normalised associated Legendre functions, we can write the recursion relations as:
[math]\tilde{P}[I]{n,\,m}(\xi) = a_{n,m}\,\xi\,\tilde{P}[I]{n-1,m}(\xi) - b[/I]{n,m}\,\tilde{P}_{n-2,m}(\xi)[/math]
for any [imath]n > m[/imath]; and
[math]\tilde{P}[I]{n,n}(\xi) = \sqrt{\frac{2\,n+1}{2\,n}}\,u\,\tilde{P}[/I]{n-1,n-1}(\xi)[/math]
where
[math]u=\sqrt{1-\xi^2}[/math]
and where
[math]a_{n,\,m} = \sqrt{\frac{4\,n^2 - 1}{n^2 - m^2}}[/math]
for [imath]n>m[/imath]; and
[math]b_{n,m} = a_{n,m} / a_{n-1,m}[/math]
for [imath]n>m+1[/imath]; and
[math]b_{n,n-1} = 0[/math]
This recursion sequence is started by setting intial values as follows:
[math]\begin{aligned}
\tilde{P}_{0,0} &= 1 \\
\tilde{P}_{1,0} &= \sqrt{3}\,\xi \\
\tilde{P}_{1,1} &= \sqrt{3}\,u
\end{aligned}[/math]
And given the recursion relations from just these three starting values, all of the values of the associated Legendre functions follow. In addition, by differentiating these recursion relations with respect to [imath]\xi[/imath], we can establish the recursion relations for [imath]\tilde{P}'_{n,m}(\xi)[/imath] as well.
Building the code
We are now in a position to build the C routine to calculate the gravitational perturbations for the EGM9 gravity model. Since we are likely to make repeated calls to the gravity perturbation function, the first thing we want to do is to read the EGM96 perturbation parameter data file and store the results somewhere for easy access. Stripping out boilerplate code, we need to something like:
Code:
file = fopen("-- path to EGM96 data file --", "r");
for (int n = 2; n <= nmax; n++) {
for (int m = 0; m <= n; m++) {
fscanf(file , "%i %i %lf %lf %lf %lf \n" , &i, &j, &C[n][m], &S[n][m], &dum1, &dum2);
}
}
fclose(file);
This code reads the data file and stores the perturbation coefficients in two arrays, C and S, up to some maximum perturbation order 'max'.
Next, and rather than calculate the perturbation coefficients each time we calculate a gravity perturbation, we do a one-off calculation of those coefficients storing these values in two arrays, 'a' and 'b':
Code:
for (int n = 1; n <= nmax; n++) {
for (int m = 0; m < n; m++) {
a[n][m] = sqrt((double)(4*n*n - 1) / (double)(n*n - m*m));
}
}
for (int n = 2; n <= nmax; n++) {
for (int m = 0; m < n-1; m++) {
b[n][m] = a[n][m] / a[n-1][m];
}
b[n][n-1] = 0.0;
}
Now, when we want to calculate a gravity perturbation, we need to apply the recursion relations set out above to calculate [MATH]\tilde{P}[I]{n,m}(\xi)[/MATH] and [MATH]\tilde{P}'[/I]{n,m}(\xi)[/MATH] for all orders [MATH]2\le n \le \text{nmax}[/MATH and [MATH]0\le m\le n[/MATH]. The code for doing this is:
Code:
r = sqrt(x*x + y*y + z*z);
t = z / r;
u = sqrt(1 - t*t);
P[0][0] = 1.0;
P[1][0] = sqrt(3.0) * t;
P[1][1] = sqrt(3.0) * u;
DP[0][0] = 0.0;
DP[1][0] = sqrt(3.0);
DP[1][1] =-sqrt(3.0) * t / u;
for (int n = 2; n <= nmax; n++) {
for (int m = 0; m < n; m++) {
P[n][m] = a[n][m] * t * P[n-1][m] - b[n][m] * P[n-2][m];
DP[n][m] = a[n][m] * P[n-1][m] + a[n][m] * t * DP[n-1][m] - b[n][m] * DP[n-2][m];
}
P[n][n] = sqrt((double)(2*n + 1) / 2 / n) * u * P[n-1][n-1];
DP[n][n] = sqrt((double)(2*n + 1) / 2 / n) * (u * DP[n-1][n-1] - t * P[n-1][n-1] / u);
}
The values and derivatives of the fully-normalised associated Legendre functions are stored in arrays 'P' and 'DP' respectively.
In addition to this recursion, we can also calculate the [MATH]\cos(m\,\phi)[/MATH] and [MATH]\sin(m\,\phi)[/MATH] recursively. The code to this is:
Code:
rho = sqrt(x*x + y*y);
cn[0] = 1.0;
sn[0] = 0.0;
for (int m = 1; m <= nmax; m++) {
cn[m] = (x * cn[m-1] - y * sn[m-1]) / rho;
sn[m] = (y * cn[m-1] + x * sn[m-1]) / rho;
}
with the results of the recursion stored in the arrays 'cn' and 'sn'.
Now, finally, we can calculate the gravitational perturbation by calculating all of the sine and cosine terms to all relevant orders. This we do with the code snippet:
Code:
double mu = 3986004.415e+8; //[M]^3[S]^-2
double R = 6378136.3; //[M]
k = mu / r / r / r / r;
phi = R / r;
q = phi;
for (int n = 2; n <= nmax; n++) {
q *= phi;
for (int m = 0; m <=n; m++) {
FC[n][m][0] = -C[n][m]*k*q*(r*((n+1)*x*cn[m]-m*pow(r/rho, 2)*y*sn[m])*P[n][m]+x*z*cn[m]*DP[n][m]);
FC[n][m][1] = -C[n][m]*k*q*(r*((n+1)*y*cn[m]+m*pow(r/rho, 2)*x*sn[m])*P[n][m]+y*z*cn[m]*DP[n][m]);
FC[n][m][2] = -C[n][m]*k*q*(r* (n+1)*z*cn[m]*P[n][m]-(x*x + y*y)*cn[m]*DP[n][m]);
FS[n][m][0] = -S[n][m]*k*q*(r*((n+1)*x*sn[m]+m*pow(r/rho, 2)*y*cn[m])*P[n][m]+x*z*sn[m]*DP[n][m]);
FS[n][m][1] = -S[n][m]*k*q*(r*((n+1)*y*sn[m]-m*pow(r/rho, 2)*x*cn[m])*P[n][m]+y*z*sn[m]*DP[n][m]);
FS[n][m][2] = -S[n][m]*k*q*(r* (n+1)*z*sn[m]*P[n][m]-(x*x + y*y)*sn[m]*DP[n][m]);
}
}
The results are stored in the arrays, 'FC' and 'FS'. Here 'mu' and 'R' are defined quantities in the EGM96 model and, for this gravity model we must use the listed values in our calculations (see
EGM96 model description file)
And, then, we can calculate sum over all of these terms to give us the total gravitational perturbation:
Code:
f[0] = 0;
f[1] = 0;
f[2] = 0;
for (int n = 2; n <= nmax; n++) {
for (int m = 0; m <= n; m++) {
f[0] += FC[n][m][0]; f[0] += FS[n][m][0];
f[1] += FC[n][m][1]; f[1] += FS[n][m][1];
f[2] += FC[n][m][2]; f[2] += FS[n][m][2];
}
}
And that is basically all there is to it.
