# How can I verify my reconstructed gravity field of Ceres from spherical harmonics?

Based on @DavidHammen's very helpful answer I've made progress reconstructing the gravity field of Ceres from the Dawn radiometric data. The question there contains further information, but suffice it to say here that I am using Version 2 of the data at https://sbn.psi.edu/pds/resource/dawn/dwncgravL2.html

Below is a Python script I've used to read the JGDWN_CER18C_SHA.TAB file and build the gravitational potential field. I am using SciPy's sph_harm which is normalized, and I am wondering how to be sure if this is the same normalization assumed by the NASA Planetary Data System's normalization in the phrase (within the JGDWN_CER18C_SHA.LBL file):

Some details describing this model are:
- The spherical harmonic coefficients are fully normalized.

The SciPy documentation says (I've just copied the html from inspection of the webpage and it magically formats here too, yay!):

$$Y^m_n(\theta,\phi) = \sqrt{\frac{2n+1}{4\pi} \frac{(n-m)!}{(n+m)!}} e^{i m \theta} P^m_n(\cos(\phi))$$

I though I would compare to a published map, and I've found the image below which is topography and gravity, but these can not be compared for several reasons, including:

1. The published map is gravitational acceleration, not potential
2. it is a plot of Bouguer anomaly which is a little (too) complicated (for me) but roughly, if I understand correctly, it means that (among other things like projection, and other correction terms) it's the gravitational acceleration evaluated on the ellipsoidal surface, rather than at a fixed radius. It's also got (at least) the monopole term removed of course.

I understand that to get an acceleration magnitude map, you need to start with the gradient of the potential, but a full blown Bouguer anomaly can have other corrections way beyond what I need to understand right now. In fact what I'm after is modeling Dawn's low periapsis orbit.

Question: Instead of comparing to this Bouguer anomaly plot, I'd like to compare my reconstruction of potential directly somehow. How can I do this? How can I check it?

"extra credit:" Do the PDS coefficients yield a reduced potential (energy per unit mass) with units of km^2/s^2 rather than m^2/s^2?

below: From Park et al. A partially differentiated interior for (1) Ceres deduced from its gravity field and shape Nature volume 537, pages 515–517 (22 September 2016) https://doi.org/10.1038/nature18955

below: I've plotted U for n ≥ 5 and 4 because the low order terms overwhelm the r = Rref plot. Of course this is (probably one of several reasons) why people use Bouguer plots and evaluate on the ellipsoid.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sph_harm

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

j = np.complex(0, 1)

fname = "JGDWN_CER18C_SHA.TAB"

with open(fname, 'r') as infile:

Rref, GM, GMerr = [float(x) for x in header_data[0:3]]
Order_0, Order_1, normalization_state = [int(x) for x in header_data[3:6]]

if normalization_state == 1:
print "coefficients are normalized"
elif normalization_state == 0:
print "coefficients are NOT normalized"
else:
print "coefficients  normalization is unclear"

h_lines = [line.split(',') for line in lines[1:]]
indices = np.array([[int(x)   for x in line[0:2]] for line in h_lines])
coeffs  = np.array([[float(x) for x in line[2:4]] for line in h_lines])

Cstars = (np.array([1, +j]) * coeffs).sum(axis=1) # make coefficient complex

ph         = np.linspace(0,  pi,    180+1)[:-1]
th         = np.linspace(0,  twopi, 360+1)[:-1]
phi, theta = np.meshgrid(ph, th, indexing='ij')

# https://docs.scipy.org/doc/scipy-1.1.0/reference/generated/scipy.special.sph_harm.html#scipy.special.sph_harm

harmonics = []

for (n, m), Cstar in zip(indices, Cstars):

Y = sph_harm(m, n, theta, phi)

harmonics.append((n, m, (Y * Cstar).real))  # 3-tuple of n, m, Y*C product

# evaluate gravitational potential
r = Rref

U_mono   = -GM/r

nmins = (5, 4)     # 5, 4, 3, 2

Us = []

for nmin in nmins:
count = 0
U = np.zeros_like(phi)
for n, m, h in harmonics:
if n >= nmin:
U += h * (Rref/r)**n
count += 1
print nmin, count
Us.append(U)

if True:
plt.figure()
for i, (nmin, U) in enumerate(zip(nmins, Us)):
plt.subplot(len(Us), 1, i+1)
plt.imshow(U, cmap='PuOr')
plt.title('U_ceres(r=' + str(round(r, 1))  + 'km), nmin = ' + str(nmin), fontsize = 16)
plt.colorbar()
plt.show()

• The scipy sph_harm package is not what you want. It uses the form of spherical harmonics favored in quantum physics. Geophysicists use a form that uses the real sine and cosine functions rather than the complex exponential, and the normalization is different. Jun 11, 2018 at 19:13
• update: I've run across the following and investigating now; Physical Geodesy, 2nd (corrected) ed. Bernhard Hofmann-Wellenhof Helmut Moritz, SpringerWienNewYork, Section 1.10 Fully normalized spherical harmonics
– uhoh
Jun 12, 2018 at 11:39
• Oh my lord... this has been a fun read man. Thank you for posting the code! As a developer this is awesome!!! Jun 12, 2018 at 18:42
• If you are trying to match gravity at the surface with spherical harmonics, good luck. It can be impossible since spherical harmonics diverge within the Brillouin sphere, that is why for landing operations polyhedron or mascons models are preferred. But maybe I have misunderstood the question. Anyway, if your final application is an orbiting body I would suggest to code manually a routine manually for the acceleration (there are analytic expression for the gradient) and check if the results matches the SciPy function (if any). Sorry, not a python user, but that is what I did in MATLAB Nov 22, 2019 at 10:25
• @uhoh Not sure I can give a full answer, but in my experience normalised gravity coefficients usually means real harmonics, 4pi normalised. I understand you want to code things yourself, but perhaps check your calculations against the SHTOOLS library which is written for potential fields work in the Earth/planetary sciences. If you have access to ellipsoid and topography info this example notebook makes it straightforward to go down the Bouguer comparison route: nbviewer.jupyter.org/github/SHTOOLS/SHTOOLS/blob/master/…
– WJB
Sep 16, 2020 at 8:56

# Formulation

Brandon A. Jones wrote an excellent summary of the different formulations and methods to compute spherical harmonics in his 2010 PhD thesis available here. Chapter 2 is especially of interest. You may also want to read Fantino & Casotto 2008 "Methods of harmonic synthesis for global geopotential models and their first-, second- and third-order gradients", which proposes a slightly simplified computation of the harmonics.

Important: the frame in which the harmonics is computed is a body fixed frame. In the case of Earth, you'll use the ECEF frame (cf. Chapter 2 of G. Xu and Y. Xu, GPS, DOI 10.1007/978-3-662-50367-6_2 -- which is open access -- for a through computation of that frame from an ECI frame using the IAU2000 framework). In the case of other bodies, and Ceres in particular, you'll need to create a body fixed frame. This GMAT documentation may be of help.

# Normalization or not

I came across a similar problem myself a few weeks ago. In terms of the normalization, and at least in the cases of the Earth, Moon, and other kernels available on PDS Geosciences Node, the SHADR files contain the normalized spherical harmonics. Similarly, the EGM2008 Tide Free model also contains normalized harmonics. This greatly simplifies the computation since you do not need to recompute harmonics with factorial math.

# Verification & validation

The verification step should be quite straight-forward, but the validation isn't. The process of verification, in my experience, consists in making sure that things are implemented correctly: one checks the coded math, and that all the block fit together. The validation step however, is more complicated. There, you want to make sure that the right thing was implemented, i.e. despite the math being right, is the correct frame used to compute the effect of the spherical harmonics on the orbiting vehicle.

As discussed above, both the SHADR and SHBDR files store the normalized coefficients. Hence, one method of validation, which I've opted for my upcoming tool, consists in rigorously comparing results of the same scenario with another known validation tool. In my case, I'm using NASA's GMAT and GMAT used STK and others. Specifically in the case of GMAT however, I do not think it's possible to import SHADR files. So you'll need to convert it to the support COF file, or find another propagation tool which also allows you to turn on and off specific dynamics (in this case you'll want harmonics turned on to a given degree and order, but only Ceres as a point mass).

• ChrisR to the rescue (again)! Okay I will give this a thorough read. Can you give me a tl;dr first though? Will I find some data points of gravitational potential at a few locations I can use as confirmation points, or is your answer that I will need to make room on my laptop, download, install under MacOS, and then learn a whole new program before I can get anything useful?
– uhoh
Jun 12, 2018 at 14:57
• Right, sorry I didn't actually get to the core of your question. If you would like to verify your implementation, I think that comparing a Bouguer plot is a good start. For validation, you would need to look at specific numerical values. There are two methods for that. The first requires that you find values from a trusted source (like a paper), ie "with this given state, we get a gravity potential of X when using degree N and order M from file Z". The second possibility is finding a propagation result, or generating it from a source like GMAT (works on Windows, Linux and Mac). Jun 12, 2018 at 16:30
• @uhoh, you may also want to reach out to the planetary scientists at Open Planetary. They are very helpful and respond quite quickly on their public Slack chat system. Jun 12, 2018 at 16:32
• Okay well I thought I'd explained why comparing to a Bouguer plot would be a bad idea for me. I've never heard of Open Planetary, what is it?
– uhoh
Jun 12, 2018 at 16:44
• Okay well I thought I'd explained why comparing to a Bouguer plot would be a bad idea for me. I've never heard of Open Planetary, what is it? Oh I see, it's a collection of head-shots of people trying to look smart, forward-thinking, and knowledgable. I think that websites which prominently display a dozen vanity head-shots rather than science tend not to be very helpful, but okay I will take a look. Thanks for that suggestion!
– uhoh
Jun 12, 2018 at 16:45