From my understanding, it should be possible to follow a similar procedure as what is done for calculating the gradient of gravity potential for the Earth (ignoring the potential corrections due to solid tides of the Moon, not sure how significant these would be, but they seem to exist judging from this article.
For that, I believe we can first look at equation 108 of this technical note, which defines the potential of an attracting body $W_a$ in spherical coordinates:
$$
W_a(r,\lambda,\varphi)=\frac{GM}{R}\sum_{\mathscr{l}=0}^{\mathscr{l}_{max}}\sum_{\mathscr{m}=0}^{\mathscr{l}}\left(\frac{R}{r}\right)^\mathscr{l}P_{\mathscr{lm}}(sin(\varphi))\left((C_{\mathscr{l}m}^Wcos(m\lambda)+S_{\mathscr{l}m}^Wsin(m\lambda)\right)
$$
Where:
- $GM$ is the gravitational parameter of the central body
- $R$ is the reference radius of the central body
- $r$ is the radial distance of the point at which the potential is to be evaluated, in spherical coordinates
- $\lambda$ is its longitude
- $\varphi$ is its latitude
- $P_{\mathscr{lm}}$ are the associated Legendre functions of degree $\mathscr{l}$ and order $\mathscr{m}$
- $C_{\mathscr{l}m}^W$ and $S_{\mathscr{l}m}^W$ are the Stokes' coefficients usually distributed with gravity models, and that I most frequently see denoted as Cnm and Snm, or simply as C and S. These are given in the .tab files you linked in comments for the Moon gravity model.
And equation 120 from the same document gives an expression for the gradient of the gravity potential:
$$
\nabla{W}=W_r\overrightarrow{e_r}+\frac{1}{r\ cos\ \varphi}W_\lambda\overrightarrow{e_\lambda}+\frac{1}{r}W_\varphi\overrightarrow{e_\varphi}
$$
Where:
- $\overrightarrow{e_r}$, $\overrightarrow{e_\lambda}$ and $\overrightarrow{e_\varphi}$ are unit vectors of the spherical coordinate system, for which I find Figure 1.1 here quite useful to understand. Just note that said figure uses polar angle instead of latitude and the symbols for each angle are also different.
- $W_r$, $W_\lambda$ and $W_\varphi$ are the partial derivatives of $W$ with respect to $r$, $\lambda$ and $\varphi$.
So we basically need these partial derivatives to calculate the gradient of potential in spherical coordinates, after which we can transform to any other system of coordinates/frame of reference. Equation 122 of the ICGEM document gives indeed expressions for each of these 3 partial derivatives:
$$
\dfrac{\partial W}{\partial r} = -\frac{GM}{r^2}\sum_{\mathscr{l}=0}^{\mathscr{l}_{max}}\left(\frac{R}{r}\right)^\mathscr{l}(\mathscr{l}+1)\sum_{m=0}^{\mathscr{l}}P_{\mathscr{l}m}(sin\ \varphi)\left(C_{\mathscr{l}m}^Wcos(m\lambda)+S_{\mathscr{l}m}^Wsin(m\lambda)\right)
$$
$$
\dfrac{\partial W}{\partial \lambda} = \frac{GM}{r}\sum_{\mathscr{l}=0}^{\mathscr{l}_{max}}\left(\frac{R}{r}\right)^\mathscr{l}\sum_{m=0}^{\mathscr{l}}mP_{\mathscr{l}m}(sin\ \varphi)\left(S_{\mathscr{l}m}^Wcos(m\lambda)-C_{\mathscr{l}m}^Wsin(m\lambda)\right)
$$
$$
\dfrac{\partial W}{\partial \varphi} = \frac{GM}{r}\sum_{\mathscr{l}=0}^{\mathscr{l}_{max}}\left(\frac{R}{r}\right)^\mathscr{l}\sum_{m=0}^{\mathscr{l}}\dfrac{\partial P_{\mathscr{l}m}(sin\ \varphi)}{\partial \varphi}\left(C_{\mathscr{l}m}^Wcos(m\lambda)+S_{\mathscr{l}m}^Wsin(m\lambda)\right)
$$
One thing I have noticed is that the partial derivative with respect to the latitude should be, I believe, multiplied by the cosine of the latitude, if we follow the chain rule for derivatives. I have asked about this here, but I am starting to believe the omission of the cosine there might be due to the fact that several recursion algorithms for the calculation of the derivative of Legendre functions are, in reality, calculating directly the product of the derivative by the cosine of latitude. I have verified that this is certainly the case for Lear's algorithm from the NASA document for example. So these values would be ready to use in the expression of that partial derivative of the gravity potential.
So the last problem to solve before plugging everything in to be able to calculate the gradient of gravity potential (asides from converting to a Cartesian system of coordinates) is calculating the values of the associated Legendre functions and their derivatives.
Regarding the normalization of Legendre functions and Stokes' coefficients, this is indeed key. The raw legendre functions and Cnm/Snm coefficients are related to their normalized versions by a normalization factor $N_{nm}$, whose expression is given here for example. There seem to be in fact, as you noted, several normalization schemes, with another one being, for example, the so-called Schmidt normalization. In the field of geodesy, however, it seems that (almost) always we use the 4pi, Kaula, geodesy or full normalization (I have seen these 4 terms being applied to the same normalization). The term 4pi normalization seems to be derived from the fact that it makes the squared spherical harmonics yield a value of 4pi when integrated over a sphere.
However, according to Vallado's Fundamentals of Astrodynamics and Applications, section 8.6, the practical reason that such normalizations were adopted in this field is that for spherical harmonics of high degree and order, the values of Stokes' coefficients become very small, while the values of the associated Legendre functions become very large. So the coefficients are divided by the mentioned normalization factor, and the Legendre functions are multiplied by it, and their normalized counterparts are denoted by a straight bar over it (I have also seen a curly bar in some sources, but this seems to denote Schmidt partial normalization). E.g., $\overline{C_{nm}}$
The consequence of this is that the product of Stokes' coefficients and Legendre functions remains unaltered regardless of the applied normalization scheme, but it is key that the same normalization is applied to both. The coefficients in the lunar gravity field seem to be fully normalized, so we need matching fully normalized Legendre functions. These can be obtained with relatively simple recursive relationships, for which several exist. But to provide a reference that I found quite good, this NASA document provides several. I am currently using Lear's algorithm from the 3 in the linked document by the way (unluckily I don't have the code in Python, but happy to provide my implementations in R and C++ if that helps; I decided to write my own implementations because it is quite simple and that way, I am completely certain that the correct normalization is being used). Also, the NASA document provides expressions with which you can verify an implementation of Lear's algorithm.
I have tested this for the case of Earth with the GGM05C model by calculating the modulo of the gradient vector and comparing it with the values reported by this tool for the same point.
Regarding units, Stokes' coefficients (both Cnm and Snm) and the Legendre functions are dimensionless, so the units will be determined by the units of $GM$ and $R$. After plugging in the 3 partial derivatives into the expression for the gradient, the dimensions of the result is in [DISTANCE]/[TIME]^2, i.e., acceleration units, which makes sense for it being the gradient of gravity.
After converting the gradient to a system of Cartesian coordinates, it can be directly used as the acceleration at that point caused by the gravity field of the central point.
Edit:
I have had a go at writing some of the required functions in Python. They can be found here. I am leaving here an overview of what the workflow would be to reconstruct the gradient of gravity in Python:
- Read Cnm and Snm coefficients from model file. It is important to make sure these are "fully normalized" coefficients (although, as far as I have seen, all gravity models of the main celestial bodies use them. But it is important to check in header files).
- Calculate the values of the normalized associated Legendre functions and their derivatives at the target latitude. This can be done with function legendre, which uses Lear's recursion algorithm and receives as input the desired degree (n) and order (m), as well as the target latitude in radians (not colatitude/polar angle!)
- Use the Cnm and Snm coefficients together with the Legendre functions and their derivatives to calculate the gradient of gravity at the desired location in spherical coordinates. This can be done with function gravityGradientSphericalCoords, which receives as input, in this order: array of Legendre functions values (Pnm), array of derivatives of Legendre functions values (dPnm), Cnm, Snm, latitude, longitude, r (radial distance of polar coordinates of the target location), the reference radius for the central body whose gravity field is being calculated (R), the gravitational parameter for said central body (GM), the maximum degree to be considered for the calculation of the gradient of potential (n) and the maximum order (m). In reality, the current implementations will use all orders up to degree n.
This will return the partial derivatives of gravity potential with respect to $r$, $\varphi$ (latitude) and $\lambda$ (longitude), in this order.
- These can then be input into equation 120 of the mentioned technical note to obtain the gravity gradient in spherical coordinates, and then further transformed into a Cartesian system of coordinates.
Edit 2
I think I finally manage to reconstruct and apply Moon's gravity model correctly.
There was, however, one more thing to solve that I did not think of before, which is the reference frame for the resulting acceleration by applying the procedure described above.
After digging into the documentation of the GRGM1200 series of models, it seems that all the GRGM models use the Moon Principal Axes (MPA) frame of reference. This seems in fact very convenient, since according to section 2.4 here, the transformation between MPA and Moon-centered ICRF (let's call it SCRF, since it is a selenocentric frame of reference) only requires the Moon libration angles, which, as described in this answer, are provided in the JPL DE ephemerides.
So basically, if we want to obtain an acceleration in SCRF due to the gravity field of the Moon, the process would be:
- Get the coordinates of the target point in SCRF frame. If we have them in GCRF (Geocentric ICRF), we can easily obtain them by calculating Moon position with the JPL DE ephemerides and subtracting.
- Calculate the rotation matrix from MPA to SCRF, and apply the inverse rotation to the position of the target point in SCRF. This yields the position in MPA frame
- Convert the MPA position to spherical coordinates. This gives a set of spherical coordinates that can be directly input in the procedure I described with my previous edit to calculate the gradient of gravity potential. As mentioned before, it will be initially a vector in spherical coordinates.
- Convert the gradient in MPA spherical coordinates to cartesian coordinates in the same MPA frame. This can be done with the equations 8-27 of section 8.6.1 of Vallado's Fundamentals of Astrodynamics and Applications, but for completeness, I am reproducing them here:
$$
a_x = \left(\frac{1}{r}\dfrac{\partial W}{\partial r} - \frac{r_z}{r^2\sqrt{r_x^2+r_y^2}}\dfrac{\partial W}{\partial \varphi}\right)r_x - \left(\frac{1}{r_x^2+r_y^2}\dfrac{\partial W}{\partial \lambda}\right)r_y
$$
$$
a_y = \left(\frac{1}{r}\dfrac{\partial W}{\partial r} - \frac{r_z}{r^2\sqrt{r_x^2+r_y^2}}\dfrac{\partial W}{\partial \varphi}\right)r_y + \left(\frac{1}{r_x^2+r_y^2}\dfrac{\partial W}{\partial \lambda}\right)r_x
$$
$$
a_z = \frac{1}{r}\dfrac{\partial W}{\partial r}r_z + \frac{\sqrt{r_x^2+r_y^2}}{r^2}\dfrac{\partial W}{\partial \varphi}
$$
Where $a_x$, $a_y$ and $a_z$ are the X, Y and Z components of the acceleration in MPA frame, and $r_x$, $r_y$ and $r_z$ are the X, Y and Z components of the target location in MPA frame as well.
- From here, applying the previously calculated rotation matrix from MPA to SCRF will yield the acceleration in the inertial SCRF frame
I have verified this approach by propagating the position of the GRAIL-A satellite and comparing it with the ephemerides provided by JPL Horizons system. The ephemerides can be obtained here, directly in SCRF frame (ICRF centered in Moon). The results are quite promising! I propagated the orbit for 1 day, outputting results every 30 minutes, since the 4th of April, 2012 at 00:00:00 (UTC). The following plot shows the difference between my calculated positions and Horizons' ephemerides. Note that propagation was performed without taking into consideration any solar radiation pressure, albedo, or infrared radiation pressure. Only gravity acceleration from the Moon, and attraction by the Sun and all planets as point masses, were considered. The black line is for the case of considering only Moon gravity spherical harmonics up to degree 10, and the red line up to degree 150:
Here we have a zoom-in the differences for degree 150:
And it gets better when considering even higher degrees!
That was fun!
Edit 3
For completion, I implemented a short Python function to convert acceleration in spherical coordinates to Cartesian (body fixed). It can be done with function accelSPHtoCART here
