An analytic solution of full-sky spherical geometry for satellite relative motions

Herein, an exact and efficient analytic solution for an unperturbed satellite relative motion was developed using a direct geometrical approach. The derivation of the relative motion geometrically interpreted the projected Keplerian orbits of the satellites on a sphere (Earth and celestial spheres) using the solutions of full-sky spherical triangles. The results were basic and computationally faster than the vector and plane geometry solutions owing to the advantages of the full-sky spherical geometry. Accordingly, the validity of the proposed solution was evaluated by comparing it with other analytic relative motion theories in terms of modeling accuracy and efficiency. The modeling accuracy showed an equivalent performance with Vadali’s nonlinear unit sphere approach, which is essentially equal to the Yan–Alfriend nonlinear theory. Moreover, the efficiency was demonstrated by the lowest computational cost compared with other models. In conclusion, the proposed modeling approach illustrates a compact and efficient closed-form solution for satellite relative motion dynamics.

In previous works, the analyst was forced into analytic approaches using vector terms rather than spherical geometry because of the singularity problem. However, these problems have been eliminated by solutions of full-sky triangles and plugged into computer programming languages for many years. This study aims to develop a solution for relative motion dynamics by direct geometrical interpretation using the advantages of full-sky spherical geometry solutions. Spherical geometry now makes computer implementation quite simple, and this has been discussed in Ref. 18 .
To derive the equations of motion, the approach geometrically interprets the projected base and the target satellite Keplerian orbits on a sphere (Earth and celestial spheres) with the spherical trigonometry solutions. The resulting equations are expressed as azimuth and elevation angles with their derivatives representing the relative angular position and velocity of the target satellite with respect to the base satellite. Subsequently, the azimuth and elevation angles are then transformed into the associated rectangular coordinates for the relative position and velocity vectors in the base satellite centered frame. The proposed solution is validated by modeling accuracy and efficiency compared with those of the exact analytic dynamic models of the satellite relative motion theories.

Keplerian orbit in spherical coordinate systems
We project a Keplerian orbit on a celestial sphere using spherical coordinates. The Keplerian orbit is commonly specified by the classical orbital elements for state representations in space 19 . The six orbital element sets are as follows: The semi-major axis, a, and the eccentricity, e, listed as the first two elements above describe the orbit size and shape, respectively, while Ω, i, and ω define the orbit plane orientation. The final classical orbital element is the true anomaly, ν, which determines the object's current angular position relative to the perigee. Figure 1 illustrates the orbit elements Ω, i, ω, and ν, which are angle-related orbit elements that describe the Keplerian orbit from the center of the Earth. A spherical coordinate system in space can be used to represent the Keplerian orbit projected on a sphere. In Fig. 2, I axis is along the vernal equinox direction and the K axis is in the direction of the north pole, and then the elevation angle, δ ′ , defines the angle between the straight line from the center of the Earth to O ′ and the projection of this line on the I J planeThe angle between this projection and the I axis is defined as the azimuth angle, α ′ . The formulas for α ′ and δ ′ can be expressed as follows in terms of the angle-related orbital elements if we represent the projected Keplerian orbit by α ′ and δ ′ : Angles α ′ and δ ′ and the radial distance, r, of the object determine the spherical coordinate system in space. The radial distance, r, is written as follows in terms of ν: If we represent the position of the object in the rectangular coordinate system, the spherical coordinate system is transformed by the following relations: www.nature.com/scientificreports/ Angles α ′ and δ ′ are expressed in the next section in terms of the angle-related orbit elements using direct geometrical interpretations.

Geometrical relative orbit modeling
In this section, we geometrically derive the relative position and the velocity vectors of a target satellite relative to the base satellite. The subscript B denotes the base satellite while subscript T denotes the target satellite. The Keplerian orbits of the two satellites are projected on a sphere for geometrical interpretation, as seen in Fig. 3. P B and P T denote the orbit poles of the satellites. The dotted lines on the sphere represent the projected Keplerian orbits of the two satellites and the solid line represents an equatorial plane. The intersection point, I P , is defined as the projected crossing point of the two orbit planes on the surface.

Note 1
In spherical geometry, intersection points always exist between great circle arcs on the sphere if they do not have the same inclinations.
The relative position of the target satellite, T, with respect to the base satellite, B, is given by the azimuth angle, α, and the elevation angle, δ. Angle α is perpendicular to angle δ through point H.  where the arc lengths φ j and θ j represent the distance from the ascending nodes, j , to the intersection point, I P , and from I P to the satellite's current angular position, respectively. For the derivation of the satellite relative motion, a key parameter is the relative inclination, i R , which is the angle between two orbit planes at I P . We use the spherical triangle B T I P to compute i R . Because i R is not equal to the difference between two inclinations of the orbits (i.e., i R = i T − i B ), we must apply the law of cosines for angles to the triangle: where, the relative ascending node, ∆Ω, is defined as follows: We first derive α and δ for the target satellite relative to the base satellite in terms of the angle-related orbit elements Ω, i, ω, and ν. The general solutions for spherical triangles are given in Ref. 18 . The spherical triangle B T I P in Fig. 3 is taken to solve φ B and φ T . Figure 4 shows a detailed view of the spherical triangle. We apply the law of sines to the spherical triangle to compute sinφ B : Applying the law of cosines for angles to the spherical triangle B T I P , we find another geometrical relationship to compute cosφ B : Dividing Eqs. (8) by (9) and applying the full-sky trigonometry solutions give To compute sinφ T , the law of sines is also applied to the spherical triangle seen in Fig. 4: Using the law of cosine for angles, we also obtain Dividing Eqs. (11) by (12) results in Now we consider the spherical triangles on the surface of the sphere with ∆Ω = 0. In this case, we construct a celestial sphere having the pole P B of the base satellite as a geographical pole, shown in Fig. 5. The celestial sphere has two spherical triangles, P B P T T and THI P . Note that the angle P B P T I P is always 90°, regardless of the inclination of the satellites. The angle T P B I P is equivalent to the angle θ T by applying the law of sines. Hence, the angle P B P T T is obtained by subtracting θ T from 90°. The arcs P T T and P B H are always 90°. Thus, the arc P B T can be found by subtracting δ from 90°.
The elevation angle δ is derived from the spherical triangle P B P T T . By applying the law of cosines for sides to the spherical triangle, we find that The azimuth angle α is found by applying the law of cosines for sides twice to the spherical triangle THI P , resulting in the following equations: and Substituting cosδ from Eqs. (16) into (17), we obtain Thus, angle α applying the full-sky trigonometry solutions is derived by Using the definition of the argument of latitudes in Eq. (5), angles α and δ are expressed as

Note 2
The quadrant ambiguity problem is avoided by using the ATAN2 built-in function in computer programming languages for φ B , φ T , α , in order to cover the full range of 0° to 360°.
For a simple analysis of satellite relative motion, angles α and δ can be directly used to determine the angular position of the target satellite with respect to the base satellite. In Eq. (20), ν is the only time-dependent variable, and the derivatives of α and δ are obtained using: The angle δ in Eq. (20) is a composite function and the derivative of δ for Eq. (21) is computed using the chain rule. The relative motion of the target satellite can be described using the previously calculated α and δ in the rectangular coordinates. The orbit radius of the base satellite is r B and the target satellite orbit radius is www.nature.com/scientificreports/ r T . We introduce the base satellite rotating frame, F R , to describe the relative motion of the target satellite with respect to the base satellite. The center of the Earth is set as the origin, and the orientation of F R is given by the unit vectors e 1 , e 2 , and e 3 . The direction of the unit vector e 1 is set to the orbit radius of the base satellite, while e 3 is perpendicular to the orbit plane of the base satellite. The unit vector e 2 is then computed following the righthand rule. F R is mathematically described by the unit vectors as follows: Superscript × denotes a skew-symmetric 3 × 3 matrix associated with a 3 × 1 column matrix. If x is a 3 × 1  20) and (21). The only assumption is that no perturbations are acting on the satellites.
The resulting solutions, we collectively call herein as geometrical relative orbit modeling (GROM), are obtained by the direct geometrical method using full-sky spherical geometry. GROM has the following explicit advantages: first, the angular values ( α, δ,α,δ ) are useful for the angle-only measurement analysis or tracking problems between satellites without any coordinate transformation; and second, the solutions provide a comparatively simpler form for the dynamic model of relative orbits. For comparison, this study introduced the USA (see Online Appendix for this approach) 10 , which nonlinearly maps orbital elements to the relative position and velocity in the same manner as GROM. The two solutions consisted of trigonometric functions as orbital element terms. Table 1 presents a comparison of the formula complexity of the solutions.
The GROM is expressed in compact form with the ATAN 2 inverse trig function, and the USA represents analytic expressions in terms of differential orbital elements, along with the direction cosine matrix of the chief and deputy satellites.

Model evaluations
We first evaluate the modeling accuracy to determine the validity of the GROM solution. To evaluate the modeling accuracy, three relative motion theories are introduced: Hill's equation 19 , the classical two-body problem (CTBP) 20 , and the USA. The USA has been evaluated to have the lowest index in terms of accuracy of satellite relative motion equivalent to the Yan-Alfriend nonlinear method 8,10 .
Model accuracy. This section evaluates the relative errors of GROM using the modeling error index which is an effective tool for evaluating the accuracy of relative motion theories 8  where R e is the Earth's radius of 6378.14 km and n is the mean motion of satellite. The modeling error index is written as follows: The modeling error index evaluates the relative errors of GROM in comparison to Hill's equation and USA, relative to the reference orbit model of CTBP. The analytic solutions of USA and CTBP describe the kinematically exact relative motion of satellites, which brings about negligible quantities of relative errors upon simulation. We use k-digit rounding arithmetic when finding the solutions, y j and − y j , to ignore computational uncertainties such as round-off errors. The k-digit rounding arithmetic is obtained by terminating the value of the solution at k decimal digits.
For numerical simulations, the orbit elements of the base satellite in Table 2 are chosen, and the orbit element differences, oe , of the target satellite are used as the following values: Figures 6 and 7 show the modeling error index using the six-digit rounding arithmetic solution with various relative distances and eccentricities. As shown in the figures, the index of GROM is exactly the same as that of USA, representing index values of 10 −6 with respect to the reference orbit model. A modeling error index of the order of 10 −3 is sufficiently small, providing a reasonable confidence regarding the modeling 21 .
Therefore, the GROM solution provides an accurate representation for all relative orbit sizes and eccentricities. In the case of the Hill's equation, the solution shows modeling indexes of nearly 10 1 at small orbit element differences, which means sufficient accuracy for small relative orbit sizes and eccentricities. As expected, however, the index values gradually grow with increasing orbit size and eccentricity.
Model efficiency. This section compares CTBP of a vector solution, and USA of a nonlinear mapping solution similar to GROM, to evaluate the computational cost of GROM. In tests, we used MATLAB (R2013b) running on an 8th Intel Core i5 8265U (3.4 GHz) with 8 GB memory. The CPU time comparison of each solution is clearly implementation-dependent on the researcher. To fairly evaluate the comparison of solutions, each solution was coded efficiently, and a simple least squares linear regression and method was applied for a more robust comparison. Generally, the relationship between iterations, ∼ N, and CPU time, ∼ T, is roughly linear. However, it is not exact linear. Along with the parameters in Table 2, five iterations, were performed with a time step 0.1 s and a time span(sec) below: The computational cost of each solution is compared by the slope m value and the least squares method is used to calculate the slope of the best fit line. Table 3 shows the coefficient m * , normalized with respect to the value of GROM. GROM is nearly 7% and 25% more efficient than USA and CTBP, respectively.
A numerical example was used to study the effect of the satellite relative motion under the influence of J 2 perturbations, through which we demonstrated the modeling efficiency of GROM. Using the time-explicit orbital elements in the analytic solutions provides a simple method of investigating the difference between unperturbed and J 2 -perturbed models for the satellite relative motion. The first-order J 2 perturbation causes secular changes (33) T k = mÑ k , k = 1, . . . , n   www.nature.com/scientificreports/ in the ascending node, Ω, argument of perigee, ω, and mean anomaly, M. The time-explicit representations of the first-order J 2 perturbation are presented as follows 10 : We use the following values in the time-explicit orbit elements: p = a(1 − e 2 ), J 2 = 0.00108263. The numerical simulation coded to iterate the solutions at each time step runs over a period of 20 days using the orbit elements for the base satellite shown in Table 2. The orbit element differences are chosen as the following values: Table 4 shows that the maximum differences between the relative position and velocity of the solutions are 3.8729 km and 0.0041 km/s over 20 days, respectively. However, each of the analytic solution results in different CPU times; thus, the GROM solution is approximately 1 ′ 40 ′′ and 5 ′ 57 ′′ faster than the USA and CTBP solutions, respectively.

Conclusions
Herein, we developed an analytic solution for the unperturbed satellite relative motion using a direct geometrical approach on a sphere. The derivation of this geometrical approach is straightforward, and the resulting solutions provide a complete analytic form for relative motion dynamics with the full-sky solutions covering the range of 0°-360°. The accuracy of the proposed GROM solution equals that of the exact relative motion theories of CTBP and USA. In particular, in terms of the modeling efficiency, GROM is approximately 7% and 25% more efficient than USA and CTBP, respectively. Consequently, using full-sky spherical geometry, the proposed GROM approach facilitates simpler equations of motion and higher computational speed compared with the other methods.