# How to make a Python function return the Moon's gravitational force (acceleration) to calculate orbits?

The Space SE question How much does it cost to land on a lunar mascon? 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:

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, click for full size.

The Bouguer anomaly plots show elevated gravity as much as about 600 mGal. Wikipedia's 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.

• Looks like somebody has already done the work, maybe not in python: pgda.gsfc.nasa.gov/products/50. But if i understand you right you could use one of the maps provided as a texture and map a texture position to the reference sphere. That could fit into a function.
– user20217
Apr 6, 2020 at 9:26
• @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 explain,...
– uhoh
Apr 6, 2020 at 11:16
• 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. Apr 6, 2020 at 21:29
• @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
Apr 8, 2020 at 19:16
• I'd really like to and will keep it in my occipetal lobe. But there is a bunch of other things before me right now ... good luck finding your spells ;-)
– user20217
Apr 8, 2020 at 19:32

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:

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, 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, 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, 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

• Thank you! Where can I get the 𝐶 and 𝑆 coefficients properly normalized to use in these expressions?
– uhoh
Jan 31 at 12:02
• Regarding the C and S coefficients properly normalized, the ones in the .tab files that you linked here 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
Jan 31 at 12:02
• I deleted most of my comment after I read more carefully. Oh, I see, I linked to that back in 2020!
– uhoh
Jan 31 at 12:05
• 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, 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
Jan 31 at 12:05
• 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
Jan 31 at 12:08