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

Question

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)]
degs, rads = 180/pi, pi/180

j = np.complex(0, 1)

fname = "JGDWN_CER18C_SHA.TAB"

with open(fname, 'r') as infile:
    lines = infile.readlines()

header_data = lines[0].split(',')

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()

Answer 1

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).