Increasing the Accuracy of Orbital Elements for a Satellite in a Low Earth Orbit under the Influence of Atmospheric Drag Using Adams-Bashforth Method

The perturbed equation of motion can be solved by using many numerical methods. Most of these solutions were inaccurate; the fourth order Adams-Bashforth method is a good numerical integration method, which was used in this research to study the variation of orbital elements under atmospheric drag influence. A satellite in a Low Earth Orbit (LEO), with altitude form perigee = 200 km, was selected during 1300 revolutions (84.23 days) and ASat / MSat value of 5.1 m 2 / 900 kg. The equations of converting state vectors into orbital elements were applied. Also, various orbital elements were evaluated and analyzed. The results showed that, for the semi-major axis, eccentricity and inclination have a secular falling discrepancy, Longitude of Ascending Node is periodic, Argument of Perigee has a secular increasing variation, while true anomaly grows linearly from 0 to 360°. Furthermore, all orbital elements, excluding Longitude of Ascending Node, Argument of Perigee, and true anomaly, were more affected by drag than other orbital elements, through their falling as the time passes. The results illustrate a high correlation as compared with literature reviews in this field.


Introduction
The main condition in orbit's determination is to find state vectors and orbital elements. During the rotation of the satellite around the Earth, it is affected by several types of perturbations, which are deviations from normal motion [1,2]. In the presence of perturbations, the ideal path of a satellite will vary from the hypothetical two-body problem [3]. They are classified into two categories: gravitational and non-gravitational. Gravitational perturbations include those of the non-spherical Earth and the third body attraction, whereas non-gravitational perturbations involve those of the atmospheric drag and solar radiation pressure. The process of orbit determination is based on the predication of the satellite's drift, which in turn is used to estimate the future orbit [4,5,6]. This process occurs by a numerically integrated perturbed equation of motion, using one of the integration methods, such as Adams-Bashforth. It is strongly affected by the orbit dynamics and, thus, quality control of the numerical methods is required prior to each orbit analysis [7].
Many researchers were interested in studying various numerical integration methods to find an exact orbit determination. Shen et al. estimated the trajectory of the satellite under the influence of J 2, J 3 , and J 4 via fourth order Runge-Kutta method [8]. Mishra et al. elaborated the method of determining a precise ephemeris for an orbiting satellite due to atmospheric drag. This method includes guessing of state vectors from a sequence of data [9]. Mohammed and Abdul-Rahman extensively studied the behavior of orbital elements using changeable values for A Sat. / M Sat. , altitudes, and eccentricities due to different perturbations [10]. Thangavel et al. studied the effects of atmospheric drag and J 2 in a closeproximity operation for LEO [11].

Methods
The following steps are used to complete the requirements of this research.

Kepler's Equation Solution
This equation is used to find eccentric anomaly E then true anomaly f. Later, both these two parameters are employed for describing the position and velocity (state vectors) at different times. Mikkola's method is one of the methods to solve Kepler's equation with an accuracy of 10 -10 . To start this method, an initial value for eccentric anomaly E o must be calculated as below [12]: where M is mean anomaly, e is eccentricity, , √ , . Hence, the final value of eccentric anomaly can be expressed as: To solve the equation above, e and M are given, then the equation must be solved directly for eccentric anomaly. The repetition will continue until the desired value of eccentric anomaly is obtained.

State Vectors Calculation
The position and velocity vectors of a satellite in its orbit can be represented as [2,13]: where μ = GM is the gravitational constant of the Earth (398600 km 3 /sec 2 ), f is the true anomaly, p is the semi latus rectum = a (1-e 2 )), h is angular momentum per unit mass, and a is semi-major axis. The elements of e, f, h and a will be explained later in the orbital elements calculation section.

Ibrahim and Saleh
Iraqi Journal of Science, 2021, Special Issue2, pp: 81-90 83 To convert the position and velocity components from the object plane to the equatorial plane, the Gaussian vector (conversion matrix) must be used as below [2]: where R -1 is the inverse of Gauss matrix, which contains Euler angles (inclination i, Argument of perigee ω and Longitude of Ascending Node Ω), as below [2]: Equations (9) and (10) can be represented as in the following: where X, Y, Z, X , Y, and Z are the state vectors of a satellite in the equatorial plane at time (t). Thus, the final illustration of state vectors can be stated as [13]: -Position is: r = -Velocity is: v = Orbital Elements Calculation  Semi-major axis can be defined as [13]:  Eccentricity can be calculated by [6,13]: ⁄ where is the distance from perigee point and is the distance from apogee point. In this study, only a semi-circular orbit is taken into account.  Inclination specifies the tilt of the orbit plane, which is represented as follows [13]: where h X , h Y , and h Z are the angular momentum per unit mass in X,Y, and Z directions, respectively. Therefore, the magnitude of h per unit mass is ( ) 0.5  Longitude of Ascending Node is [2]: Argument of perigee is stated as below [2]: where uu is the argument of latitude which can be computed by [3]: ⁄ And f is the true anomaly which can be calculated as [2]: ⁄ The input orbital elements applied in the present study are given in Table -1. The perturbed acceleration of the satellite due to atmospheric drag is [2,10]: where C D is the drag coefficient, A Sat. is the cross sectional area of the satellite perpendicular to velocity vector, M Sat. is the mass of satellite, v r is the satellite velocity vector relative to the atmosphere, and ρ is the atmospheric density. The density is calculated using the US Naval Research Laboratory mass spectrometer and incoherent scatter (NRLMSISE-00), where E indicates that the model extends from the ground through exosphere and 00 is the year of release. Table -1 gives the required values to input in equation (22). We suppose that the atmosphere rotates at the same angular speed as that of the Earth. With this hypothesis, the relative velocity vector is [10]: (23) where is the inertial velocity vector for satellite, is the Earth's rotational velocity vector, and r is the inertial satellite position vector. Further breakdown into vector components of equation (23) can be represented as [10]: The magnitude of v r is ( ) 0.5 . The perturbed equation of motion (under the atmospheric drag influence) is [2]: ̈

Fourth Order Adams-Bashforth Method
This method requires information about a series of previous points. Such methods need a boot-up to obtain enough information to start running. Adams-Bashforth is one of the multi-step methods; it is often referred to as a predictor-corrector method. There are other multi-step methods, such as Milne, Adams-Moulton, Gauss-J ackson (sum-squared), and Obrechkoff [15]. In order to predict the new positions and velocities by this method, equation (25) must be firstly integrated numerically using fourth order Runge-Kutta's method as below [2,7]: where is the initial velocity at epoch, is the, predicted velocity = ax o = ax o + = ax o + = ax o + , ax o is the acceleration at epoch, s the step time equal to . t step = T p /nn, is the sub-step number (the used value is 50), T P is the satellite's orbital period, and is the step number per one revolution. t step = T p /50000 is considered as the best value to obtain a better solution after many attempts, as compared to those reported in other literature reviews [10,14].
where is the initial position at epoch, is the predicted position, kk 1 = vx o , kk 2 = vx o + 0.5 σ step kk 1 , kk 3 = vx o + 0.5 σ step kk 2 , and kk 4 = vx o + 0.5 σ step kk 3 . Also, equations (26) and (27) can be used to calculate other positions and velocities in Y and Z directions, respectively. Hence, the final position and velocity via fourth order Adams-Bashforth method can be computed as [7]: Similarly, equations (28) and (29) can be used to calculate other positions and velocities in Y and Z directions, respectively.

Results and Discussion
A program in MATLAB was designed to analyze the variance in six orbital elements for a satellite in a LEO through 1300 revolutions. Figure -1 describes the semi-major axis; it is observed that, in ideal circumstances (no atmospheric drag influence), the value of the semi-major axis remains constant with time whereas it is reduced as the time passing under the effect of drag perturbation. On the 41 th day, the value of semi-major axis is 6563.5 km, but in the ideal condition it was 6584.5 km. Thus, drag compresses the satellite's orbit. According to this perturbation, the size of the orbit shrinks through 84.23 days. The effect of atmospheric drag on eccentricity is shown in Figure -2; on the 41 th day, the value of eccentricity is 0.0009965, while it was 0.001 for the ideal state. It is obvious that atmospheric drag not only contracts the size of the orbit, but also changes its shape. This leads to a breakdown in the satellite line.

Conclusions
The semi-major axis, eccentricity, and inclination have a secular dropping behavior, the Longitude of Ascending Node is periodic, the Argument of Perigee has a secular growing variation, and the true anomaly grows linearly from 0 to 360 °. Also, the semi-major axis, eccentricity, and inclination are more affected by atmospheric drag as the time passes than other orbital elements. Furthermore, these three elements have a fast dropping through the last revolutions. The behavior of the semi major axis, eccentricity, and inclination shows a good agreement as compared with several published studies, which used only the fourth order Runge-Kutta method [Ahmed et al., 2017, Mohammed andAbdul-Rahman, 2018]. The present study demonstrate that the modification in the fourth order Adams-Bshforth method gives more accurate values in orbital elements for a satellite in a LEO under atmospheric drag effect.