Как вычислить планеты и луны за пределами гравитационной силы Ньютона?

«Вопрос: как вычислить планеты и луны за пределами гравитационной силы Ньютона?»

Ого, ваш комментарий пригласил дополнительные источники по этому поводу. (Спасибо, кстати, за всю работу и интересные результаты, которые вы дали в своем собственном ответе.)

Вы видели, что делал Стив Мошир (SL Moshier) в начале 1990-х?

Он предоставил полную копию физической модели для (тогдашнего стандарта) численно интегрированных эфемерид JPL DE200/LE200 (использовавшихся в качестве основы для данных Солнечной системы Астрономического альманаха за 1984-2002 годы), включая полный исходный код в C вместе с подходящим интегратором и c), что также позволяет расширить 250-летний временной диапазон, опубликованный JPL для DE200. Интеграция Мошира была независимо сопоставлена ​​с интеграцией Лаборатории реактивного движения за 250-летнюю общую часть временного диапазона Дж. Шапроном из Парижской обсерватории, который обнаружил, что для пяти внешних планет «расхождения никогда не превышают 4,10 ^ -7 дуги». -секунда, которая избыточна", а в худшем случае (луна) расхождения по долготе никогда не превышали 0".008 за 250-летний интервал времени DE200.

Чтобы завершить физическую модель, чтобы она соответствовала тогдашнему стандарту, Моширу пришлось искать информацию/данные сверх того, что было опубликовано, и он признал дополнительные данные из JPL и других источников.

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

Ссылки на интеграцию Moshier DE200 (выполненную как «DE118» в системе отсчета 1950 года, а затем повернутую):

(План работы в): Мошир, С.Л. (1992), «Сравнение 7000-летней лунной эфемериды с аналитической теорией», Астрономия и астрофизика 262, 613-616: на http://adsabs.harvard.edu/abs /1992А%26А...262..613М .

Важные детали реализации с полным (C) исходным кодом отсутствуют в статье 1992 г., но все еще доступны (до написания этой статьи в марте 2018 г.) на веб-сайте автора по адресу http://www.moshier.net , особенно в файлах .

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

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

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

с комментариями на http://www.moshier.net/ssystem.html .

(Эти файлы датированы с 1993 по 2004 год, более поздние модификации заключались не в изменении модели, а в адаптации синтаксиса для дальнейших компиляторов, добавлении комментариев и разрешении дополнительных опций, таких как введение дополнительных тел в интеграцию и т. д.)

«Первичная сводная ссылка» для физической модели была:

Newhall, XX, EM Standish и JG Williams (1983), "DE102: численно интегрированная эфемерида Луны и планет, охватывающая сорок четыре века", Astronomy and Astrophysics 125, 150-167, на http://adsabs.harvard .edu/abs/1983A%26A...125..150N .

Матрица вращения для изменения системы отсчета 1950-> 2000 была взята из Стэндиша, Э.М. (1982), "Ориентация эфемерид JPL, DE200/LE200, на динамическое равноденствие J2000", Astronomy and Astrophysics 114, 297-302, at http://adsabs.harvard.edu/abs/1982A%26A...114..297S .

Независимая проверка упоминается в

Чапронт, Дж. (1995), «Представление планетарных эфемерид с помощью частотного анализа. Применение к пяти внешним планетам». Приложение по астрономии и астрофизике, т. 109, 181–192: http://adsabs.harvard.edu/abs/1995A%26AS..109..181C .

Давайте добавим приближения, чтобы учесть некоторые эффекты общей теории относительности (ОТО) — по крайней мере, для тел, вращающихся вокруг массивного Солнца, — и начнем рассматривать$J_2$мультипольный член низшего порядка для гравитационного потенциала тела после монопольного члена$-GM/r$.

Ньютон:

Ускорение тела в гравитационном поле другого тела стандартного гравитационного параметра $GM$можно написать:

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

куда$r$это вектор от тела$M$к телу, ускорение которого рассчитывается. Помните, что в ньютоновской механике ускорение каждого тела зависит только от массы другого тела , хотя сила зависит от обеих масс, потому что первая масса уравновешивается$a=F/m$.

Общая теория относительности (приблизительно):

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

• https://физика.stackexchange.com/q/313146/83380
• уравнение 1 в https://www.lpi.usra.edu/books/CometsII/7009.pdf
• уравнение 27 в https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
• уравнение 4-26 в https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
• также см. обсуждение/совет Дэвида Хаммена в этом ответе .
• уравнение 3.11 в http://adsabs.harvard.edu/full/1994AJ....107.1885S (нажмите "печатать") см. этот ответ
• уравнение A.8 в GF Rubilar1 & A. Eckart 2001 Периастральные смещения звездных орбит вблизи галактического центра , за исключением того, что$\mathbf{r}/r^3$Ньютониан включен, а не рассматривается отдельно.
• уравнение 2 в M. Parsa et al. 2017 Исследование релятивистского движения звезд вблизи сверхмассивной черной дыры в центре Галактики , за исключением того, что$\mathbf{r}/r^3$Ньютониан включен, а не рассматривается отдельно.

К ньютоновскому члену следует добавить следующее приближение:

$$\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),$$

сплюснутость ($J_2$Только):

Я просто использую математику из статьи Википедии о модели геопотенциала с очень важным для запоминания приближением; Я предполагаю, что сплюснутость находится в плоскости эклиптики — что ось вращения сплюснутого тела находится в плоскости эклиптики.$\mathbf{\hat{z}}$направление, перпендикулярное эклиптике. Не забывайте, что это приблизительно! Полные векторные уравнения более запутаны, чем это, я постараюсь вернуться и обновить это, как только буду уверен, что понял правильно. Между тем, вот приблизительное значение:

$$\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))$$

К ньютоновскому термину следует добавить следующее:

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

Приливные силы:

Это становится еще более сложным, если рассматривать термины, которые включают несферическое распределение массы в обоих телах одновременно, независимо от того, статичны они или нет. В этот момент, вероятно, необходимо поразить книги.

---

Вот пробный запуск. Я сравнил данные, загруженные с JPL Horizons . Для внешних планет я использую данные Horizons для барицентра каждой планеты, что гарантирует, что это скорость для центра масс планеты плюс все ее спутники. Я не добавлял поправку к массам планеты, но это гораздо меньший эффект, поскольку он влияет только на движение других, удаленных тел.

Для данных Земли обязательно загрузите геоцентр Земли и Луны отдельно (не барицентр Земля-Луна). Для внешних планет не забудьте загрузить барицентры.

Я интегрировал в течение года, и все находится в пределах примерно одного километра от данных Горизонтов, кроме земной Луны. Это не удивительно, учитывая все интимные приливные эффекты между этими двумя. Добавление Земли$J_2$к потенциалу, ощущаемому Луной, помогает лишь частично, это действительно неправильный способ сделать это, особенно учитывая, что ось Земли (и, следовательно, сплюснутость) так далеко от эклиптики.

Так что это в целом иллюстрация того, что чем больше физики вы добавляете, тем ближе вы можете приблизиться к действительно серьезной симуляции JPL! Это абсолютное расстояние между смоделированными позициями здесь и выводом Horizons от 2017-01-01 00:00до 2018-01-01 00:00. Далее следует скрипт Python, который я использовал для его создания.

---

Основываясь на обсуждении жесткости в комментариях ниже, вот краткий график зависимости размера шага от времени. Я думаю, что жесткость исходит из системы Земля-Луна, эти частые отклонения немного напоминают отклонения Земли и Луны. Я думаю, что я попытаюсь масштабировать проблему до безразмерных единиц, прямо сейчас числовой допуск применяется ко всем скоростям и положениям в равной степени, не очень хорошая идея.

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

Я просто хочу добавить, что релятивистский поправочный член, упомянутый в ответе uhoh, который представляет собой «постньютоновское расширение» на уровне «1PN», аппроксимирует релятивистские эффекты путем введения отталкивающего$1/r^3$срок. Выражение используется JPL, поэтому вы должны использовать его, если хотите максимально приблизиться к эфемеридам. Даже если вы правильно понимаете «аномальное смещение перигелия», вы получаете очень странные эффекты «подпрыгивания» в пределе сильного поля, поэтому выражение, вероятно, в основном работает в слабых полях нашей Солнечной системы. Я выполнил несколько симуляций ниже, зеленый кружок — это радиус Шварцшильда, а красный кружок — это радиус «самой внутренней стабильной круговой орбиты», расположенной на радиальном расстоянии в три радиуса Шварцшильда. Увиденное «подпрыгивание», очевидно, происходит из-за отталкивающего члена обратного куба. При большем начальном угловом моменте орбиты становятся менее странными .