Based on @DavidHammen’s very helpful
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

ファイルを読み、重力ポテンシャルフィールドを構築するために使用したPythonスクリプトです。 SciPyの sph_harm
を正規化しています。これがNASA惑星データシステムのフレーズの正規化によって仮定されているのと同じ正規化 (

   – 球面調和係数は完全に正規化されています。


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


  1. 公開された地図は、潜在的なものではなく重力です。
  2. これはブーゲーアノーマリのプロットであり、少し複雑です(for私が正しく理解すれば、(投射やその他の補正項の中でも)固定半径ではなく、楕円面上で重力加速度が評価されることを意味します。それはもちろん、単極項が削除されている(少なくとも)。


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?

PDS係数は、m ^ 2/s ^ 2ではなくkm ^ 2/s ^

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

enter image description here

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.

enter image description here

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"
    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

if True:
    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)


Brandon A. Jones wrote an excellent summary of the different 処方s
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


私は数週間前に同様の問題に直面しました。正規化の点で、少なくとも地球、月、およびその他のカーネルの場合は、 PDS

Verification & validation


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


メールアドレスが公開されることはありません。 * が付いている欄は必須項目です