<p>I won't go through your equations carefully nor your code, but since you haven't mentioned the exact starting positions that you used, my guess is that you haven't first <em>solved for a stable horseshoe orbit</em> state vector first.</p>

<p>The positions for the Sun and Earth are fixed in the synodic (rotating) frame. But where did you start your astroid? For a given starting position, you have to also give the asteroid the correct starting velocity (speed and direction) otherwise it may not be stable and just wander off or spiral around as you mentioned.</p>

<p>Have a look at <a href="https://astronomy.stackexchange.com/a/29596/7982">this answer</a> where I show the correct equations and then solve for stable horseshoe orbits before plotting them.</p>

<p>I start with a position on the x-axis, opposite the Earth, and a little farther or closer from the Sun than the Earth is. Then I launch it with a low velocity (in the rotating frame) exactly in the direction of the y-axis. I watch it drift towards Earth and then boomerang back to the region of the starting area. When it crosses the x-axis again, I check to see how close the velocity is to y. I use an optimization zero solver <code>brentq</code> to adjust the initial velocity until the return velocity is in exactly the opposite direction.</p>

<p>If that's true, then the orbit will be periodic and repeating, and we can call it a stale orbit under the constraints of the CR3BP.</p>

<hr>

<p>From <a href="https://astronomy.stackexchange.com/a/29596/7982">this answer</a> (there's a lot more to read there, including several references!):</p>

<blockquote>
  <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>
  
  <p><img src="https://i.sstatic.net/OkuSe.png" alt="enter image description here"></p>
  
  <p><strong>above:</strong> half-cycles of some wobbly horseshoe orbits</p>
  
  <p><img src="https://i.sstatic.net/rWKMP.png" alt="enter image description here"></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><img src="https://i.sstatic.net/cHG8s.png" alt="enter image description here"></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>
</blockquote>

<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-07-15T07:21:45.547Z
2019-07-15T07:31:23.930Z
2019-07-15T07:31:23.930Z
uhoh
<p>I'm studying a particular case of the restricted three-body problem. It's been found that some objects follow a horseshoe orbit pattern, and I'm trying to sort out something through an integration code in C. I'm following some advice in the article <a href="http://www.aanda.org/articles/aa/abs/2005/12/aa1483/aa1483.html" rel="noreferrer">Families of periodic horseshoe orbits in the restricted three-body problem</a>, which gives me ideal initial conditions and the equations in the centre of mass system. (m is the mass of Earth, and consequent position of the sun in the center of mass reference system, (x,y) are the coordinates of the third body, assumed massless (as the restricted problem requires).</p>

<p>$$ O=(x^2+y^2)/2+\dfrac{(1-m)}{r1}+\dfrac{m}{r2}+\dfrac{(1-m)m}{2} $$
$$r1^2=(x-m)^2+y^2$$<br>
$$r2^2=(x-m+1)^2+y^2$$
$$ a(x)=\dfrac{\mathrm{d}{O}}{\mathrm{d}x} + 2v(y)$$
$$a(y)=\dfrac{\mathrm{d}{O}}{\mathrm{d}y}-2v(x)$$</p>

<p>The positions of the "sun" and "earth" are fixed to (m,0) and (m-1,0), in the same reference system. (rotating reference system, assuming the earth has a circular orbit.)</p>

<p>From all this I've calculated the equations to describe the system:</p>

<p>$$a(x)=x+\dfrac{(m-1)(x-m)}{((x-m)^2+y^2)^1.5}-\dfrac{2m(x-m+1)}{((x-m+1)^2+y^2)^1.5} + 2v(y)$$
$$a(y)=y-\dfrac{y(1-m)}{((x-m)^2+y^2)^1.5} -\dfrac{2ym}{((x-m+1)^2+y^2)^1.5} -2v(x)$$</p>

<p>I've used the algorithm of Runge-Kutta 4 to integrate those equations.
(I know the code is quite mind-twisting, but I just can't use pointers and I use structs everywhere).</p>

<pre><code>#include&lt;stdio.h&gt;
#include&lt;stdlib.h&gt;
#include&lt;math.h&gt;

#define dt 0.0001
#define N 100

typedef struct{
    long double x,y;
}vec;

typedef struct{
    vec k1,k2,k3,k4;
}runge;

typedef struct{
    runge r,v;
}big;

double dS,dE,m;

double accx(double,double,double);
double accy(double,double,double);
void rad(vec);
big rungekutta(vec,vec);
vec moto(vec,runge);
double jacobi(vec);

int main(){
  vec r,v;
    big f;
    double J,t;
    int i,Num;
    FILE* s1;
    s1=fopen("HorseShoe.dat","w");

    Num=(int)N/dt;
    scanf("%Lf",&amp;r.x);
    scanf("%Lf",&amp;r.y);
    scanf("%Lf",&amp;v.x);
    scanf("%Lf",&amp;v.y);
    scanf("%lf",&amp;m);

    for(i=0;i&lt;Num;i++){
        t=(i+1)*dt;
        rad(r);
        f=rungekutta(r,v);
        r=moto(r,f.r);
        v=moto(v,f.v);
        J=jacobi(r);
        fprintf(s1,"%lf\t%Lf\t%Lf\t%Lf\t%Lf\t%lf\n",t,r.x,r.y,v.x,v.y,J);
    }
return 0;
}

void rad(vec r){
    dS=pow(r.x-m,2)+pow(r.y,2);
    dE=pow(r.x-m+1,2)+pow(r.y,2);
}

double jacobi(vec r){
    return pow(r.x,2)+pow(r.y,2)+2*(1-m)/dS+2*m/dE+m*(1-m);
}

double accx(double x,double y,double v){
    return x-(x-m)*(1-m)/pow(pow(x-m,2)+pow(y,2),1.5)-m*(x-m+1)/pow(pow(x-m+1,2)+pow(y,2),1.5)+2*v;
}

double accy(double x,double y,double v){
    return y-(1-m)*y/pow(pow(y,2)+pow(x-m,2),1.5)-m*y/pow(pow(y,2)+pow(x-m+1,2),1.5)-2*v;
}

big rungekutta(vec r,vec v){
    big f;
    f.r.k1.x=v.x;
    f.r.k1.y=v.y;
    f.v.k1.x=accx(r.x,r.y,v.y);
    f.v.k1.y=accy(r.x,r.y,v.x);
    f.r.k2.x=v.x+f.v.k1.x*dt/2;
    f.r.k2.y=v.y+f.v.k1.y*dt/2;
    f.v.k2.x=accx(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.y+f.v.k1.y*dt/2);
    f.v.k2.y=accy(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.x+f.v.k1.x*dt/2);
    f.r.k3.x=v.x+f.v.k2.x*dt/2;
    f.r.k3.y=v.y+f.v.k2.y*dt/2;
    f.v.k3.x=accx(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.y+f.v.k2.y*dt/2);
    f.v.k3.y=accy(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.x+f.v.k2.x*dt/2);
    f.r.k4.x=v.x+f.v.k3.x*dt;
    f.r.k4.y=v.y+f.v.k3.y*dt;
    f.v.k4.x=accx(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.y+f.v.k3.y*dt);
    f.v.k4.y=accy(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.x+f.v.k3.x*dt);
    return f;
}

vec moto(vec r,runge rk){
    r.x+=(rk.k1.x+2*rk.k2.x+2*rk.k3.x+rk.k4.x)*dt/6;
    r.y+=(rk.k1.y+2*rk.k2.y+2*rk.k3.y+rk.k4.y)*dt/6;
    return r;
}
</code></pre>

<p>Plotting the results I just get a spiral, while using the given inputs I should get a horseshoe orbit. I've tried many different inputs, (m=0.0001 and m=0.000003, the latter is in scale with the actual values of Earth and Sun masses (sun mass is 1-m)).</p>

<p>I just can't see what's wrong (probably everything :D), can anyone please help me?</p>

CC BY-SA 3.0
2016-12-21T22:30:39.740Z
People might be better able to help if you explain what your variables (in the equations) are: O,x,y,m,r1/2 etc. Also, your question doesn't make it clear what the problem is that you are encountering. "I can't see what's wrong" -- what are you supposed to get, and what do you get instead?
Alex
CC BY-SA 3.0
2016-12-22T03:55:38.223Z
Before proceeding, I would like to see two things: A link to the reference paper you used so as to see if you made a mistake in your math, and your inputs. I suspect your inputs might be off. This is a non-standard system of units in which the mass of the Earth $m$ should be about $3\times10^{-6}$, the mass of the Sun is 1 less this small number, and the magnitude of the position vector $(x,y)$ is about 1. Time appears to be scaled, also. This is not a place for SI units.
David Hammen
CC BY-SA 3.0
2016-12-22T08:31:38.577Z
I've just edited the question adding the details you asked (sorry if I'm a bit messy, first question in here). I used G=1, of course I couldn't use SI units ;)
Elisa
CC BY-SA 3.0
2016-12-22T08:34:40.450Z
I used many of the inputs given in the articles (at least one per family of orbits), I give you an example. (x0=0.864394016091,y0=0,vx=0,vy=0.288028401448,m=0.0001)
Elisa
CC BY-SA 3.0
2016-12-22T08:57:55.447Z
Wow, I've just found the mistake, I feel so stupid now.. Thank you anyways ;)
Elisa
CC BY-SA 3.0
2016-12-22T19:22:25.987Z
Feel free to post an answer to your own question detailing what was wrong and how you fixed it. This may prove useful to others in the future who come to this question!
zephyr
CC BY-SA 3.0
2016-12-23T10:31:47.133Z
Just a small but important remark: 4ᵗʰ order Runge-Kutta is not [symplectic](https://en.wikipedia.org/wiki/Symplectic_integrator), meaning, it changes the total orbital energy over time, slow at first but increasing exponentially in its severity. Depending on your desired integration interval, you had better ditch RK4 in favor of something more suited for celestial mechanics. See for example [this paper](https://arxiv.org/abs/1511.07476) for some background and recommendations.
Rody Oldenhuis
<gravity><n-body-simulations>
Horseshoe orbits and integration in C
Elisa
# Horseshoe orbits and integration in C

(Asked by: username_0 on 2016-12-21)

I'm studying a particular case of the restricted three-body problem. It's been found that some objects follow a horseshoe orbit pattern, and I'm trying to sort out something through an integration code in C. I'm following some advice in the article [Families of periodic horseshoe orbits in the restricted three-body problem](http://www.aanda.org/articles/aa/abs/2005/12/aa1483/aa1483.html), which gives me ideal initial conditions and the equations in the centre of mass system. (m is the mass of Earth, and consequent position of the sun in the center of mass reference system, (x,y) are the coordinates of the third body, assumed massless (as the restricted problem requires).

$$ O=(x^2+y^2)/2+\dfrac{(1-m)}{r1}+\dfrac{m}{r2}+\dfrac{(1-m)m}{2} $$
$$r1^2=(x-m)^2+y^2$$  
$$r2^2=(x-m+1)^2+y^2$$
$$ a(x)=\dfrac{\mathrm{d}{O}}{\mathrm{d}x} + 2v(y)$$
$$a(y)=\dfrac{\mathrm{d}{O}}{\mathrm{d}y}-2v(x)$$

The positions of the "sun" and "earth" are fixed to (m,0) and (m-1,0), in the same reference system. (rotating reference system, assuming the earth has a circular orbit.)

From all this I've calculated the equations to describe the system:

$$a(x)=x+\dfrac{(m-1)(x-m)}{((x-m)^2+y^2)^1.5}-\dfrac{2m(x-m+1)}{((x-m+1)^2+y^2)^1.5} + 2v(y)$$
$$a(y)=y-\dfrac{y(1-m)}{((x-m)^2+y^2)^1.5} -\dfrac{2ym}{((x-m+1)^2+y^2)^1.5} -2v(x)$$

I've used the algorithm of Runge-Kutta 4 to integrate those equations.
(I know the code is quite mind-twisting, but I just can't use pointers and I use structs everywhere).

```
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define dt 0.0001
#define N 100

typedef struct{
    long double x,y;
}vec;

typedef struct{
    vec k1,k2,k3,k4;
}runge;

typedef struct{
    runge r,v;
}big;

double dS,dE,m;

double accx(double,double,double);
double accy(double,double,double);
void rad(vec);
big rungekutta(vec,vec);
vec moto(vec,runge);
double jacobi(vec);

int main(){
  vec r,v;
    big f;
    double J,t;
    int i,Num;
    FILE* s1;
    s1=fopen("HorseShoe.dat","w");

    Num=(int)N/dt;
    scanf("%Lf",&r.x);
    scanf("%Lf",&r.y);
    scanf("%Lf",&v.x);
    scanf("%Lf",&v.y);
    scanf("%lf",&m);

    for(i=0;i<Num;i++){
        t=(i+1)*dt;
        rad(r);
        f=rungekutta(r,v);
        r=moto(r,f.r);
        v=moto(v,f.v);
        J=jacobi(r);
        fprintf(s1,"%lf\t%Lf\t%Lf\t%Lf\t%Lf\t%lf\n",t,r.x,r.y,v.x,v.y,J);
    }
return 0;
}

void rad(vec r){
    dS=pow(r.x-m,2)+pow(r.y,2);
    dE=pow(r.x-m+1,2)+pow(r.y,2);
}

double jacobi(vec r){
    return pow(r.x,2)+pow(r.y,2)+2*(1-m)/dS+2*m/dE+m*(1-m);
}

double accx(double x,double y,double v){
    return x-(x-m)*(1-m)/pow(pow(x-m,2)+pow(y,2),1.5)-m*(x-m+1)/pow(pow(x-m+1,2)+pow(y,2),1.5)+2*v;
}

double accy(double x,double y,double v){
    return y-(1-m)*y/pow(pow(y,2)+pow(x-m,2),1.5)-m*y/pow(pow(y,2)+pow(x-m+1,2),1.5)-2*v;
}

big rungekutta(vec r,vec v){
    big f;
    f.r.k1.x=v.x;
    f.r.k1.y=v.y;
    f.v.k1.x=accx(r.x,r.y,v.y);
    f.v.k1.y=accy(r.x,r.y,v.x);
    f.r.k2.x=v.x+f.v.k1.x*dt/2;
    f.r.k2.y=v.y+f.v.k1.y*dt/2;
    f.v.k2.x=accx(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.y+f.v.k1.y*dt/2);
    f.v.k2.y=accy(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.x+f.v.k1.x*dt/2);
    f.r.k3.x=v.x+f.v.k2.x*dt/2;
    f.r.k3.y=v.y+f.v.k2.y*dt/2;
    f.v.k3.x=accx(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.y+f.v.k2.y*dt/2);
    f.v.k3.y=accy(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.x+f.v.k2.x*dt/2);
    f.r.k4.x=v.x+f.v.k3.x*dt;
    f.r.k4.y=v.y+f.v.k3.y*dt;
    f.v.k4.x=accx(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.y+f.v.k3.y*dt);
    f.v.k4.y=accy(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.x+f.v.k3.x*dt);
    return f;
}

vec moto(vec r,runge rk){
    r.x+=(rk.k1.x+2*rk.k2.x+2*rk.k3.x+rk.k4.x)*dt/6;
    r.y+=(rk.k1.y+2*rk.k2.y+2*rk.k3.y+rk.k4.y)*dt/6;
    return r;
}

```

Plotting the results I just get a spiral, while using the given inputs I should get a horseshoe orbit. I've tried many different inputs, (m=0.0001 and m=0.000003, the latter is in scale with the actual values of Earth and Sun masses (sun mass is 1-m)).

I just can't see what's wrong (probably everything :D), can anyone please help me?

**Tags:** gravity, n-body-simulations

> **Comments:**
> **username_1**: People might be better able to help if you explain what your variables (in the equations) are: O,x,y,m,r1/2 etc. Also, your question doesn't make it clear what the problem is that you are encountering. "I can't see what's wrong" -- what are you supposed to get, and what do you get instead?  
> **username_2**: Before proceeding, I would like to see two things: A link to the reference paper you used so as to see if you made a mistake in your math, and your inputs. I suspect your inputs might be off. This is a non-standard system of units in which the mass of the Earth $m$ should be about $3\times10^{-6}$, the mass of the Sun is 1 less this small number, and the magnitude of the position vector $(x,y)$ is about 1. Time appears to be scaled, also. This is not a place for SI units.  
> **username_0**: I've just edited the question adding the details you asked (sorry if I'm a bit messy, first question in here). I used G=1, of course I couldn't use SI units ;)  
> **username_0**: I used many of the inputs given in the articles (at least one per family of orbits), I give you an example. (x0=0.864394016091,y0=0,vx=0,vy=0.288028401448,m=0.0001)  
> **username_0**: Wow, I've just found the mistake, I feel so stupid now.. Thank you anyways ;)  
> **username_3**: Feel free to post an answer to your own question detailing what was wrong and how you fixed it. This may prove useful to others in the future who come to this question!  
> **username_4**: Just a small but important remark: 4ᵗʰ order Runge-Kutta is not [symplectic](https://en.wikipedia.org/wiki/Symplectic_integrator), meaning, it changes the total orbital energy over time, slow at first but increasing exponentially in its severity. Depending on your desired integration interval, you had better ditch RK4 in favor of something more suited for celestial mechanics. See for example [this paper](https://arxiv.org/abs/1511.07476) for some background and recommendations.  

---

## 1 answers

### Answer by username_5 on 2019-07-15. Score: 1

I won't go through your equations carefully nor your code, but since you haven't mentioned the exact starting positions that you used, my guess is that you haven't first *solved for a stable horseshoe orbit* state vector first.

The positions for the Sun and Earth are fixed in the synodic (rotating) frame. But where did you start your astroid? For a given starting position, you have to also give the asteroid the correct starting velocity (speed and direction) otherwise it may not be stable and just wander off or spiral around as you mentioned.

Have a look at [this answer](https://astronomy.stackexchange.com/a/29596/7982) where I show the correct equations and then solve for stable horseshoe orbits before plotting them.

I start with a position on the x-axis, opposite the Earth, and a little farther or closer from the Sun than the Earth is. Then I launch it with a low velocity (in the rotating frame) exactly in the direction of the y-axis. I watch it drift towards Earth and then boomerang back to the region of the starting area. When it crosses the x-axis again, I check to see how close the velocity is to y. I use an optimization zero solver `brentq` to adjust the initial velocity until the return velocity is in exactly the opposite direction.

If that's true, then the orbit will be periodic and repeating, and we can call it a stale orbit under the constraints of the CR3BP.

---

From [this answer](https://astronomy.stackexchange.com/a/29596/7982) (there's a lot more to read there, including several references!):

> $$\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} $$
>
> ![enter image description here](https://i.sstatic.net/OkuSe.png)
>
> **above:** half-cycles of some wobbly horseshoe orbits
>
> ![enter image description here](https://i.sstatic.net/rWKMP.png)
>
> **above:** times to first x-axis crossings of the same wobbly horseshoe orbits, used to calculate half-cycle times.
>
> ![enter image description here](https://i.sstatic.net/cHG8s.png)
>
> **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()

```


---

