Abstract
Two systems of Earth-centered inertial Newtonian orbital equations for a spherical Earth and three systems of post-Newtonian nonlinear equations, derived from the second post-Newtonian approximation to the Earth Schwarzschild field, are used to carry out a performance analysis of a numerical procedure based on the Dormand-Prince method for initial value problems in ordinary differential equations. This procedure provides preliminary post-Newtonian corrections to the Newtonian trajectories of middle-size space objects with respect to space-based acquisition, pointing, and tracking laser systems, and it turns out to be highly efficient. In fact, we can show that running the standard adaptive ode45 MATLAB routine with the absolute and relative tolerance, TOLa = 10−16 and TOLr = 10−13, respectively, provides corrections that are final within the eclipses caused by the Earth and close to final during the noneclipse phases. These corrections should be taken into account to increase the pointing accuracy in implementing the space-to-space laser links required for ablation of designated objects or communications between space terminals.
1. Introduction
In a previous paper [1], the main post-Newtonian (p-N) corrections to the Newtonian trajectories of middle-size LEO space debris object D with respect to acquisition, pointing, and tracking (APT) laser system S, aimed at ablating D, were derived by means of a genuine p-N method for the Earth Schwarzschild field [2]. After the Earth-centered inertial (ECI) orbital coordinates of S and D are determined from appropriated p-N orbital equations, the main difficulty of this method is that the process to track D involves computations of the time-dependent coefficients of a third system of p-N equations for the relative motion of D with respect to S. These coefficients are line integrals of the Riemann tensor of space along the line of sight (LOS) segment that joins S and D (cf. [1]).
Despite its complexity, this method is revealed to be suitable for solving this particular problem, whose detailed description can be found in [3–9] (see also [10–12]). Namely, the method is p-N orthodox and works dependably, once the convergence of the integrals has been verified. This is the case when the distance d(S,D) from S to D is small compared to the altitude of D above the Earth surface, so that D may always be kept in the LOS from S during the action. However, this method does not work correctly in scenarios where d(S,D) is so large that, due to the difference of orbital speeds of S and D, an eclipse of D from S caused by the Earth can occur. Such a situation is laser communication between two distant space terminals. Here, the only available option to compute the p-N corrections for the relative motion of D with respect to S during the eclipse phases is to adopt the translation Euclidean rule (TER) for the respective ECI positions. These positions have to be computed from the ECI p-N geodesic orbital equations of S and D derived via the second p-N approximation to the Earth Schwarzschild field. We are entitled to adopt this rule, since the Riemann parallel transport in this approximation is the second p-N approximation to the Euclidean transport.
Therefore, in those scenarios where eclipses may occur, one has to consider the two expected alternatives, eclipse/noneclipse, by means of the “if/else” computational command, and make the respective p-N relative trajectories match. This can be done by applying the ode45 solver in a stepwise manner in order to be able to distinguish the eclipse from the noneclipse phases [13, 14]. That way, we can verify if we deal with an eclipse phase, once the complete interval of integration has been split into a number of appropriately small subintervals. From now on, this procedure is referred to as the stepwise ode45 method (SODE).
In this paper, we continue the work initiated in [14] and complete the first step in the process of estimating the p-N corrections to the Newtonian trajectory of D with respect to S. More precisely, we compute the values of the corrections by means of the TER, which can be considered as final within the eclipse phases and as preliminary for the noneclipse phases.
2. The Method
To complete the first step, we take the Newtonian ECI positions derived with the Newtonian orbital equations for S and D for a spherical Earth during the total time required for a given action. First, as in [14], we apply the standard adaptive ode45 MATLAB routine, which includes two explicit Runge-Kutta methods of order four and five. Then, we repeat the calculation using the new SODE method in order to integrate the ECI Newtonian orbital equations for S and D. For the numerical simulation, we choose decreasing values of the absolute and relative tolerances and run the MATLAB routine (on a whole time interval, as well as on a sequence of subintervals) until the accuracy requirements are satisfied. The final tolerance pair shall guarantee that the difference between the two computed ECI positions is on the subcentimeter level. In the sequel, this condition is referred to as “error condition.”
This final pair of suitable tolerances is now used in all procedures, in particular, in the SODE method including the TER applied to integrate the p-N equations for the relative motion of D with respect to S. Hence, since the TER is the only available option for the eclipse phases, we can conclude that these p-N corrections are final for the eclipse phases and close to final for the noneclipse phases.
Finally, let us stress that the p-N equations used in the SODE procedure are suitable to derive the main p-N corrections for the relative motion of D with respect to S. This is due to the fact that for the geometric model of the Earth surrounding space, we take the second p-N approximation to the Earth Schwarzschild field. The approximation in ECI coordinates, , is given by (Latin indices range from 1 to 4 and Greek from 1 to 3) where
Here, is the mass of the Earth given in seconds, , the Euclidean distance from the space object at hand to the center of the Earth, also given in seconds, and, , a small dimensionless parameter of the order of . Note that the metric has the form (1) and (2), since we set . Moreover, we specify all quantities in seconds, because time is the basic measurement quantity in accurately measuring distances in space and in finding p-N corrections in the geolocation of objects orbiting the Earth [15]. Hence, by adopting the equivalence, we take for the value and for the mean radius of the Earth, .
Then, from the geodesic principle for curved spaces, it follows that the equations suitable to describe the ECI p-N motions in the space-time given by (1) and (2) are from which we have
Thus, after omitting the expressions in (3) and (4) to facilitate the reading, we finally arrive at the following equations, where , , and :
A detailed derivation of these equations can be found in [14].
3. Results and Discussion
In order to derive the results below, the ECI orbital motions of S and D are assumed to be equatorial, so that the number of equations which we need to solve for the p-N orbit of D with respect to S reduces from eighteen to twelve. Obviously, this assumption does not cause any loss of generality, since the Earth Schwarzschild field is spherically symmetrical, but the computations become much easier.
Consequently, the inputs used to determine the ECI p-N orbits of S and D are (i)Their respective altitudes at perigees, in seconds(ii)The distance , between perigees, in seconds(iii)The distance between the initial positions of S and D, also in seconds(iv)The orbital eccentricities, and , of S and D, respectively
Note that, when the arguments of the perigees of S and D are equal, we have .
Then, it follows by the TER that the equations for the relative motion of D with respect to S are where , , and and are the p-N ECI coordinates, in seconds, of D and S, respectively.
Taking into account the equivalence between length and time, we begin our numerical simulation with the relative and absolute tolerances TOLr = 10−11 and TOLa = 10−11, as in [14]. In particular, we see that the estimates for the p-N corrections with these tolerances are not sufficiently accurate for the present aim, since the error condition is not satisfied. Therefore, the absolute and relative tolerances are gradually decreased in the course of the numerical simulation until we reach the desired accuracy. It turns out that in all considered scenarios, the error condition is satisfied for the first time for the pair of tolerances TOLr = 10−13 and TOLa = 10−16.
The inputs for the scenario in [14] were , , , , and . This scenario is quite representative of LEO-LEO configurations. For this reason, we use the same values, although here we study the long-time behavior of the solutions. The global time is (seven days), which is much larger than the one in [14] because our main interest is in observing the stability of the code. We are aware that the results most probably stop to deliver dependable corrections after two days due to overshooting effects.
Since for all pairs of tolerances the related pictures look very similar, we decided to select for the figures below the pairs TOLr = 10−13 and TOLa = 10−13 and TOLr = 10−13 and TOLa = 10−16. The results shown in the figures below are given in km, or in cm, because the units of length are more familiar in this context.
To begin with, we show in Figure 1 the virtual LOS, represented by the straight line segments from S (the outer orbit, blue) to D (the inner orbit, red), while S and D are orbiting the Earth (green). The figure shows the behavior of D with respect to S, with the true LOS beginning on the outside right, which corresponds to the shortest segment S-D. The eclipse of D from S starts when the segment S-D begins crossing the Earth, the corresponding LOS is one of the slanted segments, almost vertical on the left. The figure corresponds to the first two hours of orbital motions. It shows that already during this time, the first eclipse arises and also that it has not finished yet. In fact, the last segment S-D indicates that D is opposite to S with respect to Earth.
Figure 2 shows the p-N relative motion of D with respect to S, within the first two days as a part of the full computation of the relative position of D with respect to S for seven days. The SODE method for this experiment was executed on 3000 time subintervals. The orbit runs counter clockwise. Note that it is not periodic.
Figure 3 shows the difference between the coordinates of the relative Newtonian positions of D with respect to S computed with the ode45 MATLAB and the SODE procedure. The computation corresponds to seven days and the tolerances TOLr = 10−13 and TOLa = 10−13. It can be seen that the difference increases very slowly along time, which looks promising in terms of stability. However, the error condition is not satisfied, because the errors and become larger than one cm from the third day on.
(a)
(b)
Figure 4 shows the radial and transversal p-N corrections versus the distance from S to D, computed again with the SODE procedure and the TER but with tolerances TOLr = 10−13 and TOLa = 10−13, and Figure 5 shows the total p-N corrections versus the time with the same tolerances. Here, the overshooting effects are noticeable.
(a)
(b)
Figures 6–8 are analogous to Figures 3–5, respectively, but now, the tolerances are TOLr = 10−13 and TOLa = 10−16. It follows from Figure 6 that the values for the corrections shown in Figures 7 and 8 can be considered sufficiently accurate approximations to the true p-N corrections to the relative position of D with respect to S, since the error condition is satisfied.
(a)
(b)
(a)
(b)
Finally, Figure 9 shows that the difference between the relative Newtonian positions of D with respect to S computed with the ode45 MATLAB and the SODE method satisfies the error condition during all seven days. This means a very satisfactory stability of the code.
4. Conclusions
The p-N corrections computed with the SODE method including the TER could be implemented to increase the accuracy of the relative Newtonian trajectory of D with respect to S. In fact, this procedure is operative and ready to be used within the eclipse phases of D with respect to S, although for the time being, it does not cover all important real-life applications within the noneclipse phases, such as hypothetical APT laser system aiming at ablating space debris. In order to reach accuracies necessary for this and similar applications, the p-N corrections to the Newtonian relative motion of D with respect to S should be computed using the p-N orthodox method described in [1] for the noneclipse phases. The technical difficulty now is that these equations involve integrals (cf. [1]) whose accurate numerical approximation may be difficult. This demanding problem will be in the focus of our future work.
Note that the scenarios considered in this paper are representative for any LEO-LEO configuration and the results obtained here support the recommendation given in [16]. Since we cover scenarios with eclipse phases, we are now in the position to carry out a comprehensive analysis of LEO-MEO and LEO-GEO configurations.
Data Availability
The data used to support the findings of this study are included within the article.
Disclosure
The research of the present paper was performed as part of the employment of the authors. The data of the employers are in the front page.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.