<p>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 <a href="https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JE004755" rel="nofollow noreferrer">this article</a>.</p>
<p>For that, I believe we can first look at equation 108 of <a href="http://icgem.gfz-potsdam.de/str-0902-revised.pdf" rel="nofollow noreferrer">this technical note</a>, which defines the potential of an attracting body <span class="math-container">$W_a$</span> in spherical coordinates:</p>
<p><span class="math-container">$$
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)
$$</span></p>
<p>Where:</p>
<ul>
<li><span class="math-container">$GM$</span> is the gravitational parameter of the central body</li>
<li><span class="math-container">$R$</span> is the reference radius of the central body</li>
<li><span class="math-container">$r$</span> is the radial distance of the point at which the potential is to be evaluated, in spherical coordinates</li>
<li><span class="math-container">$\lambda$</span> is its longitude</li>
<li><span class="math-container">$\varphi$</span> is its latitude</li>
<li><span class="math-container">$P_{\mathscr{lm}}$</span> are the associated Legendre functions of degree <span class="math-container">$\mathscr{l}$</span> and order <span class="math-container">$\mathscr{m}$</span></li>
<li><span class="math-container">$C_{\mathscr{l}m}^W$</span> and <span class="math-container">$S_{\mathscr{l}m}^W$</span> 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.</li>
</ul>
<p>And equation 120 from the same document gives an expression for the gradient of the gravity potential:</p>
<p><span class="math-container">$$
\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}
$$</span></p>
<p>Where:</p>
<ul>
<li><span class="math-container">$\overrightarrow{e_r}$</span>, <span class="math-container">$\overrightarrow{e_\lambda}$</span> and <span class="math-container">$\overrightarrow{e_\varphi}$</span> are unit vectors of the spherical coordinate system, for which I find <a href="https://iopscience.iop.org/book/978-1-64327-470-6/chapter/bk978-1-64327-470-6ch1" rel="nofollow noreferrer">Figure 1.1 here</a> quite useful to understand. Just note that said figure uses polar angle instead of latitude and the symbols for each angle are also different.</li>
<li><span class="math-container">$W_r$</span>, <span class="math-container">$W_\lambda$</span> and <span class="math-container">$W_\varphi$</span> are the partial derivatives of <span class="math-container">$W$</span> with respect to <span class="math-container">$r$</span>, <span class="math-container">$\lambda$</span> and <span class="math-container">$\varphi$</span>.</li>
</ul>
<p>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:</p>
<p><span class="math-container">$$
\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)
$$</span>
<span class="math-container">$$
\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)
$$</span>
<span class="math-container">$$
\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)
$$</span></p>
<p>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. <a href="https://earthscience.stackexchange.com/questions/23480/expression-for-the-gradient-of-gravity-potential-in-spherical-coordinates">I have asked about this here</a>, 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.</p>
<p>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.</p>
<p>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 <span class="math-container">$N_{nm}$</span>, whose expression is given <a href="http://mitgcm.org/%7Emlosch/geoidcookbook/node11.html" rel="nofollow noreferrer">here</a> 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.</p>
<p>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., <span class="math-container">$\overline{C_{nm}}$</span></p>
<p>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, <a href="https://ntrs.nasa.gov/api/citations/20110023121/downloads/20110023121.pdf" rel="nofollow noreferrer">this NASA document</a> 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.</p>
<p>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 <a href="http://icgem.gfz-potsdam.de/calcgrid" rel="nofollow noreferrer">by this tool</a> for the same point.</p>
<p>Regarding units, Stokes' coefficients (both Cnm and Snm) and the Legendre functions are dimensionless, so the units will be determined by the units of <span class="math-container">$GM$</span> and <span class="math-container">$R$</span>. 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.</p>
<p>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.</p>
<p><strong>Edit</strong>:</p>
<p>I have had a go at writing some of the required functions in Python. <a href="https://github.com/Rafael-Ayala/astroUtils" rel="nofollow noreferrer">They can be found here</a>. I am leaving here an overview of what the workflow would be to reconstruct the gradient of gravity in Python:</p>
<ol>
<li>Read Cnm and Snm coefficients from model file. It is important to make sure these are &quot;fully normalized&quot; 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).</li>
<li>Calculate the values of the normalized associated Legendre functions and their derivatives at the target latitude. This can be done with <a href="https://github.com/Rafael-Ayala/astroUtils/blob/main/legendre.py" rel="nofollow noreferrer">function legendre</a>, 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!)</li>
<li>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 <a href="https://github.com/Rafael-Ayala/astroUtils/blob/main/gravityGradient.py" rel="nofollow noreferrer">function gravityGradientSphericalCoords</a>, 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 <span class="math-container">$r$</span>, <span class="math-container">$\varphi$</span> (latitude) and <span class="math-container">$\lambda$</span> (longitude), in this order.</li>
<li>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.</li>
</ol>
<p><strong>Edit 2</strong></p>
<p>I think I finally manage to reconstruct and apply Moon's gravity model correctly.</p>
<p>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.</p>
<p>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 <a href="https://iopscience.iop.org/article/10.3847/1538-3881/abd414&amp;lang=en" rel="nofollow noreferrer">according to section 2.4 here</a>, 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 <a href="https://astronomy.stackexchange.com/questions/48264/for-how-many-solar-system-bodies-do-recent-jpl-development-ephemerides-include-l">in this answer</a>, are provided in the JPL DE ephemerides.</p>
<p>So basically, if we want to obtain an acceleration in SCRF due to the gravity field of the Moon, the process would be:</p>
<ul>
<li>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.</li>
<li>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</li>
<li>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.</li>
<li>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:
<span class="math-container">$$
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
$$</span>
<span class="math-container">$$
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
$$</span>
<span class="math-container">$$
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}
$$</span></li>
</ul>
<p>Where <span class="math-container">$a_x$</span>, <span class="math-container">$a_y$</span> and <span class="math-container">$a_z$</span> are the X, Y and Z components of the acceleration in MPA frame, and <span class="math-container">$r_x$</span>, <span class="math-container">$r_y$</span> and <span class="math-container">$r_z$</span> are the X, Y and Z components of the target location in MPA frame as well.</p>
<ul>
<li>From here, applying the previously calculated rotation matrix from MPA to SCRF will yield the acceleration in the inertial SCRF frame</li>
</ul>
<p>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 <a href="https://ssd.jpl.nasa.gov/horizons/app.html#/" rel="nofollow noreferrer">here</a>, 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:
<a href="https://i.sstatic.net/pDhLV.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/pDhLV.png" alt="enter image description here" /></a></p>
<p>Here we have a zoom-in the differences for degree 150:
<a href="https://i.sstatic.net/zLXlY.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zLXlY.png" alt="enter image description here" /></a></p>
<p>And it gets better when considering even higher degrees!</p>
<p>That was <strong>fun</strong>!</p>
<p><strong>Edit 3</strong></p>
<p>For completion, I implemented a short Python function to convert acceleration in spherical coordinates to Cartesian (body fixed). It can be done with <a href="https://github.com/Rafael-Ayala/astroUtils/blob/main/accelSPHtoCART.py" rel="nofollow noreferrer">function accelSPHtoCART here</a></p>

CC BY-SA 4.0
2022-01-31T12:02:14.200Z
Thank you! Where can I get the  and  coefficients properly normalized to use in these expressions?
uhoh
CC BY-SA 4.0
2022-01-31T12:02:30.120Z
Regarding the C and S coefficients properly normalized, the ones in the .tab files that you linked [here](https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/) are already properly normalized for use in these expressions, as indicated by the fact that the corresponding .lbl files state "The spherical harmonic coefficients are fully normalized." Yes, I have tested my solution for the case of Earth. Sorry, I see, I didn't focus it correctly! I will try to edit and add pseudocode more explicitly
Rafa
CC BY-SA 4.0
2022-01-31T12:05:19.900Z
I deleted most of my comment after I read more carefully. Oh, I see, I [linked to that](https://earthscience.stackexchange.com/questions/19566/how-to-make-a-python-function-return-the-moons-gravitational-force-acceleratio/23481#comment44693_19566) back in 2020!
uhoh
CC BY-SA 4.0
2022-01-31T12:05:55.613Z
The key is to use an algorithm that generates Pnm and their derivatives with the same normalization applied. I tried to follow up to the source of scipy's  [lpmv](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.lpmv.html), which is called by scipy's sph_harm, but it leads to FORTRAN code, which I cannot read unluckily. However, Lear's algorithm is relatively straightforward to implement following the instructions on the document by NASA
Rafa
CC BY-SA 4.0
2022-01-31T12:08:27.977Z
Indeed! The model that you linked seems to provide the correct coefficients. I remember also finding at some point ephemerides for some satellite orbiting the Moon. I was planning on using those ephemerides for testing a reconstruction of the gradient of gravity of Moon with this procedure! But unluckily cannot find the ephemerides source again just now... will try to get it again!
Rafa
CC BY-SA 4.0
2022-01-31T16:51:08.097Z
@uhoh I've had a go at coding some of the steps in Python. Hopefully that helps!
Rafa
CC BY-SA 4.0
2022-02-16T11:42:03.867Z
@uhoh I've updated with some additional steps and application results!
Rafa
CC BY-SA 4.0
2022-02-16T22:04:21.990Z
Yes you certainly have! Okay I will dig in today and go through it carefully; I need titer exactly the right amount of coffee to take it in :-)
uhoh
CC BY-SA 4.0
2022-01-31T11:44:03.603Z
2022-02-18T02:52:18.020Z
2022-02-18T02:52:18.020Z
Rafa
<p>The Space SE question <a href="https://space.stackexchange.com/q/43273/12102">How much does it cost to land on a lunar mascon?</a> includes the background information below. I am sure it's not too hard to find a link where spherical harmonic coefficients for some model of the Moon's gravity can be found, but the hard part is to turn them into an actual vector force field in units of m/s<sup>2</sup>.</p>
<p>Likely the coefficients will be for the gravitational potential, and there will be some method to take the gradient to obtain a field, but I have never been able to make sense of this kind of data. For example see:</p>
<ul>
<li><a href="https://space.stackexchange.com/q/27689/12102">Ceres gravity from spherical harmonics from Dawn, how to get the coefficients, definitions and potential?</a></li>
<li><a href="https://space.stackexchange.com/q/27779/12102">How can I verify my reconstructed gravity field of Ceres from spherical harmonics?</a></li>
</ul>
<hr />
<p><a href="https://i.sstatic.net/CEllE.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/CEllE.jpg" alt="These maps of the moon show gravity anomalies measured by NASA's GRAIL mission. (Image credit: NASA/JPL-Caltech/CSM) " /></a></p>
<blockquote>
<p>These maps of the moon show gravity anomalies measured by NASA's GRAIL mission. (Image credit: NASA/JPL-Caltech/CSM)</p>
</blockquote>
<p>From Space.com's <a href="https://www.space.com/21364-moon-gravity-mascons-mystery.html" rel="nofollow noreferrer">Mystery of Moon's Lumpy Gravity Explained</a>, click for full size.</p>
<p>The <a href="https://en.wikipedia.org/wiki/Bouguer_anomaly" rel="nofollow noreferrer">Bouguer anomaly plots</a> show elevated gravity as much as about 600 mGal. Wikipedia's <a href="https://en.wikipedia.org/wiki/Gal_(unit)" rel="nofollow noreferrer">Gal (unit)</a> explains what this means; 0.6 Gal is 0.006 m/s<sup>2</sup> or almost a half of a percent of the Moon's average gravity, and they extend over quite large distances.</p>
<blockquote>
<p>The gal (symbol: Gal), sometimes called galileo after Galileo Galilei, is a unit of acceleration used extensively in the science of gravimetry. The gal is defined as 1 centimeter per second squared (1 cm/s<sup>2</sup>). The milligal (mGal) and microgal (µGal) are respectively one thousandth and one millionth of a gal.</p>
</blockquote>

CC BY-SA 4.0
2020-04-06T11:16:49.110Z
@a_donda orbits are in 3D, I need to calculate everywhere from the surface up to hundreds of kilometers above. In free space gravity obeys the Laplace equation ($\nabla^2 \phi=0)$ so the spherical harmonics work everywhere. I suppose that if I had the potential on one smooth surface I could reconstruct the coefficients for the rest of space as answers to my old question [Calculating the potential on a surface from the potential on another surface](https://physics.stackexchange.com/q/278664/83380) explain,...
uhoh
CC BY-SA 4.0
2020-04-06T11:17:03.067Z
@a_donda ...but they do provide the raw coefficients https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/gggrx_1200a_sha.tab Now all I need is the specific functional form they go into. There are several ways these can be normalized (see my linked questions) so I'll need the correct form. It's the exact equation with the Bessel functions in it that I need. I can calculate orbits no problem, I just need the gradient of the potential field.
uhoh
CC BY-SA 4.0
2020-04-06T11:33:19.850Z
Actually the `.tab` files [here](https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/) contain the coefficients and the `.lbl` files contain the explanations. It's all starting to come back to me now, but I think the question of how the coefficients are normalized still needs to be pinned down, there are several different conventions. See https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-2-edr-v1/grail_0001/document/dpsis.pdf
uhoh
CC BY-SA 4.0
2020-04-06T12:31:07.857Z
@a_donda please take a moment and have a look at the the two questions I've linked to in my question and their answers to see what I mean my normalization. This is a question about math and Earth Scientists do indeed use math!
uhoh
CC BY-SA 4.0
2020-04-06T21:29:33.960Z
With variations as small as half a percent, the gravity vector is almost perfectly radially inward if you put the moon's center at the coordinate origin. Then, if you have the magnitude of the gravity at some particular height datum (as in the data sets mentioned above), you just need to scale it with $1/r^2$ to get the value at some height above the surface.
Wolfgang Bangerth
CC BY-SA 4.0
2020-04-06T22:05:23.257Z
@WolfgangBangerth no that's not sufficient for me. We've already found spherical harmonic coefficients necessary to build the potential field for all points above the Moon's surface and the gradient will provide the acceleration vectors. I just need to find the correct equations and the normalizations for those coefficients.
uhoh
CC BY-SA 4.0
2020-04-08T19:16:21.357Z
@a_donda No. The perturbation from a mass concentration will certainly be important out to a distance of several times it's diameter, and some of these are several hundreds of kilometers in diameter. They were discovered by their substantial perturbations on objects orbiting the Moon and they even resulted in one crashing into the Moon I think. I recommend you take some time to read up on them if you are interested.
uhoh
<open-data><moon><gravity>
How to make a Python function return the Moon's gravitational force (acceleration) to calculate orbits?
uhoh
# How to make a Python function return the Moon's gravitational force (acceleration) to calculate orbits?

(Asked by: username_0 on 2020-04-06)

The Space SE question [How much does it cost to land on a lunar mascon?](https://space.stackexchange.com/q/43273/12102) includes the background information below. I am sure it's not too hard to find a link where spherical harmonic coefficients for some model of the Moon's gravity can be found, but the hard part is to turn them into an actual vector force field in units of m/s2.

Likely the coefficients will be for the gravitational potential, and there will be some method to take the gradient to obtain a field, but I have never been able to make sense of this kind of data. For example see:

* [Ceres gravity from spherical harmonics from Dawn, how to get the coefficients, definitions and potential?](https://space.stackexchange.com/q/27689/12102)
* [How can I verify my reconstructed gravity field of Ceres from spherical harmonics?](https://space.stackexchange.com/q/27779/12102)

---

[![These maps of the moon show gravity anomalies measured by NASA's GRAIL mission. (Image credit: NASA/JPL-Caltech/CSM) ](https://i.sstatic.net/CEllE.jpg)](https://i.sstatic.net/CEllE.jpg)

> These maps of the moon show gravity anomalies measured by NASA's GRAIL mission. (Image credit: NASA/JPL-Caltech/CSM)

From Space.com's [Mystery of Moon's Lumpy Gravity Explained](https://www.space.com/21364-moon-gravity-mascons-mystery.html), click for full size.

The [Bouguer anomaly plots](https://en.wikipedia.org/wiki/Bouguer_anomaly) show elevated gravity as much as about 600 mGal. Wikipedia's [Gal (unit)](https://en.wikipedia.org/wiki/Gal_(unit)) explains what this means; 0.6 Gal is 0.006 m/s2 or almost a half of a percent of the Moon's average gravity, and they extend over quite large distances.

> The gal (symbol: Gal), sometimes called galileo after Galileo Galilei, is a unit of acceleration used extensively in the science of gravimetry. The gal is defined as 1 centimeter per second squared (1 cm/s2). The milligal (mGal) and microgal (µGal) are respectively one thousandth and one millionth of a gal.

**Tags:** open-data, moon, gravity

> **Comments:**
> **username_0**: @a_donda orbits are in 3D, I need to calculate everywhere from the surface up to hundreds of kilometers above. In free space gravity obeys the Laplace equation ($\nabla^2 \phi=0)$ so the spherical harmonics work everywhere. I suppose that if I had the potential on one smooth surface I could reconstruct the coefficients for the rest of space as answers to my old question [Calculating the potential on a surface from the potential on another surface](https://physics.stackexchange.com/q/278664/83380) explain,...  
> **username_0**: @a_donda ...but they do provide the raw coefficients https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/gggrx_1200a_sha.tab Now all I need is the specific functional form they go into. There are several ways these can be normalized (see my linked questions) so I'll need the correct form. It's the exact equation with the Bessel functions in it that I need. I can calculate orbits no problem, I just need the gradient of the potential field.  
> **username_0**: Actually the `.tab` files [here](https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/) contain the coefficients and the `.lbl` files contain the explanations. It's all starting to come back to me now, but I think the question of how the coefficients are normalized still needs to be pinned down, there are several different conventions. See https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-2-edr-v1/grail_0001/document/dpsis.pdf  
> **username_0**: @a_donda please take a moment and have a look at the the two questions I've linked to in my question and their answers to see what I mean my normalization. This is a question about math and Earth Scientists do indeed use math!  
> **username_1**: With variations as small as half a percent, the gravity vector is almost perfectly radially inward if you put the moon's center at the coordinate origin. Then, if you have the magnitude of the gravity at some particular height datum (as in the data sets mentioned above), you just need to scale it with $1/r^2$ to get the value at some height above the surface.  
> **username_0**: @username_1 no that's not sufficient for me. We've already found spherical harmonic coefficients necessary to build the potential field for all points above the Moon's surface and the gradient will provide the acceleration vectors. I just need to find the correct equations and the normalizations for those coefficients.  
> **username_0**: @a_donda No. The perturbation from a mass concentration will certainly be important out to a distance of several times it's diameter, and some of these are several hundreds of kilometers in diameter. They were discovered by their substantial perturbations on objects orbiting the Moon and they even resulted in one crashing into the Moon I think. I recommend you take some time to read up on them if you are interested.  

---

## 1 answers

### Answer by username_2 on 2022-01-31. Score: 6

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](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2014JE004755).

For that, I believe we can first look at equation 108 of [this technical note](http://icgem.gfz-potsdam.de/str-0902-revised.pdf), 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](https://iopscience.iop.org/book/978-1-64327-470-6/chapter/bk978-1-64327-470-6ch1) 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](https://earthscience.stackexchange.com/questions/23480/expression-for-the-gradient-of-gravity-potential-in-spherical-coordinates), 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](http://mitgcm.org/%7Emlosch/geoidcookbook/node11.html) 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](https://ntrs.nasa.gov/api/citations/20110023121/downloads/20110023121.pdf) 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](http://icgem.gfz-potsdam.de/calcgrid) 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](https://github.com/Rafael-Ayala/astroUtils). I am leaving here an overview of what the workflow would be to reconstruct the gradient of gravity in Python:

1. 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).
2. Calculate the values of the normalized associated Legendre functions and their derivatives at the target latitude. This can be done with [function legendre](https://github.com/Rafael-Ayala/astroUtils/blob/main/legendre.py), 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!)
3. 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](https://github.com/Rafael-Ayala/astroUtils/blob/main/gravityGradient.py), 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.
4. 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](https://iopscience.iop.org/article/10.3847/1538-3881/abd414&lang=en), 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](https://astronomy.stackexchange.com/questions/48264/for-how-many-solar-system-bodies-do-recent-jpl-development-ephemerides-include-l), 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](https://ssd.jpl.nasa.gov/horizons/app.html#/), 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:
[![enter image description here](https://i.sstatic.net/pDhLV.png)](https://i.sstatic.net/pDhLV.png)

Here we have a zoom-in the differences for degree 150:
[![enter image description here](https://i.sstatic.net/zLXlY.png)](https://i.sstatic.net/zLXlY.png)

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](https://github.com/Rafael-Ayala/astroUtils/blob/main/accelSPHtoCART.py)

> **Comments:**
> **username_0**: Thank you! Where can I get the  and  coefficients properly normalized to use in these expressions?  
> **username_2**: Regarding the C and S coefficients properly normalized, the ones in the .tab files that you linked [here](https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/) are already properly normalized for use in these expressions, as indicated by the fact that the corresponding .lbl files state "The spherical harmonic coefficients are fully normalized." Yes, I have tested my solution for the case of Earth. Sorry, I see, I didn't focus it correctly! I will try to edit and add pseudocode more explicitly  
> **username_0**: I deleted most of my comment after I read more carefully. Oh, I see, I [linked to that](https://earthscience.stackexchange.com/questions/19566/how-to-make-a-python-function-return-the-moons-gravitational-force-acceleratio/23481#comment44693_19566) back in 2020!  
> **username_2**: The key is to use an algorithm that generates Pnm and their derivatives with the same normalization applied. I tried to follow up to the source of scipy's  [lpmv](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.lpmv.html), which is called by scipy's sph_harm, but it leads to FORTRAN code, which I cannot read unluckily. However, Lear's algorithm is relatively straightforward to implement following the instructions on the document by NASA  
> **username_2**: Indeed! The model that you linked seems to provide the correct coefficients. I remember also finding at some point ephemerides for some satellite orbiting the Moon. I was planning on using those ephemerides for testing a reconstruction of the gradient of gravity of Moon with this procedure! But unluckily cannot find the ephemerides source again just now... will try to get it again!  
> **username_2**: @username_0 I've had a go at coding some of the steps in Python. Hopefully that helps!  
> **username_2**: @username_0 I've updated with some additional steps and application results!  
> **username_0**: Yes you certainly have! Okay I will dig in today and go through it carefully; I need titer exactly the right amount of coffee to take it in :-)  

---

