Ceresの再構成された重力場を球面調和関数から検証するにはどうすればよいですか?

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

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

このモデルについて説明するいくつかの詳細は次のとおりです:
   – 球面調和係数は完全に正規化されています。

SciPyのマニュアルには、(私はちょうどウェブページの検査からhtmlをコピーしました。それもここでも魔法のようにフォーマットされています)

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

私は出版された地図と比較するつもりですが、私は地形と重力である下の画像を見つけましたが、これらは以下を含むいくつかの理由で比較することはできません:

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

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"
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()
ベストアンサー

処方

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

正規化の有無

私は数週間前に同様の問題に直面しました。正規化の点で、少なくとも地球、月、およびその他のカーネルの場合は、 PDS
Geosciencesノード
の場合、SHADRファイルには正規化された球面調波が含まれます。同様に、EGM2008タイドフリーモデルには正規化された高調波も含まれています。これは、階乗計算で高調波を再計算する必要がないので、計算が大幅に簡素化されます。

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

返信を残す

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