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:
- The published map is gravitational acceleration, not potential
- 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()
