Какие переменные необходимы для расчета времени обращения в простой подкове?

Это зависит от точности, с которой вы хотите работать. Для нулевого порядка, как указано в Мюррее и Дермотте, «Динамика Солнечной системы», глава 3, вы можете сделать следующее:

• Контуры нулевой скорости, нанесенные на ваше изображение, не будут совпадать с орбитами частиц с бесконечной точностью, но они являются хорошим приближением нулевого порядка для объектов с низким эксцентриситетом относительно звезды ($v_{r}/v_{\theta} \ll 1$)
• Частица на этих орбитах находится на регулярной кеплеровской орбите радиуса$r_{\rm H}$, вне досягаемости гравитационного влияния планеты на$r_{\rm P}$. Таким образом, чтобы получить самый большой «кусок» времени прохождения одной части подковообразной орбиты, будь то внутренней или внешней части, вы можете работать с относительными скоростями и допущением кеплеровских скоростей.
• Требуется осторожность в отношении того, какое орбитальное время вас интересует: если$a$— расстояние Земли до Солнца, а подкова —$d$от идеально круговой орбиты, таким образом чередуя расстояния$a\pm d$, то время обращения в системе покоя Солнца будет$v_{-}=\sqrt{\frac{GM}{a-d}}$, на большое количество витков до близкого сближения, а после этого$v_{+}=\sqrt{\frac{GM}{a+d}}$.
• Следовательно, упрощенный случай, когда спутник находится на собственной кеплеровской орбите, в большинстве случаев верен. Вооружившись этими знаниями, мы можем аппроксимировать планетоцентрическое время повторения как$t_{\rm rec}=\frac{2\pi a}{v_{K} - v_{-}}$как простое время наверстывания между объектами на разных орбитах.
• Относительная скорость$v_{\rm rel} = v_{K} - v_{-} = \sqrt{\frac{GM}{a}}-\sqrt{\frac{GM}{a-d}}$может быть расширен в пределах$d/a\ll1$в$v_{\rm rel} \approx \frac{1}{2} \sqrt\frac{GM}{a} \frac{d}{a}$и поэтому я исхожу из этого$t_{rec}\approx \frac{4\pi}{\sqrt{GM}} \frac{a^{5/2}}{d}$. Как и ожидалось, время цикла расходится для$d=0$, как и в случае совместной орбиты с Землей, это время должно быть бесконечным. Верхний предел на$d$не может быть получено из этого, нужно было бы обратиться к полному решению для этого.

Просто из любопытства я подставил несколько значений в эту формулу и быстро написал что-то на питоне:

import numpy as np

#Basic physics quantities
G      = 6.678e-8 #cgs units
pi     = 3.141592
navo   = 6e23 # particles per mole
sigma  = 5.67e-5   #erg cm-2 s-1 K-4
kb     = 1.38e-16  #erg/K
km     = 1e5 #kilometers in cm
mearth = 5.98e27  #g
msun   = 2.0e33   #g
au     = 1.49e13  #cm
yr   = 365*24*3600
rearth   = 6370e5
rjupiter = 74000*km

#
# Returns the approximate horseshoe-cycle time in the reduced 3body problem
# Masses of bodies: m0>>m1>>m2
# Semimajor axis distance is from m0 to m1, radial distance is a(m0->m1)-a(m0->m2)
#
def hs_cycle(mcentral, semimajor_axis, radial_distance):
    return 4*pi/np.sqrt(G*mcentral)*semimajor_axis**(5./2.)/radial_distance/yr

#
# https://en.wikipedia.org/wiki/(419624)_2010_SO16 around the Sun
#
# Quoted cycle time ~350 years, with d=0.004 AU
#
print("Predicted 2010_SO16 cycletime [years] = ", hs_cycle(msun, 1*au,0.004*au), " predicted = 350 yrs")

#
# Janus/Epimetheus around Saturn
#
# a = 151410 km, d = 25 km, as stated in https://en.wikipedia.org/wiki/Epimetheus_(moon)
# Quoted cycle time = 8 years (from comments)
#
print("Predicted Janus/Epimetheus cycletime [years] = ", hs_cycle(95*mearth, 151410*km,50*km), " predicted = 4 yrs")


#
# 3753 Cruithne
#
# a = 1 AU and semimajor axis difference from https://en.wikipedia.org/wiki/3753_Cruithne
# Quoted cycle time = 770 years
#
print("Predicted 3753 Cruithne cycletime [years] = ",hs_cycle(msun, 1*au, (1.0-0.99774)*au), " predicted = 770 yrs")

и результаты, которые я получаю,

Predicted 2010_SO16 cycletime [years] =  495.7747141830971  predicted = 350 yrs
Predicted Janus/Epimetheus cycletime [years] =  11.542076781209305  predicted = 8 yrs
Predicted 3753 Cruithne cycletime [years] =  877.4773702355546  predicted = 770 yrs

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

Любой, кто найдет более согласованные значения, может комментировать.

Книга Эрнеста В. Брауна « О новом семействе периодических орбит в задаче трех тел»: (табл. 6, 7.) в MNRAS, 71, (5), стр. 438–454, опубликованная 10 марта 1911 г. , кажется, там, где были подковообразные орбиты. первым предложил. (Доступно и здесь ). Начинается:

Известно четыре астероида, которые, по-видимому, колеблются вокруг одной или другой вершины двух равносторонних треугольников, основанием которых является линия, соединяющая Юпитер и Солнце. Эти вершины являются хорошо известными положениями относительного равновесия. Гелиоцентрический вектор одного из этих астероидов, по-видимому, может отклоняться от своего положения равновесия на 17° *). Поэтому колебания нельзя считать очень малыми. Естественно возникает вопрос, возможны ли колебания такого рода в дугах еще большей протяженности; и если да, то каким образом наиболее удобно получить орбиты.

* Л. Дж. Линдерс, Архивфор мат., Аст. оч Фыс., Со. Вет. Ак. я Стокгольм, Bd. 4, № 20.

---

---

Я сделаю несколько подковообразных орбит в формализме круговой ограниченной задачи трех тел и построю их на Python, а затем сравним с оценкой синодического периода, описанной в ответе @AtmosphericPrisonEscape .

Вкратце: Хорошее качественное совпадение, без сюрпризов!

---

Краткое изложение математики CR3BP в безразмерных единицах. Расстояние между двумя телами равно 1, как и гравитационная постоянная. Они вращаются вокруг общего центра масс по круговым орбитам с периодом$2 \pi$. Это легче визуализировать и вычислить, если вы сделаете это во вращающейся системе координат, так что две массы будут фиксированными. Третье тело на позиции$x, y, z$считается, что гравитационное воздействие на первые два не оказывает.

$$\mu = \frac{m_2}{m_1 + m_2}$$

$$x_1 = -\mu$$ $$x_2 = 1-\mu$$

$$r_1 = \sqrt{(x-x1)^2 + y^2 + z^2}$$ $$r_2 = \sqrt{(x-x2)^2 + y^2 + z^2}$$

Якоби Энерджи$C$— сохраняющаяся величина в этой вращающейся системе отсчета:

$$C = x^2 + y^2 + 2\frac{1-\mu}{r_1} + 2\frac{\mu}{r_2} - (\dot{x}^2 + \dot{y}^2 + \dot{z}^2)$$

где$x^2 + y^2$является псевдопотенциалом. Если вы установите условия, зависящие от скорости$(\dot{x}^2 + \dot{y}^2 + \dot{z}^2)$к нулю, вы получаете поверхность нулевой скорости , эту поверхность, которая вставлена ​​во многие/большинство вопросов об орбитах трех тел. Эти графики неприменимы, когда объект движется, поэтому вы не можете накладывать на них орбиты!

Ускорение, ощущаемое третьим телом в этой вращающейся системе отсчета, имеет ожидаемое$1/r^2$силы и псевдосила , зависящая от скорости, которая не является реальной, но объясняет тот факт, что система отсчета вращается, а не инерционна.

$$\ddot{x} = x + 2\dot{y} - \frac{(1-\mu)(x+\mu)}{r_1^3} - \frac{\mu(x-1+\mu)}{r_2^3}$$ $$\ddot{y} = y - 2\dot{x} - \frac{(1-\mu)y}{r_1^3} - \frac{\mu y}{r_2^3}$$ $$\ddot{z} = -\frac{(1-\mu) z}{r_1^3} - \frac{\mu z}{r_2^3}$$

---

Вот некоторые расчеты. Я выбрал$\mu = 0.001$что довольно близко к ситуации Юпитера и Солнца. Я выбрал массив начальных точек в точке, противоположной$m_2$примерно в$x=-1$ но это не то, что я действительно сделал. Что я действительно сделал, так это выбрал кучу начальных скоростей.$-0.08 < \dot{y} < 0.08$и для каждого я рассчитал позицию на$x$ось рядом$x=-1$где ускорение в$x$направление было нулевым.

Это придает решениям небольшую начальную симметрию, но гало-орбиты неровные, извилистые и не всегда такие стабильные, так что в этом нет необходимости.

Я распространял каждую орбиту, пока она не возвращалась в ту же область, и останавливал ее, когда она пересекала ось x, создавая семейство полупериодов.

Короче говоря, метод, показанный в ответе @AtmosphericPrisonEscape для оценки времени цикла путем вычисления синодического периода в инерциальной системе отсчета, довольно хорошо согласуется с этими гало-орбитами, и это не должно быть очень удивительно!

вверху: полупериоды некоторых шатких подковообразных орбит

выше: время до первого пересечения оси X тех же шатких подковообразных орбит, используемое для расчета времени полупериода.

выше: время цикла из этого расчета (черные точки) по сравнению с синодическим методом оценки периода (красные точки). Хорошее качественное согласие. Также начальные скорости y в каждой начальной точке по x.

ниже: скрипт Python для этих графиков.

def x_acc(x, ydot):
    r1    = np.abs(x-x1)
    r2    = np.abs(x-x2)
    xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
    return xddot

def C_calc(x, y, z, xdot, ydot, zdot):
    r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
    r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
    C = (x**2 + y**2 + 2.*(1-mu)/r1 + 2.*mu/r2 - (xdot**2 + ydot**2 + zdot**2))
    return C

def deriv(X, t): 
    x, y, z, xdot, ydot, zdot = X
    r1 = np.sqrt((x-x1)**2 + y**2 + z**2)
    r2 = np.sqrt((x-x2)**2 + y**2 + z**2)
    xddot = x + 2*ydot  -  ((1-mu)/r1**3)*(x+mu) - (mu/r2**3)*(x-(1-mu))
    yddot = y - 2*xdot  -  ((1-mu)/r1**3)*y      - (mu/r2**3)*y
    zddot =             -  ((1-mu)/r1**3)*z      - (mu/r2**3)*z
    return np.hstack((xdot, ydot, zdot, xddot, yddot, zddot))

# http://cosweb1.fau.edu/~jmirelesjames/hw4Notes.pdf

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

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

mu = 0.001

x1 = -mu
x2 = 1. - mu

x = np.linspace(-1.4, 1.4, 1201)
y = np.linspace(-1.4, 1.4, 1201)

Y, X = np.meshgrid(y, x, indexing='ij')
Z    = np.zeros_like(X)

xdot, ydot, zdot = [np.zeros_like(X) for i in range(3)]

C = C_calc(X, Y, Z, xdot, ydot, zdot)
C[C>8] = np.nan

if True:
    plt.figure()
    plt.imshow(C)
    plt.colorbar()
    levels = np.arange(2.9, 3.2, 0.04) 
    CS = plt.contour(C, levels,
                 origin='lower',
                 linewidths=2) 
    plt.show()

ydot0s   = np.linspace(-0.08, 0.08, 20)
x0ydot0s = []
for ydot0 in ydot0s:
    x0, infob =  brentq(x_acc, -1.5, -0.5, args=(ydot0), xtol=1E-11, rtol=1E-11,
                           maxiter=100, full_output=True, disp=True)
    x0ydot0s.append((x0, ydot0))

states = [np.array([x0, 0, 0, 0, ydot0, 0]) for (x0, ydot0) in x0ydot0s]

times  = np.arange(0, 150, 0.01)

results = []
for X0 in states:
    answer, info = ODEint(deriv, X0, times, atol = 1E-11, full_output=True)
    results.append(answer.T.copy())

resultz = []
for x0ydot0, thing in zip(x0ydot0s, results):
    y     = thing[1]
    check = y[2:]*y[1:-1] < 0
    zc    = np.argmax(y[2:]*y[1:-1] < 0) + 1
    if zc > 10:
        resultz.append((thing, zc, x0ydot0))

if True:
    plt.figure()
    hw = 1.6
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x, y = thing[:2,:zc]
        plt.plot(x, y)
    plt.xlim(-hw, hw)
    plt.ylim(-hw, hw)
    plt.plot([x1], [0], 'ok')
    plt.plot([x2], [0], 'ok')
    plt.show()

if True:
    plt.figure()
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x, y = thing[:2]
        plt.plot(times[:zc], y[:zc])
    plt.show()

if True:
    plt.figure()
    for j, (thing, zc, x0ydot0) in enumerate(resultz):
        x0, ydot0 = x0ydot0
        cycle_time = 2. * times[zc] / twopi
        ratio = abs(x0/x2)
        T_simple_model = twopi * abs(x0/x2)**1.5
        T_synodic_simple_model = 1. / (1. - twopi/T_simple_model) # https://astronomy.stackexchange.com/a/25002/7982
        plt.subplot(2, 1, 1)
        plt.plot(x0, cycle_time, 'ok')
        plt.plot(x0, abs(T_synodic_simple_model), 'or')
        plt.subplot(2, 1, 2)
        plt.plot(x0, ydot0, 'ok')
    plt.subplot(2, 1, 1)
    plt.xlabel('x0', fontsize=16)
    plt.ylabel('cycle times (periods)', fontsize=16)
    plt.subplot(2, 1, 2)
    plt.xlabel('x0', fontsize=16)
    plt.ylabel('ydot0', fontsize=16)
    plt.show()