Basado en la muy útil respuesta de @DavidHammen, he progresado en la reconstrucción del campo de gravedad de Ceres a partir de los datos radiométricos de Dawn. La pregunta allí contiene más información, pero basta con decir aquí que estoy usando la versión 2 de los datos en https://sbn.psi.edu/pds/resource/dawn/dwncgravL2.html

A continuación se muestra un script de Python que he usado para leer el JGDWN_CER18C_SHA.TABarchivo y construir el campo de potencial gravitacional. Estoy usando SciPy, sph_harmque está normalizado, y me pregunto cómo estar seguro de si esta es la misma normalización asumida por la normalización del Sistema de datos planetarios de la NASA en la frase (dentro del JGDWN_CER18C_SHA.LBLarchivo):

Algunos detalles que describen este modelo son:
- Los coeficientes armónicos esféricos están completamente normalizados.

La documentación de SciPy dice (acabo de copiar el html de la inspección de la página web y también se formatea mágicamente aquí, ¡sí!):

Pensé que lo compararía con un mapa publicado, y encontré la imagen a continuación que es topografía y gravedad, pero estos no se pueden comparar por varias razones, que incluyen:

1. El mapa publicado es aceleración gravitacional , no potencial.
2. es una gráfica de la anomalía de Bouguer que es un poco (demasiado) complicada (para mí) pero aproximadamente, si entiendo correctamente, significa que (entre otras cosas como proyección y otros términos de corrección) es la aceleración gravitacional evaluada en el elipsoidal superficie, en lugar de un radio fijo. También tiene (al menos) el término monopolo eliminado, por supuesto.

Entiendo que para obtener un mapa de magnitud de aceleración , debe comenzar con el gradiente del potencial, pero una anomalía de Bouguer en toda regla puede tener otras correcciones más allá de lo que necesito entender en este momento. De hecho, lo que busco es modelar la órbita periáptica baja de Dawn .

Pregunta: En lugar de comparar con este diagrama de anomalías de Bouguer, me gustaría comparar mi reconstrucción de potencial directamente de alguna manera. ¿Cómo puedo hacer esto? ¿Cómo puedo comprobarlo?

" Crédito adicional: " ¿Los coeficientes PDS producen un potencial reducido (energía por unidad de masa) con unidades de km^2/s^2 en lugar de m^2/s^2?

abajo: De Park et al. Un interior parcialmente diferenciado para (1) Ceres deducido de su campo de gravedad y forma Nature volumen 537, páginas 515–517 (22 de septiembre de 2016) https://doi.org/10.1038/nature18955

a continuación: tracé U para n ≥ 5 y 4 porque los términos de orden bajo abruman la gráfica r = Rref. Por supuesto, esta es (probablemente una de varias razones) por qué las personas usan diagramas de Bouguer y evalúan en el elipsoide.

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