# Numerical Solution for the Single-Impulse Flyby Co-Orbital Spacecraft Problem
## 3 Mathematical Formulation of the Problem  
### Nonlinear Equations for Terminal Constraints  
For each spacecraft in its orbital coordinate system, its position can be obtained using the following formula:  
\[r_{i,t_{k}}=\frac{p_{i}}{1+e_{i}\cos\theta_{i,t_{k}}}(i,j=0,1,2;k=s,0,1,2) \tag{2}\]  
where \(\theta_{i,t_{k}}\) is the true anomaly of its orbit, \(p_{i}\) is the semi-latus rectum of the current orbit, and \(e_{i}\) is the eccentricity of the current orbit. All variables and their descriptions can be obtained in Glossary. Simultaneously, since all three orbits are located in the same plane, there is a conversion relationship between the true anomaly of the same position in various orbits. As shown in Figure 2, the angle between the interception's first position and the perigee of the chaser's orbit can therefore be found as follows:  
\[\theta_{0,t_{1}}=\theta_{1,t_{1}}+\Delta\omega_{o_{2},o_{1}} \tag{3}\]  
where \(\Delta\omega_{o_{0},o_{1}}\) is the angle between the perigee of the target's orbit and the perigee of the chaser's orbit. This can be provided by the eccentricity vector of the chaser's orbit and target's orbit since  
\[\Delta\omega_{o_{2},o_{1}}=\arccos\biggl{(}\frac{\mathbf{e_{2}}\cdot\mathbf{e _{1}}}{|\mathbf{e_{2}}|\cdot|\mathbf{e_{1}}|}\biggr{)} \tag{4}\]  
Subsequently, considering the relationship between the eccentricity vector and the argument of periapsis, it can also be calculated using  
\[\Delta\omega_{o_{2},o_{1}}=\omega_{o_{1}}-\omega_{o_{2}} \tag{5}\]  
Meanwhile, Equation (2) may be used to write the two crossings of the transfer orbit and the target orbit. \[ \begin{cases}r_{0,t_{1}}=\frac{p_{2}}{1+e_{2}\cos\theta_{0,t_{1}}}\\ r_{1,t_{1}}=\frac{p_{1}}{1+e_{1}\cos\theta_{t_{1}}}\end{cases}  \tag{6}\]  
Substituting Equation (3) into the first equation of Equation (6) and combining it with Equation (1) gives  
\[\frac{p_{1}}{1+e_{1}\cos\theta_{1,t_{1}}}=\frac{p_{2}}{1+e_{2}\cos(\theta_{1, t_{1}}+\Delta\omega_{o_{2},o_{1}})} \tag{7}\]  
Organizing this gives  
Figure 2: Geometric interpretation of flyby positions. \[(e_{1}p_{2}-e_{2}p_{1}\cos\Delta\omega_{o_{2},o_{1}})\cos\theta_{1,t_{1}}+e_{2}p_{ 1}\sin\Delta\omega_{o_{2},o_{1}}\sin\theta_{1,t_{1}}=p_{1}-p_{2} \tag{8}\]  
Simplifying this expression using trigonometric functions, when \(\mathbf{e}_{2}\
eq 0\) and \(\Delta\omega_{o_{2},o_{1}}\
eq k\pi(k=0,1)\), the equation is obtained as follows:  
\[\left\{ \begin{array}{l}A\cos\theta_{1,t_{1}}+B\sin\theta_{1,t_{1}}=C\\ A=e_{1}p_{2}-e_{2}p_{1}\cos\Delta\omega_{o_{2},o_{1}}\\ B=e_{2}p_{1}\sin\Delta\omega_{o_{2},o_{1}}\\ C=p_{1}-p_{2}\end{array} \right. \tag{9}\]  
Then, two solutions of Equation (9), which represent the true anomaly of the two interception positions in the target's orbit, can be obtained for the intersecting ellipse equations. Then, considering the laws of trigonometric functions, the logic between those two solutions can be expressed more concisely as follows:  
\[\left\{ \begin{array}{l}\theta_{1,t_{1}}+\theta_{2,t_{2}}=\pi\quad\quad\quad \quad A=0,B\cdot C>0\\ \theta_{1,t_{1}}+\theta_{2,t_{2}}=3\pi\quad\quad\quad\quad A=0,B\cdot C<0\\ \theta_{1,t_{1}}+\theta_{2,t_{2}}=2\pi\quad\quad\quad B=0\\ \theta_{1,t_{1}}+\theta_{2,t_{2}}=\pi-2d\quad\quad A\cdot B\
eq 0,C>0\\ \theta_{1,t_{1}}+\theta_{2,t_{2}}=3\pi-2d\quad\quad A\cdot B\
eq 0,C<0\end{array} \right. \tag{10}\]  
where  
\[d=\arctan\frac{A}{B},d\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right] \tag{11}\]  
The relationship between \(\theta_{2,t_{2}}\) and \(\theta_{1,t_{1}}\) is the same as the relationship shown in Equation (12). Thus, once the true anomaly of the first flyby target in the target orbit is determined, the true anomaly of the second flyby target is also obtained. Therefore, the mean anomaly can be calculated using  
\[\left\{ \begin{array}{l}\tan\frac{F}{2}=\sqrt{\frac{1-\varepsilon}{1+ \varepsilon}}\tan\frac{\theta}{2}\\ M=E-e\sin E\end{array} \right.