¿Qué variables se necesitan para calcular tiempos de órbita de herradura simple?

Eso depende de la precisión con la que quieras trabajar. En orden cero, como se describe en Murray & Dermott, "Dinámica del sistema solar", Capítulo 3, puede hacer lo siguiente:

• Los contornos de velocidad cero que se trazan en su imagen no coincidirán con las órbitas de las partículas con una precisión infinita, pero son una buena aproximación de orden cero para objetos con baja excentricidad respecto a la estrella ($v_{r}/v_{\theta} \ll 1$)
• Una partícula en esas órbitas está en una órbita Kepleriana regular de radio$r_{\rm H}$, fuera del alcance de la influencia gravitatoria del planeta en$r_{\rm P}$. Por lo tanto, para obtener el "trozo" más grande de tiempo de viaje de una parte de la órbita de herradura, ya sea la parte interna o externa, puede trabajar con velocidades relativas y la suposición de velocidades keplerianas.
• Se requiere precaución en cuanto a qué tiempo orbital le interesa: si$a$es la distancia de la Tierra al sol, y la herradura es$d$lejos de una órbita perfectamente circular, alternando así entre distancias$a\pm d$, entonces el tiempo orbital en el marco de reposo del sol será$v_{-}=\sqrt{\frac{GM}{a-d}}$, para un gran número de órbitas hasta el encuentro cercano, y después de eso$v_{+}=\sqrt{\frac{GM}{a+d}}$.
• Por lo tanto, el caso simplificado de que el satélite está en su propia órbita kepleriana es cierto la mayor parte del tiempo. Armados con este conocimiento, podemos aproximar el tiempo de recurrencia planetocéntrica como$t_{\rm rec}=\frac{2\pi a}{v_{K} - v_{-}}$como un simple tiempo de recuperación entre objetos en diferentes órbitas.
• La velocidad relativa$v_{\rm rel} = v_{K} - v_{-} = \sqrt{\frac{GM}{a}}-\sqrt{\frac{GM}{a-d}}$puede expandirse en el límite de$d/a\ll1$en$v_{\rm rel} \approx \frac{1}{2} \sqrt\frac{GM}{a} \frac{d}{a}$y así derivo de esto$t_{rec}\approx \frac{4\pi}{\sqrt{GM}} \frac{a^{5/2}}{d}$. Como era de esperar, el tiempo de ciclo diverge para$d=0$, como en el caso de la coorbita con la Tierra este tiempo debe ser infinito. Un límite superior en$d$no se puede derivar de esto, uno tendría que buscar la solución completa para eso.

Solo por curiosidad, ingresé algunos valores en esta fórmula y escribí algo rápido en python:

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

y los resultados que obtengo son

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

Por lo tanto, la fórmula puede tener un error de hasta un factor de ~2. Esto es seguro simplemente porque la realidad es más compleja que una simple aproximación de órbita circular, pero también debido a la calidad de los valores utilizados. Wikipedia no es muy conocida por investigar bien valores particulares. Tomé los que encontré allí. Para SO16, eso fue particularmente confuso de seleccionar, así que tomé los dos que se mencionaron en la misma línea de texto, con la esperanza de que vinieran de la misma fuente.

Cualquiera que encuentre valores más consistentes, es libre de comentar.

Ernest W. Brown's On a New Family of Periodic Orbits in the Problem of Three Bodies: (Plates 6, 7.) en MNRAS, 71, (5), pp 438–454 publicado el 10 de marzo de 1911 parece ser donde las órbitas de herradura estaban primero propuesto. (Disponible aquí también). Comienza:

Hay cuatro asteroides conocidos que parecen oscilar alrededor de uno u otro de los vértices de los dos triángulos equiláteros que tienen como base la línea que une a Júpiter y al Sol. Estos vértices son las conocidas posiciones de equilibrio relativo. El vector heliocéntrico de uno de estos asteroides aparentemente puede alejarse hasta 17° de su posición de equilibrio.* Por lo tanto, las oscilaciones no pueden considerarse muy pequeñas. Naturalmente, uno se pregunta si son posibles oscilaciones de este tipo en arcos de extensión aún mayor; y si es así, de qué manera se pueden obtener las órbitas más convenientemente.

*LJ Linders, Arhivfor Mat., Ast. och Fys., So. Veterinario. Alaska. yo Estocolmo, Bd. 4, núm. 20.

---

---

Haré algunas órbitas de herradura en el formalismo del problema circular restringido de tres cuerpos y las trazaré en Python, luego las compararé con la estimación del período sinódico descrita en la respuesta de @AtmosphericPrisonEscape .

tl;dr: Hay un buen acuerdo cualitativo, ¡sin sorpresas!

---

Un breve resumen de las matemáticas CR3BP en unidades adimensionales. La distancia entre los dos cuerpos es igual a 1, al igual que la constante gravitatoria. Orbitan alrededor de un centro de masa común en órbitas circulares, con un período de$2 \pi$. Es más fácil de visualizar y calcular si lo haces en un marco giratorio, por lo que las dos masas están fijas. El tercer cuerpo en la posición$x, y, z$se considera que no tiene efecto gravitatorio sobre los dos primeros,

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

La Energía de Jacobi$C$es una cantidad conservada en este marco giratorio:

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

donde el$x^2 + y^2$es el pseudopotencial. Si establece los términos dependientes de la velocidad$(\dot{x}^2 + \dot{y}^2 + \dot{z}^2)$a cero, obtienes una superficie de velocidad cero , esa superficie que se pega en muchas/la mayoría de las preguntas sobre tres órbitas corporales. Estos gráficos no se aplican cuando un objeto se está moviendo, por lo que no puede superponer órbitas encima de ellos.

La aceleración que siente el tercer cuerpo en este marco giratorio tiene tanto la esperada$1/r^2$fuerzas y una pseudofuerza dependiente de la velocidad que no es real, pero explica el hecho de que el marco gira y no es inercial.

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

---

Aquí hay algunos cálculos. Elegí$\mu = 0.001$que está bastante cerca de la situación de Júpiter y el Sol. Elegí una matriz de puntos de partida en el punto opuesto de$m_2$a aproximadamente$x=-1$ pero eso no es lo que realmente hice. Lo que realmente hice fue elegir un montón de velocidades iniciales$-0.08 < \dot{y} < 0.08$y para cada uno calculé la posición en el$x$eje cerca$x=-1$donde la aceleración en el$x$dirección era cero.

Eso le da a las soluciones un poco de simetría inicial, pero las órbitas de halo son irregulares y onduladas y no siempre tan estables, por lo que este esfuerzo no es realmente necesario.

Propagué cada órbita hasta que volvió a la misma área y la detuve cuando cruzó el eje x, produciendo una familia de semiciclos.

Para abreviar, el método que se muestra en la respuesta de @AtmosphericPrisonEscape de estimar el tiempo del ciclo calculando el período sinódico en el marco inercial está bastante de acuerdo con estas órbitas de halo, ¡y eso no debería ser muy sorprendente!

arriba: semiciclos de algunas órbitas de herradura tambaleantes

arriba: tiempos hasta los primeros cruces del eje x de las mismas órbitas de herradura tambaleantes, utilizadas para calcular tiempos de medio ciclo.

arriba: tiempos de ciclo de este cálculo (puntos negros) frente al método de estimación del período sinódico (puntos rojos). Buen acuerdo cualitativo. También las velocidades iniciales y en cada punto inicial en x.

a continuación: secuencia de comandos de Python para estas parcelas.

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