Основываясь на очень полезном ответе @DavidHammen, я добился прогресса в восстановлении гравитационного поля Цереры по радиометрическим данным Dawn. Вопрос содержит дополнительную информацию, но здесь достаточно сказать, что я использую версию 2 данных по адресу https://sbn.psi.edu/pds/resource/dawn/dwncgravL2.html .

Ниже приведен скрипт Python, который я использовал для чтения JGDWN_CER18C_SHA.TABфайла и построения поля гравитационного потенциала. Я использую SciPy, sph_harmкоторый нормализован, и мне интересно, как убедиться, что это та же нормализация , которую предполагает нормализация системы планетарных данных НАСА во фразе (внутри JGDWN_CER18C_SHA.LBLфайла):

Некоторые детали, описывающие эту модель:
- Коэффициенты сферических гармоник полностью нормализованы.

В документации SciPy говорится (я только что скопировал html из проверки веб-страницы, и здесь он тоже волшебным образом форматируется, ура!):

Хотя я бы сравнил с опубликованной картой, и я нашел изображение ниже, на котором показаны топография и гравитация, но их нельзя сравнивать по нескольким причинам, в том числе:

1. Опубликованная карта представляет собой гравитационное ускорение, а не потенциальное
2. это график аномалии Бугера , который немного (слишком) сложен (для меня), но грубо, если я правильно понимаю, это означает, что (среди прочего, как проекция и другие поправочные термины) это гравитационное ускорение, оцененное на эллипсоиде поверхности, а не по фиксированному радиусу. Конечно, у него также (по крайней мере) удален термин «монополь».

Я понимаю, что для того, чтобы получить карту величины ускорения , вам нужно начать с градиента потенциала, но полномасштабная аномалия Бугера может иметь другие поправки, выходящие далеко за рамки того, что мне нужно понять прямо сейчас. На самом деле то, что мне нужно, это смоделировать низкую перицентрическую орбиту Dawn .

Вопрос: Вместо того, чтобы сравнивать с этим графиком аномалии Бугера, я хотел бы как-то напрямую сравнить мою реконструкцию потенциала. Как я могу это сделать? Как я могу это проверить?

" дополнительный балл: " Дают ли коэффициенты PDS уменьшенный потенциал (энергия на единицу массы) с единицами измерения км ^ 2 / с ^ 2, а не м ^ 2 / с ^ 2?

ниже: Из Park et al. Частично дифференцированный интерьер для (1) Цереры, выведенный из ее гравитационного поля и формы , том Nature 537, страницы 515–517 (22 сентября 2016 г.) https://doi.org/10.1038/nature18955

ниже: я построил график U для n ≥ 5 и 4, потому что члены младшего порядка подавляют график r = Rref. Конечно, это (вероятно, одна из нескольких причин), почему люди используют графики Бугера и оценивают их на эллипсоиде.

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