How to know when a body completes an orbit?

Imagine a vector from the central body to the orbiting body. An orbit has been completed when this vector points in the same direction as when the simulation started. The orbital period is the time it takes for the objects to return to their original orientation with respect to each other.

In programming terms, let's say the planet has a position $vec{x}_p$ and the star has a position $vec{x}_s$. Define a vector $vec{r} = vec{x}_p - vec{x}_s$. The angle of $vec{r}$ is given by $theta = textrm{atan2}(r_y, r_x)$. An orbit has been completed when $theta$ returns to its original value.

Depending on what information you start with, another way to predict the period is to use Kepler's third law: $$T = sqrt{frac{4pi^2 a^3}{G(M+m)}}$$ where $T$ is the orbital period, $a$ is the semi-major axis of the orbit (half the larger width of the ellipse of the orbit, or the average of the closest and farthest distances between the bodies), $G$ is the gravitational constant, and $M$ and $m$ are the masses of the orbiting bodies.

If your starting point is just the position and velocity of the orbiting body, then you'll need the more complicated math found in Brian Opatosky's answer.

For the two-body central force problem, the radial equation is:

$$r(theta)=frac{p}{1+e,cos(theta - phi)}$$

where $$phi$$ is the phase angle, $$theta$$ is the angle of the satellite object in its orbit, $$p, equiv , frac{L^2}{Gm^2M}$$

(L is the angular momentum, m is the mass of the satellite, M is the mass of the central body, and G is Newton's gravitational constant), and $$e$$ is the eccentricity for the orbit. For an elliptical orbit:

$$ e, =, sqrt{frac{2,E_0,m,p^2}{{L_0}^2},+,1}$$

Note that the angular momentum is:

$$overrightarrow{L_{net}}=overrightarrow{r} times overrightarrow{p}=$$

$$r,hat{r} times m, (dot r hat{r} + r,dottheta hat{theta})=$$

$$m,r^2,dottheta,hat{z}$$

Note that, for strictly central force problems:

$$dot{overrightarrow{L_{net}}},=,tau_{net},=, overrightarrow{r} timessum_i, overrightarrow{F_i},=, r,hat{r}, times, sum_i, F_i, hat{r}, =, overrightarrow{0} $$

This means that, for a central force problem, the angular momentum is constant. Thus:

$$dottheta, =, frac{L_0}{m,r^2}$$

where $$L_0$$ is that inital angular momentum (or any convenient time).

Thus:

$$frac{dtheta}{dt},=, frac{L_0}{m,r^2}$$

So:

$$frac{m,r^2, dtheta}{L_0}, =, dt$$

$$frac{m}{L_0}int_{theta_0}^{theta_0,+,2pi},r^2,dtheta,=,int_{0}^{T}, dt$$

where T is the period and the limits of the integral on the left were chosen from any arbitrary angle and the upper limit was is because of the $$2pi$$ periodicity of the radial function (at the top).

Plugging in the radial equation:

$$frac{m}{L_0}int_{theta_0}^{theta_0,+,2pi},({frac{p}{1,+,e,cos(theta-phi)}})^2,dtheta,=,int_{0}^{T}, dt$$

$$frac{m,p^2}{L_0}int_{theta_0}^{theta_0,+,2pi},[1,+,e,cos(theta-phi)]^{-2},dtheta,=,T$$

Using a change of variables:

$$u ,equiv, (theta, -, phi),,, Rightarrow,,,du,=, dtheta, $$and$$, u_{initial},=, (theta_0,-,phi) , $$and$$u_{final},=,(theta_0,+,2pi,-,phi)$$

The equation then becomes:

$$frac{m,p^2}{L_0}int_{(theta_0,-,phi)}^{(theta_0,+,2pi,-,phi)},[1,+,e,cos(u)]^{-2},du,=,T$$

That integral is still pretty nasty... here is the result without the going through it:

$$left.left(frac{e,sin(u)}{(e^2 -1),(e,cos(u),+,1)},-,frac{2,arctanhleft(frac{(e-1),tan(frac{u}{2})}{sqrt{e^2 -1}}right)}{(e^2 -1)^{frac{3}{2}}}right)rightvert_{theta_0,-,phi}^{theta_0,+,2pi,-,phi},=,,,T$$

While this solution is exact... there's a problem. The "function":

$$arctanh(z)$$

is not truly a function, because it is not single-valued. It is a multivalued analytic function. Basically, while the hyperbolic tangent function is a true function, its inverse is not, as single inputs have multiple outputs. To overcome this, something called branch cuts are made and principal values are defined. Put simply, the domain and range of arctanh is restricted so that it behaves like a function. In doing this, a variety of values are neglected. So, in this particular problem, plugging in the arbitrary bounds always yields:

$$T, =, 0$$

which is the solution of how long it takes to reach the same point... without moving, which is useless. The more useful period is lost because of the branch cutting. All hope isn't lost though! Instead of using the exact, analytic solution, we can just define the previous integral function as our "period function". The true, useful period can be computed via numerical integration (would highly recommend doing it computationally). I ran a few tests and found Simpson's Rule integration and the Trapezoid method to be highly accurate, though Monte Carlo integration was fairly accurate as well. All that is left is to recognize that the initial angle $$theta_0$$ and phase angle $$phi$$ are arbitrary when plugged into the solution, as they will cancel out. Thus, written simpler:

$$T,=,frac{m,p^2}{L_0}int_{0}^{2pi},[1,+,e,cos(u)]^{-2},du$$

I hope this helped! Good luck with the project; it sounds awesome.