## 3 Orbit Determination with ISL Measurements  
Let us assume to have a system of \(N>1\) orbiters in a Halo orbit. This system can be easily described by exploiting the Circular Restricted Three-Body Problem (CRTBP) where two massive bodies are considered orbiting around their mutual center of mass, and the spacecraft with infinitesimal masses experience forces due to the gravitational influence of both bodies simultaneously. All the bodies are approximated as point masses. Moreover, to simplify the notation, let us assume to use a reference frame centered in the two celestial body center of gravity and keep the \(x\)-axis pointed towards the second celestial body, obtaining a frame rotating with constant angular velocity, equal to the mean motion of the system. Finally, let us scale the system in such a way that the distance between the two main bodies, as well as the angular velocity, are set to 1 (see Figure 3). In this coordinate frame, it is possible to define the state vector for any spacecraft (containing its position and velocity) as:  
\[\mathbf{X}=[X,Y,Z,\dot{X},\dot{Y},\dot{Z}]^{T} \tag{14}\]  
or, in the general case of \(N\) spacecraft as:  
\[\mathbf{X}=[X_{1},Y_{1},Z_{1},\dot{X}_{1},\dot{Y}_{1},\dot{Z}_{1},\ldots,X_{N},Y_{ N},Z_{N},\dot{X}_{N},\dot{Y}_{N},\dot{Z}_{N}] \tag{15}\]Moreover, in this reference frame for the generic spacecraft \(i\), the equations take the following form: [30]:  
\[ \begin{split}\ddot{X}_{i}&=\frac{\partial U_{i}}{ \partial X_{i}}+2\dot{Y}_{i}\\ \ddot{Y}_{i}&=\frac{\partial U_{i}}{\partial Y_{i}}- 2\dot{X}_{i}\\ \ddot{Z}_{i}&=\frac{\partial U_{i}}{\partial Z_{i}} \end{split}  \tag{16}\]  
where \(U_{i}\) is the potential function:  
\[U_{i}=(X_{i}^{2}+Y_{i}^{2})/2+(1-\mu)/d_{1,i}+\mu/d_{2,i} \tag{17}\]  
and  
\[ \begin{split}\mu&=m_{2}/(m_{1}+m_{2})\\ d_{1,i}&=\sqrt{(X_{i}+\mu)^{2}+Y_{i}^{2}+Z_{i}^{2}} \\ d_{2,i}&=\sqrt{(X_{i}+\mu-1)^{2}+Y_{i}^{2}+Z_{i}^{2} }\end{split}  \tag{18}\]  
and \(m_{1}\) and \(m_{2}\) are the celestial body masses. In this particular frame, a batch-weighted least square orbit determination algorithm [30] can be used to estimate the initial state of all the spacecraft. The differential equation of the satellite motion, \(\dot{X}=\mathbf{F}(\mathbf{X},t)\), can be derived from the motion law and can be used to propagate the initial state to any possible time \(t\) by integration. It follows that, knowing the initial state of the system, it is possible to derive the state, \(\mathbf{X}(t)\), for any possible time \(t\). In the same way, having a reference solution (that is a guess with an error of the initial state), it is possible to also compute a reference state for the time \(t\), \(\mathbf{X}^{*}(t)\). Finally, for any time \(t\), it is also possible to define the expected measurement vector \(\mathbf{Y}^{*}(t)\) for the reference solution and the real measurement vector for \(\mathbf{Y}(t)\). The vector \(\mathbf{Y}(t)\) contains the set of the range and range-rate measurements taken at time \(t\) as defined in Equation (3) and in Equation (6). Making the hypothesis conduct the measurement process with an update time of \(\Delta t\) for a Time span of \(T\), it is possible to compute the state deviation, \(\mathbf{x}(t_{k})\), and the observation residual, \(\mathbf{y}(t_{K})\), for each measurement instant \(k=1\dots M\):  
Figure 3: The reference frame used to solve the OD problem.