International Journal of Aerospace Engineering

Volume 2019, Article ID 3723018, 7 pages

https://doi.org/10.1155/2019/3723018

## Numerical Approach for the Computation of Preliminary Post-Newtonian Corrections for Laser Links in Space

^{1}Gregorio Millan Institute, Univ. Carlos III de Madrid, 28911 Leganés, Madrid, Spain^{2}I.E.S. Alpajés, 28300 Aranjuez, Madrid, Spain^{3}Institute for Analysis and Scientific Computing, Vienna University of Technology, A-1040 Vienna, Austria

Correspondence should be addressed to Jose M. Gambi; se.m3cu.htam@ibmag

Received 30 August 2018; Accepted 17 October 2018; Published 9 January 2019

Academic Editor: Franco Bernelli-Zazzera

Copyright © 2019 Jose M. Gambi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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, TOL_{a} = 10^{−16} and TOL_{r} = 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 TOL_{r} = 10^{−11} and TOL_{a} = 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 TOL_{r} = 10^{−13} and TOL_{a} = 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 TOL_{r} = 10^{−13} and TOL_{a} = 10^{−13} and TOL_{r} = 10^{−13} and TOL_{a} = 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.