Kann „Prolateness“ durch einen J2J2J_2-Koeffizienten wie „Olateness“ charakterisiert werden?

Dies ist eine Erweiterung meines Kommentars.

Ganz grob gesagt,$J_2$ist ein Maß für Folgendes: wenn ein Punkt über dem Körperpol liegt und ein anderer Punkt über dem Äquator im gleichen Abstand liegt$R$von der Körpermitte (nicht von der Oberfläche!), wie viel ist der erste Punkt "bergauf" im Vergleich zum zweiten Punkt in Bezug auf die Schwerkraft des Körpers?

Betrachten wir zunächst den Fall zweier gleicher Punktmassen. An jedem Punkt, der nicht auf derselben Linie mit ihnen und nicht auf ihrer Symmetrieebene liegt, zeigt die Gesamtgravitation aufgrund dieser Massen nicht auf ihren Mittelpunkt, sondern weicht etwas in Richtung der nahen Masse ab (siehe Bild 1 im beigefügten Bild:$A$Und$B$sind die Massen;$O$ist ihr Mittelpunkt;$X$ist der Punkt, an dem wir die Schwerkraft messen; sorry für die schlechte Qualität). Für die Summe der Erdbeschleunigungen$\mathbf{g}_A$Und$\mathbf{g}_B$deuten auf$O$, das Verhältnis ihrer Größen$\frac{g_A}{g_B}$müßte das gleiche sein wie$\frac{|XA|}{|XB|}$(Bild 2), dh,$g_A$müsste kleiner sein als$g_B$. Aber$g_A$ist tatsächlich größer als$g_B$, also die Summe$\mathbf{g}_A+\mathbf{g}_B$zeigt auf einen Punkt auf dem Segment$OA$.

Das bedeutet, wenn Sie sich entlang eines Viertelkreises mit dem Mittelpunkt bei bewegen$O$von einem Punkt über dem "Äquator" zu einem Punkt über dem "Pol" (Abb. 3), ist der Winkel zwischen der Richtung Ihrer Bewegung und der Richtung der lokalen Schwerkraft immer (außer an den Endpunkten) kleiner als$90^\circ$, dh Sie bewegen sich "bergab". Das bedeutet es$J_2$ist in diesem Fall negativ.

Wenn wir nun eine flache Scheibe betrachten, können wir sie in Paare kleiner Abschnitte unterteilen, die symmetrisch zur Mitte angeordnet sind (Bild 4). Für je zwei solcher Schnitte zeigt die Summe der ihnen zuzuschreibenden Erdbeschleunigungen nicht auf die Mitte der Scheibe, sondern weicht zum nahen Schnitt hin ab. Dadurch zeigt die Gesamtschwerkraft der Scheibe auf einen Punkt innerhalb der nahen Halbscheibe. Wenn Sie sich also wieder entlang eines Viertelkreises mit dem Mittelpunkt in der Mitte der Scheibe von einem Punkt über dem „Äquator“ zu einem Punkt über dem „Pol“ bewegen (Bild 5), wird der Winkel zwischen der Richtung Ihrer Bewegung und die Richtung der lokalen Schwerkraft ist größer als$90^\circ$, dh Sie bewegen sich "bergauf". So$J_2$ist in diesem Fall positiv.

Zur Vereinfachung der Berechnung simulieren wir einen gestreckten Körper als zentrales Potential plus zwei kleinere Massen nördlich und südlich des Zentrums

Und simulieren wir einen abgeflachten Körper als zentrales Potential plus einen äquatorialen Ring aus Punktmassen.

Wenn wir setzen$G=1$der Einfachheit halber und berechnen Sie das Gravitationspotential$GM/r$Wir können sehen, dass es für den abgeflachten Fall am Äquator seinen Höhepunkt erreicht und für den prolaten Fall am Äquator ein Minimum hat.

Diese Wackelbewegungen sind keine reinen Sinuskurven, aber wir können sehen, dass ihr Verhalten größtenteils ähnlich ist$\sin^2(\theta)$Und$-\sin^2(\theta)$und da$J_2$ist der Koeffizient vor einem solchen Term (aus dem Geopotentialmodell zwischen den Gleichungen 9 und 10):

wir können das für einen Oblatenkörper sehen$J_2$positiv und für einen gestreckten Körper negativ.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 


def phi_central(x, m):
    """central potential"""
    x0 = np.zeros(3)
    r = np.sqrt(((x - x0)**2).sum(axis=-1))
    return m / r, x0

def phi_polar(x, m, h):
    """quadrupole potential (axial)"""
    zhat = np.array([0, 0, 1])
    x1, x2 = h * zhat, -h * zhat
    r1 = np.sqrt(((x - x1)**2).sum(axis=-1)) # top
    r2 = np.sqrt(((x - x2)**2).sum(axis=-1)) # bottom
    return m * (1/r1 + 1/r2), x1, x2
                 
def phi_equatorial(x, m, r, N=100):
    """quadrupole potential (equatorial (xy))"""
    theta = np.linspace(0, 2*np.pi, N+1)[:-1]
    # ring in the xy plane
    ring = np.vstack([r * f(theta) for f in (np.cos, np.sin, np.zeros_like)])
    r = np.sqrt(((x[..., None] - ring)**2).sum(axis=-2)) # equator
    return m / r.mean(axis=-1), r, ring


mcen, mpol, meq = 1, 0.1, 0.3
hpol, req = 0.3, 0.9

theta = np.linspace(0, 2*np.pi, 361)
R = 1.0
N = 1000

xz_plane = np.stack([R * f(theta) for f in (np.sin, np.zeros_like, np.cos)], axis=-1)

phi_c, xcen = phi_central(xz_plane, mcen)

phi_p, xp1, xp2 = phi_polar(xz_plane, mpol, hpol)

phi_e, r, ring = phi_equatorial(xz_plane, meq, req, N=N)
                    
fig, ax = plt.subplots(1, 1)
ax.plot(np.degrees(theta), phi_c)
ax.plot(np.degrees(theta), phi_p)
ax.plot(np.degrees(theta), phi_e)
ax.set_ylim(0, 1.1)
ax.set_xlabel('theta (deg)')
ax.set_ylabel('gravitational potential')
plt.show()

fig = plt.figure()
ax  = fig.add_subplot(1, 1, 1, projection='3d', proj_type = 'ortho')
ax.plot([0], [0], [0], 'ok')
x, y, z = zip(xp1, xp2)
ax.plot(x, y, z, 'or')
x, y, z = ring
ax.plot(x, y, z, '-b')
x, y, z = xz_plane.T
ax.plot(x, y, z, '-g')
plt.show()