Equations

The equations of motion

Each spacecraft dynamics is described, between two manoeuvres, by the following set of Ordinary Differential Equations (ODEs):

$$ \begin{array}{l} \ddot {\mathbf x} = - \frac{\mu x}{r^3} \left\{ 1 + \frac 32 J_2 \left( \frac {r_{eq}}r \right)^2 \left( 1 - 5 \frac {z^2}{r^2} \right) \right\}\\ \ddot {\mathbf y} = - \frac{\mu y}{r^3} \left\{ 1 + \frac 32 J_2 \left( \frac {r_{eq}}r \right)^2 \left( 1- 5 \frac {z^2}{r^2} \right) \right\}\\ \ddot {\mathbf z} = - \frac{\mu z}{r^3} \left\{ 1 + \frac 32 J_2 \left( \frac {r_{eq}}r \right)^2 \left( 3 - 5 \frac {z^2}{r^2} \right) \right\}\\ \end{array} $$ that describe a Keplerian motion perturbed by main effects of an oblate Earth, i.e. \(J_2\). Note that between an arrival and a departure event the spacecarft is co-orbiting with the debris piece and hence its position and velocity is not described by the above equations, rather it coincides with the orbiting debris'.

The orbital ephemerides

Each debris orbit is defined by the values \(t_0, a, e, i, \Omega_0, \omega_0, M_0\) as read from the file debris.csv part of the competition data. At any given epoch \(t\) the position and velocity vectors of each debris piece must be computed by updating its osculating Keplerian elements using the mean motion and the precession rates and then converting, as in the Keplerian case, the updated osculating elements to position and velocity. Note that by doing so we are neglecting the velocity component deriving from \(\dot \Omega\) and \(\dot \omega\), for the purpose of this competition this is deemed as appropriate as it removes complexity from the equations without introducing any significant change on the search space landscape.

The procedure detailed below (assumes consistent units everywhere) shows all necessary equations.

1 - Compute the osculating Keplerian parameters

After having defined the mean motion \(n = \sqrt{\frac{\mu}{a^3}}\), the semilatus rectum \(p = a(1-e^2)\) and the precession rates: $$ \dot \Omega =-\frac 32 J_2 \left( \frac{r_{eq}} {p} \right)^2 n \cos i $$ $$ \dot \omega = \frac 34 J_2 \left( \frac{r_{eq}} {p} \right)^2 n (5\cos^2 i - 1) $$ compute the right ascension of the ascending node \(\Omega\) from: $$ \Omega - \Omega_0 = \dot \Omega (t - t_0), $$ the argument of perigee \(\omega\) from: $$ \omega - \omega_0 = \dot \omega (t - t_0), $$ and the mean anomaly from: $$ M - M_0 = n (t - t_0) $$

2 - Compute position and velocity as in the Keplerian case

The Kepler's equation is used to compute the eccentric anomaly \(E\) from the mean anomaly: $$ E - e\sin E = M $$ while the true anomaly \(\theta\) can be obtained from the relation: $$ \tan \frac E2 = \sqrt{\frac{1-e}{1+e}} \tan\frac \theta 2 , $$ Compute the flight path angle \(\gamma\) from: $$ \tan \gamma = \frac{e \sin\theta}{1 + e\cos\theta}, $$ the norm of the radius vector from: $$ r = \frac{a(1-e^2)}{1+e\cos\theta} $$ and the velocity norm from: $$ v = \sqrt{\frac {2\mu}r - \frac{\mu}a}, $$

The cartesian coordinates of the position vector \(\mathbf r\) and velocity vector \(\mathbf v\) can then be computed from:

$$ \begin{array}{l} x = r[\cos(\theta + \omega) \cos \Omega - \sin(\theta + \omega) \cos i \sin \Omega] \\ y = r[\cos(\theta + \omega) \sin \Omega + \sin(\theta + \omega) \cos i \cos \Omega] \\ z = r[\sin(\theta + \omega) \sin i] \\ v_x = v[- \sin(\theta + \omega - \gamma) \cos \Omega - \cos(\theta + \omega - \gamma) \cos i \sin \Omega] \\ v_y = v[- \sin(\theta + \omega - \gamma) \sin \Omega + \cos(\theta + \omega - \gamma) \cos i \cos \Omega] \\ v_z = v[\cos(\theta + \omega - \gamma) \sin i] \\ \end{array} $$