Improvement of the Accuracy of the Perturbed Orbital Elements for LEO Satellite by Improving 4th Order Runge–Kutta’s Method

Background/objectives: To study the motion equation under all perturbations effect for Low Earth Orbit (LEO) satellite. Predicting a satellite’s orbit is an important part of mission exploration. Methodology: Using 4th order Runge–Kutta’s method this equation was integrated numerically. In this study, the accurate perturbed value of orbital elements was calculated by using sub-steps number m during one revolution, also different step numbers nnn during 400 revolutions. The predication algorithm was applied and orbital elements changing were analyzed. The satellite in LEO influences by drag more than other perturbations regardless nnn through semi-major axis and eccentricity reducing. Findings and novelty/ improvement: The results demonstrated that when m for Runge– Kutta’s method is large; the perturbed value for orbital element considers more acceptable. Furthermore, as nnn increases the step will reduce.


Introduction
The satellite orbit is determined by using perturbed equation of motion and perturbation effect to match the perturbed orbit with Kepler's orbit. Predicting a satellite's orbit is an important part of mission exploration, it is the beginning point in preparation whether a proposed mission is possible or not and how a satellite is taken in the consideration. The problem of calculating the orbits of satellites is not straight forward and unsettled due to a large number of perturbations affecting the basic of two-body motion [1][2]. The two-body motion considers the Earth is a perfectly symmetrical spherical mass that means the only force is centripetal force resulting from satellite motion [3][4]. The given orbital elements for a single point at initial time, make the estimation of the future position becomes relatively straight forward as the mean anomaly moves uniformly with time [5]. The main perturbations affecting the orbit of a satellite are as follows: the non-spherical of Earth, atmospheric drag, the gravitational pull of the Sun and Moon besides the solar radiation pressure. There are other effects of perturbations, but not take in the consideration in this study such as electromagnetic forces and the gravitational pull of near big plants like Jupiter. These perturbations depend on the types of satellite orbits.
For an example, satellites in LEO are intensely affected by the non-spherical of the Earth besides atmospheric drag [6][7]. The effect of these perturbations on a body depends on its distance from the center of the Earth. That leads to the model on which Kepler's orbit determination are not precise. It should not be supposed that perturbations are always small. For an example, disregarding the effect of the oblateness of the Earth on an artificial satellite would produce completely failed in the prediction of its state vectors as the time passes [1].

Methodology
The following steps are required to design a program and to achieve the requirements of this study.

The State Vectors Calculation
• The position vector is a function of time, r = Xi + Yj + Zk, its magnitude is given by: r • The velocity vector is also a function of time, v = v X i + v Y j + v Z k, its magnitude is given by: v=

The Orbital Elements Calculation
• The semi-major axis (a) is defined as [8]: where µ is the gravitational constant for the Earth.
• The eccentricity (e) is represented as [8]: where p r is the distance from perigee, a r is the distance from apogee.
• The inclination (i) is stated as follows [8]: Where h is angular momentum vector can be stated as [3]: Hence, x h yz zy • The Argument of Perigee (AP) is defined as [3]: Where u is the argument of the latitude can be written as [3]: f is the true anomaly can be represented as follows [7]: where E is the Eccentric anomaly can be obtained by solving Kepler's equation using Newton's-Raphson Method [9]. Table 1 explains the orbital elements for Cartosat-2B set at epoch obtained from NORD that available on Two Line Elements TLE, these are the orbital elements that used in our study.

Equation of Motion
The equation of motion with perturbation is named Cowell's equation, which can be defined as [6]: Where a p is the sum of all accelerations vectors, can be represented as: a p = a Drag + a J2 + a SRP + a 3rd-body (10) where a Drag is the atmospheric drag acceleration vector, a J2 is the non-spherical acceleration vector of the Earth, a SRP is the solar radiation pressure acceleration vector, and a 3rd-body is the third body gravity acceleration due to the Sun and Moon.

The Perturbations Types
The used perturbations in this study can be explained as below:

Atmospheric Drag
It is a perturbation force caused by impact between the particles in the Earth's upper atmosphere and surface of the satellite, atmospheric drag acceleration vector can be stated as [10]: where C D is the drag coefficient between 2 and 3, the used value is 2.2, A is the area of the satellite the used value is 5.1 m 2 , M is the mass of satellite the used value is 900 kg, ñ is the air density at satellite altitude, which depends on the satellite's altitude, that is obtained from NLRMSISE-00, and r v is the relative speed between satellite and the atmosphere.

Non-Spherical of the Earth (J 2 )
The Earth, similar to other planets in the solar system has high rotational rates, bulges out at the equator. The Earth's equatorial radius is 21 km larger than the polar radius. This flattening at the poles is called oblateness (non-spherical of the Earth), the perturbing acceleration vectors effect on a satellite are given as [8,[11][12]: Their magnitude is given by: a a a a = + + where 2 J = 0.001082, the other factors of the Earth oblateness like 3 4 5 6 J , J , J , J have small effect on the satellite orbit, (R Earth ) is the mean radius of the Earth, J2X a and J2Y a are symmetric due to the equatorial bulge.

Solar Radiation Pressure (SRP)
It is a force on the satellite caused by momentum flux from the Sun. The perturbing acceleration vector can be stated as [11]: where C R is the radiation pressure coefficient the used value 1. ) has very small variation through the orbit.

Third Body Attraction
It refers to any other body in space besides the Earth which could have a gravitational influence on the satellite. The equation that describes the perturbing acceleration vector due to the Sun can be defined as follows [11]: where µ Sun is the gravitational constant for the Sun, r Sun is the position vector from the Sun, r is the position of the satellite, r Sun is the Sun distance, and r Sun-sat. = r Sun -r sat. is the Sunsatellite distance. There is another equation describes the effect of third body attraction due to the Moon can be represented as follows [11]: where v Xo is the initial velocity at epoch, v X1 is the predicated velocity, k 1 = ax o, k 2 =  (17) and (18) can be used to calculate the other velocities and positions components by the same way.

Results and Discussion
The input data from Cartosat-2B were used to compare the results, as shown in Table 1. and orbital period during one revolution and for all perturbations; the period of this revolution is about 5317.5 second, each step represents 5.3 second. h step for Range Kutta's solution equals to 5.3/m, which was used to determine the accurate perturbed value. According to Figure 1, the semi-major axis has a periodic variation behavior and one peak, as (m) increasing leads to small amplitude, it changes between about (6584.69-6584.71) km, as compared to other (m) their changing is larger. Eccentricity has a secular variation behavior with too slight value, constant, and identical for all (m). In this figure, one can note the eccentricity has not noticeable changing in one revolution, because it is too small. Inclination and Longitude of Ascending Node have a periodic variation behavior, and two peaks, as (m) increasing makes the amplitude become too small that considers more appropriate. Argument of Perigee has a secular variation behavior, and constant, but not identical, as (m) increasing this value be more suitable as compared to the ideal value. The true anomaly, which represents the motion of the satellite with respect to the Earth; it grows frequently from (0-360)°. Figure 2 shows the orbital period (T p ) variation for the satellite between about (88.6245-88.6264) minute, as (m) increasing orbital period is more stable. Table 3 gives the orbital elements at the beginning and the end of orbital period; the variations are very small less than 0.00002 km for semi-major axis, the difference for eccentricity is not noticeable during one revolution, 0.0000001° for inclination and Longitude of Ascending Node, and 0.00000001° for argument of perigee. Figures 3-9 show the behavior of orbital elements and orbital period during 400 revolutions. The semi-major axis in Figure 3 has a secular behavior. It reduces about 6 km for nnn = 1800; it's dropping too slowly as compared to small (nnn). The reducing of the orbit continues as the time passes. From Figure 4, one can note the eccentricity has also a secular behavior with too slight values, for nnn = 1800 the eccentricity has a stable and slight dropping with noticeable changing as compared with one revolution that produces a shrinking in the orbit shape. Inclination in Figure 5 has a periodic variation behavior, and one peak, in case of nnn = 1800, the amplitude becomes too small. It is changed about 0.0012°, while others have a variation more than 0.0012° they are considered as unacceptable variation. According to Figure 6, Longitude of Ascending Node has a secular variation behavior, the dropping is too slow for nnn = 1800. It changes about 0.005°. The Argument of Perigee, which has a secular variation behavior, it grows suddenly from  (120-131)°, as shown in Figure 7. In general, the behavior of perturbed orbital elements becomes more stable as nnn increases. The true anomaly in Figure 8 grows from (0-360)°. Finally, Figure 9 represents the variation of satellite's orbital period. That has a secular variation behavior when nnn increasing leads to the satellite's dropping too slowly, it reduces about (0.1) minute as compared to others.

Conclusions
From this research, we can conclude the following: 1. The results clearly showed the satellite at LEO is more affected by the drag more than other perturbations through the decaying of semi-major axis and the eccentricity that described the shape and the size of the orbit. The variation of orbital elements affected on the coordinates system of the satellite, which in turn led to lose the satellite. 2. As (m) of Runge-Kutta's solution for one revolution became too large the value of orbital element was more stable. The semi-major axis, inclination, longitude of ascending node and time period behave periodically through the time passing, but their magnitudes were different. Eccentricity and argument of perigee behaved secularly with time and also their magnitudes were different, while true anomaly grows from (0-360)°. 4. By (nnn) increasing the step it would reduce and the value of perturbed orbital elements which became more stable. Besides, the lifetime of the satellite, the semi-major axis, eccentricity, longitude of ascending node, argument of perigee, and time period behaved secularly with time while inclination behaved periodically as the time passing. The true anomaly remains in its path from (0-360)°.