<p>That depends on the accuracy you want to work with. To zeroth order, as outlined in Murray &amp; Dermott, "Solar System dynamics", Chapt 3., you can do the following:</p>

<ul>
<li>The zero-velocity contours that are plotted in your image will not be coinciding with particle orbits to infinite precision, but they're a good zeroth order approximation for objects with low eccentricity w.r.t the star (<span class="math-container">$v_{r}/v_{\theta} \ll 1$</span>)</li>
<li>A particle on those orbits is on a regular Keplerian orbit of radius <span class="math-container">$r_{\rm H}$</span>, out of reach of the gravitational influence of the planet at <span class="math-container">$r_{\rm P}$</span>. Thus, to get the largest "chunk" of traveltime of one part of the horseshoe orbit, either inner or outer part, you can work out with relative velocities and the assumption of keplerian velocities.</li>
<li>Caution is required as to which orbital time you're interested in: If <span class="math-container">$a$</span> is Earth's distance to the sun, and the horseshoe is <span class="math-container">$d$</span> away from a perfectly circular orbit, thus alternating between distances <span class="math-container">$a\pm d$</span>, then the orbital time in the rest-frame of the sun will be <span class="math-container">$v_{-}=\sqrt{\frac{GM}{a-d}}$</span>, for a large number of orbits until close encounter, and after that <span class="math-container">$v_{+}=\sqrt{\frac{GM}{a+d}}$</span>.</li>
<li>Therefore, the simplified case that the satellite is on its own Keplerian orbit is true for most of the time. Armed with this knowledge, we can approximate the Planetocentric recurrence time as <span class="math-container">$t_{\rm rec}=\frac{2\pi a}{v_{K} - v_{-}}$</span> as a simple catch-up time between objects on different orbits.</li>
<li>The relative velocity <span class="math-container">$v_{\rm rel} = v_{K} - v_{-} = \sqrt{\frac{GM}{a}}-\sqrt{\frac{GM}{a-d}}$</span> can be expanded in the limit of <span class="math-container">$d/a\ll1$</span> into <span class="math-container">$v_{\rm rel} \approx \frac{1}{2} \sqrt\frac{GM}{a} \frac{d}{a}$</span> and thus I derive from this <span class="math-container">$t_{rec}\approx \frac{4\pi}{\sqrt{GM}} \frac{a^{5/2}}{d}$</span>. As expected, the cycle time diverges for <span class="math-container">$d=0$</span>, as in the co-orbiting case with Earth this time must be infinite. An upper limit on <span class="math-container">$d$</span> cannot be derived from this, one would need to turn for the full solution for that.</li>
</ul>

<p>Just out of curiosity, I've plugged in a few values into this formula and written something quick in python:</p>

<pre><code>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&gt;&gt;m1&gt;&gt;m2
# Semimajor axis distance is from m0 to m1, radial distance is a(m0-&gt;m1)-a(m0-&gt;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")
</code></pre>

<p>and the results I get are</p>

<pre><code>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
</code></pre>

<p>So the formula may be off by up to a factor of ~2. This is for sure simply because reality is more complex than a simple circular-orbit approximation, but also due to the quality of the values used.
Wikipedia isn't well known for researching particular values well. I took those that I found there. For SO16 that was particularly confusing to select, so I took those two that were mentioned in the same line of text, hoping they would come from the same source.</p>

<p>Anyone finding more consistent values, is free to comment.</p>

CC BY-SA 4.0
2019-02-12T18:12:54.113Z
While I understand that accuracy is an issue, to some extent, I'm having trouble reconciling this equation with the examples I provided of Earth+Cruithne and Janus+Epimetheus.  I'd understand a deviation of a few decades on a 770 year cycle, or several months or a  year on an 8 year cycle.  But I get values of hundreds of thousands of years on one, and a couple of weeks on the other.  What am I missing here?
Harthag
CC BY-SA 4.0
2019-02-12T19:10:48.347Z
@Dalila: I'm not getting those values. Mine seem much more reasonable, but it's still an approximation. You have to be careful on what is the semimajor axis and the value d, which is the difference of semimajor axis and the 2nd massive body. The central object is the Sun for Cruithne and Saturn for JEpi. But I've added my code, so you can check for yourself, and of course dig into the literature if you wish.
AtmosphericPrisonEscape
CC BY-SA 4.0
2019-02-12T19:35:09.327Z
I think this is where I ran in to trouble: "You have to be careful on what is the semimajor axis and the value d, which is the difference of semimajor axis and the 2nd massive body."  Can you clarify that value?  did you mean it's the difference between the semi-major axes of the two co-orbital bodies?  I just want to make sure I understood that last comment correctly
Harthag
CC BY-SA 4.0
2019-02-12T20:05:04.860Z
@Dalila: If we denote the bodies by 0,1,2, and the masses are $m_2 \gg m_1\ gg m_0$, then d is the difference of semimajor axes $a_{01}-a_{02}$. Note that this shouldn't exclude Janus and Epimetheus, as they obey the same physics as a regular Horseshoe-orbit: Each moon individually orbits Saturn unperturbed, for most of the time. So I'm not sure what's wrong with the value I get for JEpi.
AtmosphericPrisonEscape
CC BY-SA 4.0
2019-02-12T20:11:12.400Z
Thanks for that clarification.  I understood the JEpi cycle to be 8 years, total, not 4.  4 years was the time between the two approaches in the cycle, so 11 vs 8 is actually a reasonable margin of error, at least for my purposes, on that calculation.  Did I misunderstand something else on that one?
Harthag
CC BY-SA 4.0
2019-02-12T20:12:48.683Z
@Dalila: I don't think you misunderstood that. I don't know what the actual cycle time is, I just quickly looked up the values that seemed correctly place on wiki ;) So probably my mistake if it's actually 8 years.
AtmosphericPrisonEscape
CC BY-SA 4.0
2019-02-12T07:56:43.810Z
2019-02-13T23:34:24.243Z
2019-02-13T23:34:24.243Z
AtmosphericPrisonEscape
<p>Ernest W. Brown's <a href="https://academic.oup.com/mnras/article/71/5/438/972133" rel="nofollow noreferrer">On a New Family of Periodic Orbits in the Problem of Three Bodies: (Plates 6, 7.)</a> in MNRAS, 71, (5), pp 438–454 <strong>published on 10 March 1911</strong> seems to be where horseshoe orbits were first proposed. (Available <a href="http://adsabs.harvard.edu/abs/1911MNRAS..71..438B" rel="nofollow noreferrer">here</a> as well). It begins:</p>
<blockquote>
<p>There are four known asteroids which appear to oscillate about one or other of the vertices of the two equilateral triangles which have the line joining Jupiter and the Sun as base. These vertices are the well-known positions of relative equilibrium. The heliocentric vector of one of these asteroids can apparently move as far as 17° away from its equilibrium position.* The oscillations cannot therefore be considered very small. One naturally asks whether oscillations of this kind in arcs of still greater extent are possible ; and if so, in what manner the orbits may he most conveniently obtained.</p>
<p>*L. J. Linders, Arhivfor Mat., Ast. och Fys., So. Vet. Ak. i Stockholm, Bd. 4, No. 20.</p>
</blockquote>
<hr />
<hr />
<p>I'll make some horseshoe orbits in the Circular Restricted Three-Body Problem formalism and plot them in Python, then compare to the synodic period estimation described in @AtmosphericPrisonEscape's
<a href="https://astronomy.stackexchange.com/a/29561/7982">answer</a>.</p>
<p><strong>tl;dr:</strong> There's good qualitative agreement, no surprises!</p>
<hr />
<p>A brief summary of CR3BP math in dimensionless units. The distance between the two bodies is equal to 1, as is the gravitational constant. They orbit around a common center of mass in circular orbits, with a period of <span class="math-container">$2 \pi$</span>. It's easier to visualize and calculate if you do it in a rotating frame, so the two masses are fixed. The third body at position <span class="math-container">$x, y, z$</span> is considered to have no gravitational effect on the first two,</p>
<p><span class="math-container">$$\mu = \frac{m_2}{m_1 + m_2}$$</span></p>
<p><span class="math-container">$$x_1 = -\mu $$</span>
<span class="math-container">$$x_2 = 1-\mu $$</span></p>
<p><span class="math-container">$$r_1 = \sqrt{(x-x1)^2 + y^2 + z^2}$$</span>
<span class="math-container">$$r_2 = \sqrt{(x-x2)^2 + y^2 + z^2}$$</span></p>
<p>The Jacobi Energy <span class="math-container">$C$</span> is a conserved quantity in this rotating frame:</p>
<p><span class="math-container">$$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)$$</span></p>
<p>where the <span class="math-container">$x^2 + y^2$</span> is the pseudopotential. If you set the velocity dependent terms <span class="math-container">$(\dot{x}^2 + \dot{y}^2 + \dot{z}^2)$</span> to zero, you get a <em>zero velocity surface</em>, that surface that's pasted into many/most questions about three body orbits. These plots do not apply when an object is moving, and so you can't superimpose orbits on top of them!</p>
<p><a href="https://i.sstatic.net/fDq7l.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/fDq7l.png" alt="random zero-velocity curve" /></a></p>
<p>The acceleration felt by the third body in this rotating frame has both the expected <span class="math-container">$1/r^2$</span> forces and a velocity-dependent <em>pseudoforce</em> which is not real, but accounts for the fact that the frame is rotating and not inertial.</p>
<p><span class="math-container">$$\ddot{x} = x + 2\dot{y} - \frac{(1-\mu)(x+\mu)}{r_1^3} - \frac{\mu(x-1+\mu)}{r_2^3}$$</span>
<span class="math-container">$$\ddot{y} = y - 2\dot{x} - \frac{(1-\mu)y}{r_1^3} - \frac{\mu y}{r_2^3}$$</span>
<span class="math-container">$$\ddot{z} =  -\frac{(1-\mu) z}{r_1^3} - \frac{\mu z}{r_2^3} $$</span></p>
<hr />
<p>Here are some calculations. I chose <span class="math-container">$\mu = 0.001$</span> which is pretty close to the situation Jupiter and the Sun. I chose an array of starting points at the opposite point from <span class="math-container">$m_2$</span> at about <span class="math-container">$x=-1$</span> <em>but that's not what I really did.</em> What I really did is choose a bunch of starting velocities <span class="math-container">$-0.08 &lt; \dot{y} &lt; 0.08$</span> and for each I calculated the position on the <span class="math-container">$x$</span> axis near <span class="math-container">$x=-1$</span> where the acceleration in the <span class="math-container">$x$</span> direction was zero.</p>
<p>That gives the solutions a tiny bit of starting symmetry, but halo orbits are bumpy and wiggly and not always so stable, so this effort isn't really necessary.</p>
<p>I propagated each orbit until it came back around to the same area and stopped it when it crossed the x-axis, producing a family of half-cycles.</p>
<p>To make a long story short, the method shown in @AtmosphericPrisonEscape's
<a href="https://astronomy.stackexchange.com/a/29561/7982">answer</a> of estimating the cycle time by calculating the synodic period in the inertial frame is in pretty good agreement with these halo orbits, and that shouldn't be very surprising!</p>
<p><a href="https://i.sstatic.net/OkuSe.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/OkuSe.png" alt="enter image description here" /></a></p>
<p><strong>above:</strong> half-cycles of some wobbly horseshoe orbits</p>
<p><a href="https://i.sstatic.net/rWKMP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/rWKMP.png" alt="enter image description here" /></a></p>
<p><strong>above:</strong> times to first x-axis crossings of the same wobbly horseshoe orbits, used to calculate half-cycle times.</p>
<p><a href="https://i.sstatic.net/cHG8s.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/cHG8s.png" alt="enter image description here" /></a></p>
<p><strong>above:</strong> cycle times from this calculation (black dots) versus from the synodic period estimation method (red dots). Good qualitative agreement. Also the starting y velocities at each starting point in x.</p>
<p><strong>below:</strong> Python script for these plots.</p>
<pre><code>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&gt;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] &lt; 0
    zc    = np.argmax(y[2:]*y[1:-1] &lt; 0) + 1
    if zc &gt; 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()
</code></pre>

CC BY-SA 4.0
2019-02-14T13:50:43.990Z
Hah, that's nice. Yesterday I was also working on my own three-body solver to test the formula, but you beat me to it.
AtmosphericPrisonEscape
CC BY-SA 4.0
2019-02-14T14:29:24.983Z
@AtmosphericPrisonEscape I've been dabbling in this for a very [long time](https://space.stackexchange.com/q/20590/12102) but have a long way to go.
uhoh
CC BY-SA 4.0
2019-02-14T08:42:00.963Z
2019-02-14T14:32:46.093Z
2020-06-17T09:47:40.653Z
uhoh
<p><strong>EDIT</strong> This was NOT a duplicate of <a href="https://astronomy.stackexchange.com/questions/29270/horseshoe-orbit-cycle-times">Horseshoe orbit cycle times</a>.</p>

<p>But <strong>that other question has been deleted, regardless</strong>. My original question asked for the answer to the equation(s), and provides many of the variables.  But since no one was apparently willing/able to answer it, this new question asks only for the form of the equation(s) itself and a list of variables involved, so that I can either provide more information on a future question, or try and solve the equations myself (though I doubt my math ability to do that, hence the original question, if I know what equations and variables are involved I can at least make an attempt myself).</p>

<p>What variables need to be known in order to calculate a horseshoe orbit cycle time?</p>

<p><img src="https://upload.wikimedia.org/wikipedia/commons/d/d6/Lagrange_Horseshoe_Orbit.jpg" alt=""></p>

<p>In other words: What is the equation, and what do the variables in that equation represent (speed? mass? semi major axes? etc.), to calculate the time it takes for an object to go from Point A on that image, through Points B, C, D, and E and back to A, along the light blue line in the image above (image also available <a href="https://upload.wikimedia.org/wikipedia/commons/d/d6/Lagrange_Horseshoe_Orbit.jpg" rel="nofollow noreferrer">here</a>)?</p>

<p>For example, Earth and Cruithne complete their cycle in 770 years, while Janus and Epimetheus complete theirs in 8 years.  I'm confident that semi-major axes are key factors in the equation, but I don't know what other variables are included, or how the variables relate to each other in the calculation.  </p>

<p>Also, I understand that these orbits are not stable, in the long term, and that my example of Cruithne is especially complex, as it will actually change orbit types periodically, and will likely be ejected from the solar system entirely, or impact the Sun or Jupiter, at some point.  But those complexities are for another time.  For this question, I just want to know the equation in its simplest form. </p>

CC BY-SA 4.0
2019-02-09T02:26:34.343Z
That's not really a plot of a horseshoe orbit, it's a equipotential line for a pseudopotential in a rotating frame for an object with *zero velocity* in that frame. If the object were in a real horseshoe orbit, those potential lines would change. I wonder if you can find an image of an actual horseshoe orbit?
uhoh
CC BY-SA 4.0
2019-02-09T02:28:56.950Z
See also [Horseshoe orbit cycle times](https://astronomy.stackexchange.com/q/29270/7982) which also incorrectly uses that image and currently has no answer. Also unanswered is [Horseshoe orbits and integration in C](https://astronomy.stackexchange.com/q/19476/7982)
uhoh
CC BY-SA 4.0
2019-02-09T02:31:54.290Z
Wait, [Horseshoe orbit cycle times](https://astronomy.stackexchange.com/q/29270/7982) is your question, which you've also posted in [Physics SE](https://physics.stackexchange.com/q/455742/83380). How is this question different than those?
uhoh
CC BY-SA 4.0
2019-02-11T14:44:37.123Z
@uhoh  My original question asks for the answer to the equation with many of the variables already provided.  This one asks only for the equation(s) to make sure I know what all I need to provide on my original question, or maybe take a crack at solving it myself (though I doubt my math ability to complete the calculation myself, since no one is willing/able to answer my original question, I figured this would be my best shot at getting closer to the answer myself).
Harthag
CC BY-SA 4.0
2019-02-11T15:06:17.423Z
@uhoh  Also, I understand that the image is not (exactly) a "horseshoe" orbit specifically, but it IS a perfect representation of exactly what I actually want to know.  I want to know how to calculate the time for the complete cycle illustrated, how long it takes for one of the orbiting objects to make two approaches (one from each side) to the other orbiting object, how long it takes to complete one of those cycles.
Harthag
CC BY-SA 4.0
2019-02-11T15:11:18.740Z
that makes sense. I think the issue of the cycle time of horseshoe orbits is really interesting. I'd like to propose a course of action. Why don't you *temporarily* delete your other question to address any concerns of duplication, and let me do a bit of a rewrite on this question. Then I'll vote to re-open it. I'm planing on visiting the library tomorrow and I can do a search on this and see what analytical models there have been published.
uhoh
CC BY-SA 4.0
2019-02-11T15:15:46.233Z
When you delete a post, it doesn't really ever get deleted. For the first few months it still shows up on your list of questions. Here's mine, and you can see the link at the bottom for recently deleted ones  https://i.sstatic.net/G3Wh9.png You can undelete them any time. If you copy the link to your question and save it somewhere, then you can access it any time in the future.
uhoh
CC BY-SA 4.0
2019-02-11T15:34:27.803Z
These are all unstable orbits...
peterh
CC BY-SA 4.0
2019-02-11T19:10:19.647Z
@uhoh  Deleted.  I'm not worried about losing that other question, I have notes to be able to ask it again, if needed.
Harthag
CC BY-SA 4.0
2019-02-11T19:15:17.643Z
@peterh  If/when I figure out how to do these calculations, I'll use that information in a way that doesn't require long term stability in geological timescales.  I'm aware of their inherent instability and have taken that in to consideration for my final goals.  But thanks for the concern and pointing it out.
Harthag
CC BY-SA 4.0
2019-02-13T06:34:30.370Z
@Dalila I'm almost done, will post later today. This turns out to be a *really interesting problem!*
uhoh
CC BY-SA 4.0
2019-02-14T08:46:48.973Z
@Dalila answer posted, please let me know if you have further questions.
uhoh
<orbit><natural-satellites><orbital-mechanics><celestial-mechanics>
What variables are needed to calculate simple horseshoe orbit times?
Harthag
# What variables are needed to calculate simple horseshoe orbit times?

(Asked by: username_0 on 2019-02-09)

**EDIT** This was NOT a duplicate of [Horseshoe orbit cycle times](https://astronomy.stackexchange.com/questions/29270/horseshoe-orbit-cycle-times).

But **that other question has been deleted, regardless**. My original question asked for the answer to the equation(s), and provides many of the variables. But since no one was apparently willing/able to answer it, this new question asks only for the form of the equation(s) itself and a list of variables involved, so that I can either provide more information on a future question, or try and solve the equations myself (though I doubt my math ability to do that, hence the original question, if I know what equations and variables are involved I can at least make an attempt myself).

What variables need to be known in order to calculate a horseshoe orbit cycle time?

![](https://upload.wikimedia.org/wikipedia/commons/d/d6/Lagrange_Horseshoe_Orbit.jpg)

In other words: What is the equation, and what do the variables in that equation represent (speed? mass? semi major axes? etc.), to calculate the time it takes for an object to go from Point A on that image, through Points B, C, D, and E and back to A, along the light blue line in the image above (image also available [here](https://upload.wikimedia.org/wikipedia/commons/d/d6/Lagrange_Horseshoe_Orbit.jpg))?

For example, Earth and Cruithne complete their cycle in 770 years, while Janus and Epimetheus complete theirs in 8 years. I'm confident that semi-major axes are key factors in the equation, but I don't know what other variables are included, or how the variables relate to each other in the calculation.

Also, I understand that these orbits are not stable, in the long term, and that my example of Cruithne is especially complex, as it will actually change orbit types periodically, and will likely be ejected from the solar system entirely, or impact the Sun or Jupiter, at some point. But those complexities are for another time. For this question, I just want to know the equation in its simplest form.

**Tags:** orbit, natural-satellites, orbital-mechanics, celestial-mechanics

> **Comments:**
> **username_1**: That's not really a plot of a horseshoe orbit, it's a equipotential line for a pseudopotential in a rotating frame for an object with *zero velocity* in that frame. If the object were in a real horseshoe orbit, those potential lines would change. I wonder if you can find an image of an actual horseshoe orbit?  
> **username_1**: See also [Horseshoe orbit cycle times](https://astronomy.stackexchange.com/q/29270/7982) which also incorrectly uses that image and currently has no answer. Also unanswered is [Horseshoe orbits and integration in C](https://astronomy.stackexchange.com/q/19476/7982)  
> **username_1**: Wait, [Horseshoe orbit cycle times](https://astronomy.stackexchange.com/q/29270/7982) is your question, which you've also posted in [Physics SE](https://physics.stackexchange.com/q/455742/83380). How is this question different than those?  
> **username_0**: @username_1  My original question asks for the answer to the equation with many of the variables already provided.  This one asks only for the equation(s) to make sure I know what all I need to provide on my original question, or maybe take a crack at solving it myself (though I doubt my math ability to complete the calculation myself, since no one is willing/able to answer my original question, I figured this would be my best shot at getting closer to the answer myself).  
> **username_0**: @username_1  Also, I understand that the image is not (exactly) a "horseshoe" orbit specifically, but it IS a perfect representation of exactly what I actually want to know.  I want to know how to calculate the time for the complete cycle illustrated, how long it takes for one of the orbiting objects to make two approaches (one from each side) to the other orbiting object, how long it takes to complete one of those cycles.  
> **username_1**: that makes sense. I think the issue of the cycle time of horseshoe orbits is really interesting. I'd like to propose a course of action. Why don't you *temporarily* delete your other question to address any concerns of duplication, and let me do a bit of a rewrite on this question. Then I'll vote to re-open it. I'm planing on visiting the library tomorrow and I can do a search on this and see what analytical models there have been published.  
> **username_1**: When you delete a post, it doesn't really ever get deleted. For the first few months it still shows up on your list of questions. Here's mine, and you can see the link at the bottom for recently deleted ones  https://i.sstatic.net/G3Wh9.png You can undelete them any time. If you copy the link to your question and save it somewhere, then you can access it any time in the future.  
> **username_2**: These are all unstable orbits...  
> **username_0**: @username_1  Deleted.  I'm not worried about losing that other question, I have notes to be able to ask it again, if needed.  
> **username_0**: @username_2  If/when I figure out how to do these calculations, I'll use that information in a way that doesn't require long term stability in geological timescales.  I'm aware of their inherent instability and have taken that in to consideration for my final goals.  But thanks for the concern and pointing it out.  
> **username_1**: @Dalila I'm almost done, will post later today. This turns out to be a *really interesting problem!*  
> **username_1**: @Dalila answer posted, please let me know if you have further questions.  

---

## 2 answers

### Answer by username_3 on 2019-02-12. Score: 4

That depends on the accuracy you want to work with. To zeroth order, as outlined in Murray & Dermott, "Solar System dynamics", Chapt 3., you can do the following:

The zero-velocity contours that are plotted in your image will not be coinciding with particle orbits to infinite precision, but they're a good zeroth order approximation for objects with low eccentricity w.r.t the star ($v_{r}/v_{\theta} \ll 1$)
A particle on those orbits is on a regular Keplerian orbit of radius $r_{\rm H}$, out of reach of the gravitational influence of the planet at $r_{\rm P}$. Thus, to get the largest "chunk" of traveltime of one part of the horseshoe orbit, either inner or outer part, you can work out with relative velocities and the assumption of keplerian velocities.
Caution is required as to which orbital time you're interested in: If $a$ is Earth's distance to the sun, and the horseshoe is $d$ away from a perfectly circular orbit, thus alternating between distances $a\pm d$, then the orbital time in the rest-frame of the sun will be $v_{-}=\sqrt{\frac{GM}{a-d}}$, for a large number of orbits until close encounter, and after that $v_{+}=\sqrt{\frac{GM}{a+d}}$.
Therefore, the simplified case that the satellite is on its own Keplerian orbit is true for most of the time. Armed with this knowledge, we can approximate the Planetocentric recurrence time as $t_{\rm rec}=\frac{2\pi a}{v_{K} - v_{-}}$ as a simple catch-up time between objects on different orbits.
The relative velocity $v_{\rm rel} = v_{K} - v_{-} = \sqrt{\frac{GM}{a}}-\sqrt{\frac{GM}{a-d}}$ can be expanded in the limit of $d/a\ll1$ into $v_{\rm rel} \approx \frac{1}{2} \sqrt\frac{GM}{a} \frac{d}{a}$ and thus I derive from this $t_{rec}\approx \frac{4\pi}{\sqrt{GM}} \frac{a^{5/2}}{d}$. As expected, the cycle time diverges for $d=0$, as in the co-orbiting case with Earth this time must be infinite. An upper limit on $d$ cannot be derived from this, one would need to turn for the full solution for that.

Just out of curiosity, I've plugged in a few values into this formula and written something quick in 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")

and the results I get are
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

So the formula may be off by up to a factor of ~2. This is for sure simply because reality is more complex than a simple circular-orbit approximation, but also due to the quality of the values used.
Wikipedia isn't well known for researching particular values well. I took those that I found there. For SO16 that was particularly confusing to select, so I took those two that were mentioned in the same line of text, hoping they would come from the same source.
Anyone finding more consistent values, is free to comment.

> **Comments:**
> **username_0**: While I understand that accuracy is an issue, to some extent, I'm having trouble reconciling this equation with the examples I provided of Earth+Cruithne and Janus+Epimetheus.  I'd understand a deviation of a few decades on a 770 year cycle, or several months or a  year on an 8 year cycle.  But I get values of hundreds of thousands of years on one, and a couple of weeks on the other.  What am I missing here?  
> **username_3**: @Dalila: I'm not getting those values. Mine seem much more reasonable, but it's still an approximation. You have to be careful on what is the semimajor axis and the value d, which is the difference of semimajor axis and the 2nd massive body. The central object is the Sun for Cruithne and Saturn for JEpi. But I've added my code, so you can check for yourself, and of course dig into the literature if you wish.  
> **username_0**: I think this is where I ran in to trouble: "You have to be careful on what is the semimajor axis and the value d, which is the difference of semimajor axis and the 2nd massive body."  Can you clarify that value?  did you mean it's the difference between the semi-major axes of the two co-orbital bodies?  I just want to make sure I understood that last comment correctly  
> **username_3**: @Dalila: If we denote the bodies by 0,1,2, and the masses are $m_2 \gg m_1\ gg m_0$, then d is the difference of semimajor axes $a_{01}-a_{02}$. Note that this shouldn't exclude Janus and Epimetheus, as they obey the same physics as a regular Horseshoe-orbit: Each moon individually orbits Saturn unperturbed, for most of the time. So I'm not sure what's wrong with the value I get for JEpi.  
> **username_0**: Thanks for that clarification.  I understood the JEpi cycle to be 8 years, total, not 4.  4 years was the time between the two approaches in the cycle, so 11 vs 8 is actually a reasonable margin of error, at least for my purposes, on that calculation.  Did I misunderstand something else on that one?  
> **username_3**: @Dalila: I don't think you misunderstood that. I don't know what the actual cycle time is, I just quickly looked up the values that seemed correctly place on wiki ;) So probably my mistake if it's actually 8 years.  

---

### Answer by username_1 on 2019-02-14. Score: 2

Ernest W. Brown's On a New Family of Periodic Orbits in the Problem of Three Bodies: (Plates 6, 7.) in MNRAS, 71, (5), pp 438–454 published on 10 March 1911 seems to be where horseshoe orbits were first proposed. (Available here as well). It begins:

There are four known asteroids which appear to oscillate about one or other of the vertices of the two equilateral triangles which have the line joining Jupiter and the Sun as base. These vertices are the well-known positions of relative equilibrium. The heliocentric vector of one of these asteroids can apparently move as far as 17° away from its equilibrium position.* The oscillations cannot therefore be considered very small. One naturally asks whether oscillations of this kind in arcs of still greater extent are possible ; and if so, in what manner the orbits may he most conveniently obtained.
*L. J. Linders, Arhivfor Mat., Ast. och Fys., So. Vet. Ak. i Stockholm, Bd. 4, No. 20.



I'll make some horseshoe orbits in the Circular Restricted Three-Body Problem formalism and plot them in Python, then compare to the synodic period estimation described in @username_3's
answer.
tl;dr: There's good qualitative agreement, no surprises!

A brief summary of CR3BP math in dimensionless units. The distance between the two bodies is equal to 1, as is the gravitational constant. They orbit around a common center of mass in circular orbits, with a period of $2 \pi$. It's easier to visualize and calculate if you do it in a rotating frame, so the two masses are fixed. The third body at position $x, y, z$ is considered to have no gravitational effect on the first two,
$$\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}$$
The Jacobi Energy $C$ is a conserved quantity in this rotating frame:
$$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)$$
where the $x^2 + y^2$ is the pseudopotential. If you set the velocity dependent terms $(\dot{x}^2 + \dot{y}^2 + \dot{z}^2)$ to zero, you get a zero velocity surface, that surface that's pasted into many/most questions about three body orbits. These plots do not apply when an object is moving, and so you can't superimpose orbits on top of them!

The acceleration felt by the third body in this rotating frame has both the expected $1/r^2$ forces and a velocity-dependent pseudoforce which is not real, but accounts for the fact that the frame is rotating and not inertial.
$$\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} $$

Here are some calculations. I chose $\mu = 0.001$ which is pretty close to the situation Jupiter and the Sun. I chose an array of starting points at the opposite point from $m_2$ at about $x=-1$ but that's not what I really did. What I really did is choose a bunch of starting velocities $-0.08 < \dot{y} < 0.08$ and for each I calculated the position on the $x$ axis near $x=-1$ where the acceleration in the $x$ direction was zero.
That gives the solutions a tiny bit of starting symmetry, but halo orbits are bumpy and wiggly and not always so stable, so this effort isn't really necessary.
I propagated each orbit until it came back around to the same area and stopped it when it crossed the x-axis, producing a family of half-cycles.
To make a long story short, the method shown in @username_3's
answer of estimating the cycle time by calculating the synodic period in the inertial frame is in pretty good agreement with these halo orbits, and that shouldn't be very surprising!

above: half-cycles of some wobbly horseshoe orbits

above: times to first x-axis crossings of the same wobbly horseshoe orbits, used to calculate half-cycle times.

above: cycle times from this calculation (black dots) versus from the synodic period estimation method (red dots). Good qualitative agreement. Also the starting y velocities at each starting point in x.
below: Python script for these plots.
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()

> **Comments:**
> **username_3**: Hah, that's nice. Yesterday I was also working on my own three-body solver to test the formula, but you beat me to it.  
> **username_1**: @username_3 I've been dabbling in this for a very [long time](https://space.stackexchange.com/q/20590/12102) but have a long way to go.  

---

