Wie berechnet man die Planeten und Monde jenseits der Gravitationskraft von Newton?

"Frage: Wie berechnet man die Planeten und Monde jenseits der Gravitationskraft von Newton?"

Uhoh, Ihr Kommentar hat weitere Quellen dazu eingeladen. (Ein großes Lob übrigens für all die Arbeit und die interessanten Ergebnisse, die Sie in Ihrer eigenen Antwort gegeben haben.)

Haben Sie gesehen, was Steve Moshier (SL Moshier) Anfang der 1990er gemacht hat?

Er gab eine vollständige Replikation des Physikmodells für die (damals standardmäßige) JPL numerisch integrierte Ephemeride DE200/LE200 (verwendet als Grundlage der Sonnensystemdaten des Astronomical Almanac in den Jahren 1984-2002), einschließlich des vollständigen Quellcodes C zusammen mit einem geeigneten Integrator &c), wodurch auch eine Erweiterung des 250-Jahre-Zeitbereichs ermöglicht wird, der dann für DE200 von JPL veröffentlicht wurde. Moshiers Integration wurde unabhängig mit der JPL-Integration über den 250-jährigen gemeinsamen Teil des Zeitbereichs von J Chapront am Pariser Observatorium verglichen, der herausfand, dass für die fünf äußeren Planeten „die Diskrepanzen niemals über 4,10^-7 Bogen hinausgehen -Sekunde, was überreichlich ist", und im schlimmsten Fall (Mond) überstiegen Längenunterschiede nie 0,008 über das 250-Jahre-Zeitintervall von DE200.

Um das Physikmodell zu vervollständigen, damit es dem damaligen Standard entspricht, musste Moshier Informationen/Daten suchen, die über das hinausgingen, was veröffentlicht worden war, und er bestätigte die zusätzlichen Daten von JPL und anderswo.

Soweit mir bekannt ist, ist dies die einzige standardmäßige Ephemeriden-Integration des Sonnensystems, für die eine vollständige und funktionsfähige Implementierung öffentlich zugänglich gemacht wurde, und dies scheint sie zu einer bemerkenswerten und sogar historisch bedeutsamen Arbeit zu machen.

Die Verweise auf Moshiers DE200-Integration (ausgeführt als „DE118“ im Referenzrahmen von 1950 und dann gedreht) sind:

(Überblick über die Arbeit in): Moshier, SL (1992), "Comparison of a 7000-year lunar ephemeris with analytics theory", Astronomy and Astrophysics 262, 613-616: unter http://adsabs.harvard.edu/abs /1992A%26A...262..613M .

Die wichtigen Implementierungsdetails mit vollständigem (C-)Quellcode sind nicht in dem Papier von 1992 enthalten, aber noch verfügbar (bis zu diesem Schreiben im März 2018) auf der Website des Autors unter http://www.moshier.net , insbesondere in Dateien

http://www.moshier.net/de118i.zip ;

http://www.moshier.net/de118i-1.zip ;

und http://www.moshier.net/de118i-2.zip ;

mit Kommentar in http://www.moshier.net/system.html .

(Diese Dateien stammen von 1993 bis 2004, die späteren Modifikationen sollten das Modell nicht ändern, aber die Syntax für weitere Compiler aufnehmen, Kommentare hinzufügen und zusätzliche Optionen wie die Einführung weiterer Körper in die Integration usw. ermöglichen.)

Die "primäre zusammenfassende Referenz" für das Physikmodell war:

Newhall, XX, EM Standish und JG Williams (1983), „DE102: eine numerisch integrierte Ephemeride des Mondes und der Planeten über vierundvierzig Jahrhunderte“, Astronomy and Astrophysics 125, 150–167, unter http://adsabs.harvard .edu/abs/1983A%26A...125..150N .

Die Rotationsmatrix zum Ändern des Referenzrahmens 1950 -> 2000 stammt von Standish, EM (1982), "Orientation of the JPL Ephemerides, DE200/LE200, to the Dynamical Equinox of J2000", Astronomy and Astrophysics 114, 297-302, at http://adsabs.harvard.edu/abs/1982A%26A...114..297S .

Die unabhängige Überprüfung wird in erwähnt

Chapront, J. (1995), "Darstellung planetarer Ephemeriden durch Frequenzanalyse. Anwendung auf die fünf äußeren Planeten." Astronomy and Astrophysics Suppl., v.109, 181-192: unter http://adsabs.harvard.edu/abs/1995A%26AS..109..181C .

Lassen Sie uns eine Annäherung hinzufügen, um einige der Effekte der Allgemeinen Relativitätstheorie (GR) zu berücksichtigen – zumindest für Körper, die in der Nähe der massereichen Sonne kreisen – und beginnen, sie zu betrachten$J_2$der Multipolterm niedrigster Ordnung für das Gravitationspotential eines Körpers jenseits des Monopolterms$-GM/r$.

Newton:

Die Beschleunigung eines Körpers im Gravitationsfeld eines anderen Körpers mit Standard-Gravitationsparametern $GM$kann geschrieben werden:

$$\mathbf{a_{Newton}} = -GM \frac{\mathbf{r}}{|r|^3},$$

wo$r$ist der Vektor vom Körper$M$an den Körper, dessen Beschleunigung berechnet wird. Denken Sie daran, dass in der Newtonschen Mechanik die Beschleunigung jedes Körpers nur von der Masse des anderen Körpers abhängt, obwohl die Kraft von beiden Massen abhängt, da sich die erste Masse durch aufhebt$a=F/m$.

Allgemeine Relativitätstheorie (ungefähr):

Obwohl ich mit GR nicht vertraut bin, werde ich eine Gleichung empfehlen, die gut zu funktionieren scheint und von mehreren Links unterstützt zu werden scheint. Es ist eine ungefähre relativistische Korrektur der Newtonschen Schwerkraft, die in Simulationen der Orbitalmechanik verwendet wird. Unter den folgenden Links sehen Sie verschiedene Formen, die meisten mit zusätzlichen Begriffen, die hier nicht aufgeführt sind:

• https://physics.stackexchange.com/q/313146/83380
• Gl. 1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf
• Gl. 27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
• Gl. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
• siehe auch David Hammens Diskussion/Ratschlag in dieser Antwort .
• Gl. 3.11 in http://adsabs.harvard.edu/full/1994AJ....107.1885S (klicken Sie auf "Drucken"), siehe diese Antwort
• Gl. A.8 in GF Rubilar1 & A. Eckart 2001 Periastron-Verschiebungen von Sternbahnen in der Nähe des Galaktischen Zentrums , außer dass die$\mathbf{r}/r^3$Newton ist enthalten und wird nicht separat behandelt.
• Gl. 2 in M. Parsa et al. 2017 Untersuchung der relativistischen Bewegung der Sterne in der Nähe des supermassiven Schwarzen Lochs im galaktischen Zentrum , außer dass die$\mathbf{r}/r^3$Newton ist enthalten und wird nicht separat behandelt.

Die folgende Annäherung sollte dem Newtonschen Term hinzugefügt werden:

$$\mathbf{a_{GR}} = GM \frac{1}{c^2 |r|^3}\left(4 GM \frac{\mathbf{r}}{|r|} - (\mathbf{v} \cdot \mathbf{v}) \mathbf{r} + 4 (\mathbf{r} \cdot \mathbf{v}) \mathbf{v} \right),$$

Oblate ($J_2$nur):

Ich verwende nur die Mathematik aus dem Wikipedia-Artikel über das Geopotentialmodell mit einer sehr wichtigen Annäherung, die man sich merken sollte. Ich gehe davon aus, dass die Abflachung in der Ebene der Ekliptik liegt – dass die Rotationsachse des abgeflachten Körpers in der Ebene liegt$\mathbf{\hat{z}}$Richtung, senkrecht zur Ekliptik. Vergessen Sie nicht, dass dies eine Annäherung ist! Die vollständigen Vektorgleichungen sind chaotischer als diese, ich werde versuchen, zurückzukommen und dies zu aktualisieren, sobald ich sicher bin, dass ich es richtig gemacht habe. In der Zwischenzeit hier eine Annäherung:

$$\mathbf{r} = x \mathbf{\hat{x}} + y \mathbf{\hat{y}} + z \mathbf{\hat{z}}$$

$$a_x = J_2 \frac{x}{|r|^7} (6z^2 - 1.5(x^2+y^2))$$

$$a_y = J_2 \frac{y}{|r|^7} (6z^2 - 1.5(x^2+y^2))$$

$$a_z = J_2 \frac{z}{|r|^7} (3z^2 - 4.5(x^2+y^2))$$

Dem Newtonschen Term ist folgendes hinzuzufügen :

$$\mathbf{a_{J2}} = a_x \mathbf{\hat{x}} + a_y \mathbf{\hat{y}} + a_z \mathbf{\hat{z}}$$

Gezeitenkräfte:

Noch komplizierter wird es, wenn man Terme betrachtet, die nicht-sphärische Massenverteilungen in beiden Körpern gleichzeitig beinhalten, egal ob sie statisch sind oder nicht. An diesem Punkt ist es wahrscheinlich notwendig, die Bücher zu schlagen.

---

Hier ist ein Testlauf. Ich habe mit heruntergeladenen Daten von JPLs Horizons verglichen . Für die äußeren Planeten verwende ich die Horizons-Daten für das Baryzentrum jedes Planeten, was sicherstellt, dass es die Geschwindigkeit für das Massenzentrum des Planeten plus all seiner Monde ist. Ich habe die Korrektur nicht zu den Massen des Planeten hinzugefügt, aber das ist ein viel kleinerer Effekt, da er nur die Bewegung anderer, entfernter Körper beeinflusst.

Stellen Sie für die Erddaten sicher, dass Sie das Geozentrum der Erde und den Mond separat herunterladen (nicht das Baryzentrum Erde-Mond). Denken Sie für die äußeren Planeten daran, die Schwerpunkte herunterzuladen.

Ich habe ein Jahr lang integriert, und alles ist innerhalb von etwa einem Kilometer von den Horizons-Daten entfernt, mit Ausnahme des Erdmonds. Das ist keine Überraschung angesichts all der intimen Gezeiteneffekte zwischen diesen beiden. Erde hinzufügen$J_2$Das vom Mond gefühlte Potenzial hilft nur teilweise, es ist wirklich nicht der richtige Weg, besonders wenn man bedenkt, dass die Erdachse (und damit die Abflachung) so weit von der Ekliptik entfernt ist.

Dies ist also insgesamt ein Beispiel dafür, dass Sie einer wirklich ernsthaften JPL-Simulation umso näher kommen, je mehr Physik Sie einsetzen! Dies ist der absolute Abstand zwischen den hier simulierten Positionen und der Horizons-Ausgabe von 2017-01-01 00:00bis 2018-01-01 00:00. Darauf folgt das Python-Skript, mit dem ich es generiert habe.

---

Basierend auf der Diskussion der Steifheit in den Kommentaren unten ist hier ein kurzes Diagramm der Schrittgröße im Vergleich zur Zeit. Ich denke, die Steifheit kommt vom Erde-Mond-System, diese häufigen Exkursionen sehen ein bisschen wie die Fehlerexkursionen der Erde und des Mondes aus. Ich denke, ich werde versuchen, das Problem auf dimensionslose Einheiten umzuskalieren, im Moment gilt die numerische Toleranz für alle Geschwindigkeiten und Positionen gleichermaßen, keine gute Idee.

def deriv_Newton_Only(X, t):
    x,  v  = X.reshape(2, -1)
    xs, vs = x.reshape(-1, 3), v.reshape(-1, 3)
    things = zip(bodies, xs, vs)

    accs, vels = [], []
    for a, xa, va in things:
        acc_a = np.zeros(3)
        for b, xb, vb in things:
            if b != a:
                r = xa - xb
                acc_a += -b.GM * r * ((r**2).sum())**-1.5
        accs.append(acc_a)
        vels.append(va)

    return np.hstack((np.hstack(vels), np.hstack(accs)))

def deriv_sunGRJ2EarthJ2(X, t):
    x,  v  = X.reshape(2, -1)
    xs, vs = x.reshape(-1, 3), v.reshape(-1, 3)
    things = zip(bodies, xs, vs)

    accs, vels = [], []
    for a, xa, va in things:
        acc_a = np.zeros(3)
        for b, xb, vb in things:
            if b != a:
                r = xa - xb
                acc_a += -b.GM * r * ((r**2).sum())**-1.5

        if a.do_SunGR and not a.name == 'Sun':

            a.flag0 = True

            r   = xa - xs[0]
            v   = va - vs[0]
            rsq = (r**2).sum()
            rm3 = rsq**-1.5
            rm1 = rsq**-0.5

            # https://physics.stackexchange.com/q/313146/83380
            # Eq.    1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf
            # Eq.   27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
            # Eq. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
            # Eq. 3.11 in http://adsabs.harvard.edu/full/1994AJ....107.1885S
            
            term_0 = Sun.GM / (clight**2) * rm3
            term_1 = 4.*Sun.GM * r * rm1
            term_2 =   -np.dot(v, v) * r
            term_3 = 4.*np.dot(r, v) * v

            accGR = term_0*(term_1 + term_2 + term_3)
            acc_a += accGR
            
        if a.do_SunJ2 and not a.name == 'Sun':

            a.flag1 = True

            r = xa - xs[0] # position relative to Sun
            x,   y,   z   = r
            xsq, ysq, zsq = r**2

            rsq = (r**2).sum()
            rm7 = rsq**-3.5

            # https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
            accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))

            accJ2  = J2s * np.hstack((accJ2x, accJ2y, accJ2z))
            acc_a += accJ2
            
        if a.do_EarthJ2 and not a.name == 'Earth':

            a.flag2 = True

            r = xa - xs[3] # position relative to Earth
            
            x,   y,   z   = r
            xsq, ysq, zsq = r**2

            rsq = (r**2).sum()
            rm7 = rsq**-3.5

            # https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
            accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))

            accJ2  = J2e * np.hstack((accJ2x, accJ2y, accJ2z))
            acc_a += accJ2
            
        accs.append(acc_a)
        vels.append(va)

    return np.hstack((np.hstack(vels), np.hstack(accs)))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

names = ['Sun', 'Mercury', 'Venus',
         'Earth', 'Moon', 'Mars',
         'Ceres', 'Pallas', 'Vesta',
         'Jupiter', 'Saturn', 'Uranus',
         'Neptune']

GMsDE430 = [1.32712440040944E+20, 2.203178E+13,  3.24858592E+14,
        3.98600435436E+14,    4.902800066E+12,  4.2828375214E+13,
        6.28093938E+10,       1.3923011E+10,    1.7288009E+10, 
        1.267127648E+17,      3.79405852E+16,   5.7945486E+15,
        6.83652719958E+15 ] # https://ipnpr.jpl.nasa.gov/progress_report/42-178/178C.pdf

# for masses also see ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck

# https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/
# https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de436s.bsp.lbl
# https://astronomy.stackexchange.com/questions/13488/where-can-i-find-visualize-planets-stars-moons-etc-positions
# https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/jup310.cmt
# ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck

GMs = GMsDE430

clight = 2.9979E+08 # m/s

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]

# J2 values
J2_sun = 2.110608853272684E-07  # unitless
R_sun  = 6.96E+08 # meters
J2s    = J2_sun * (GMs[0] * R_sun**2)   # is there a minus sign?

J2_earth = 1.08262545E-03  # unitless
R_earth  = 6378136.3 # meters
J2e      = J2_earth * (GMs[3] * R_earth**2)   # is there a minus sign?

JDs, positions, velocities, linez = [], [], [], []

use_outer_barycenters = True

for name in names:

    fname = name + ' horizons_results.txt'

    if use_outer_barycenters:
        if name in ['Jupiter', 'Saturn', 'Uranus', 'Neptune']:
            fname = name + ' barycenter horizons_results.txt'

    with open(fname, 'r') as infile:

        lines = infile.read().splitlines()

        iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
        iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]

        # print name, iSOE, iEOE, lines[iSOE], lines[iEOE]

        lines = lines[iSOE+1:iEOE]

        lines = [line.split(',') for line in lines]
        JD  = np.array([float(line[0]) for line in lines])
        pos = np.array([[float(item) for item in line[2:5]] for line in lines])
        vel = np.array([[float(item) for item in line[5:8]] for line in lines])

        JDs.append(JD)
        positions.append(pos * 1000.)   # km   to m
        velocities.append(vel * 1000.)  # km/s to m/s
        linez.append(lines)

JD = JDs[0] # assume they are identical

class Body(object):
    def __init__(self, name):
        self.name = name

bodies = []
for name, GM, pos, vel in zip(names, GMs, positions, velocities):
    
    body = Body(name)
    bodies.append(body)
    
    body.GM = GM
    
    body.daily_positions  = pos
    body.daily_velocities = vel
    
    body.initial_position = pos[0]
    body.initial_velocity = vel[0]

x0s = np.hstack([b.initial_position for b in bodies])
v0s = np.hstack([b.initial_velocity for b in bodies])

X0  = np.hstack((x0s, v0s))

ndays   = 365
nperday = 144

time = np.arange(0, ndays*24*3600+1, 24*3600./nperday)
days = time[::nperday]/(24*3600.)

for body in bodies:
    body.do_SunGR   = False
    body.do_SunJ2   = False
    body.do_EarthJ2 = False
    body.flag0      = False
    body.flag1      = False
    body.flag2      = False

Sun, Mercury, Venus, Earth, Moon, Mars = bodies[:6]
Ceres, Pallas, Vesta = bodies[6:9]
Jupiter, Saturn, Uranus, Neptune = bodies[9:]

Mercury.do_SunGR = True
Venus.do_SunGR   = True
Earth.do_SunGR   = True
Moon.do_SunGR    = True
Mars.do_SunGR    = True
Ceres.do_SunGR   = True
Pallas.do_SunGR  = True
Vesta.do_SunGR   = True

Mercury.do_SunJ2 = True

Moon.do_EarthJ2  = True

rtol = 1E-12   # this is important!!!

answer, info = ODEint(deriv_sunGRJ2EarthJ2, X0, time,
                      rtol = rtol, full_output=True)

print answer.shape

nbodies = len(bodies)

xs, vs = answer.T.reshape(2, nbodies, 3, -1)

for body, x, v in zip(bodies, xs, vs):
    body.x = x
    body.v = v
    body.x_daily = body.x[:, ::nperday]
    body.v_daily = body.v[:, ::nperday]
    body.difference = np.sqrt(((body.x_daily - body.daily_positions.T)**2).sum(axis=0))

if True:
    for body in bodies[:6]:
        print body.name, body.flag0, body.flag1, body.flag2

if True:
    plt.figure()
    for i, body in enumerate(bodies[:12]):  # Sorry Neptune!!!
        plt.subplot(4, 3, i+1)
        plt.plot(days, 0.001*body.difference)
        plt.title(body.name, fontsize=14)
        plt.xlim(0, 365)
    plt.suptitle("calc vs JPL Horizons (km vs days)", fontsize=16)
    plt.show()

Ich möchte nur hinzufügen, dass der in der Antwort von uhoh erwähnte relativistische Korrekturterm, bei dem es sich um die "post-Newtonsche Erweiterung" auf der "1PN" -Ebene handelt, relativistische Effekte durch Einführung einer Abstoßung annähert$1/r^3$Begriff. Der Ausdruck wird von der JPL verwendet, also müssen Sie ihn verwenden, wenn Sie so nah wie möglich an die Ephemeriden herankommen wollen. Obwohl Sie die "anomale Perihelverschiebung" richtig verstehen, erhalten Sie sehr seltsame Effekte des "Aufpralls" in der starken Feldgrenze, sodass der Ausdruck wahrscheinlich hauptsächlich in den schwachen Feldern unseres Sonnensystems funktioniert. Ich habe unten einige Simulationen durchgeführt, der grüne Kreis ist der Schwarzschild-Radius und der rote Kreis ist der Radius der "innersten stabilen Kreisbahn", die sich in einem radialen Abstand von drei Schwarzschild-Radien befindet. Das sichtbare "Aufprallen" ist offensichtlich auf den abstoßenden inversen Würfelterm zurückzuführen. Mit mehr anfänglichem Drehimpuls werden die Bahnen weniger seltsam .