Numerical Model of Perturbated Earth’s Satellite Orbit

The study aimed to develop a two-dimensional numerical model of a perturbed Earth’s satellite orbit under the influence of the Moon. The first step was to model, numerically, the Earth-satellite orbit. The interaction was assumed to be first order. The basis of the model was that for two-dimensional motion, influence in the radial direction does not affect the motion in the tangential direction and vice versa. Based on this, the satellite’s motion was decomposed into radial and tangential directions. The trajectory was segmented into time intervals and the curve swept over each interval was approximated as a straight line with the assumption that acceleration in each interval was constant. Equations of constant accelerated motion were used to describe the motion of the satellite over each interval. When the model results were compared with the exact solution, for an elliptical orbit, they matched perfectly well over the entire orbit with a maximum relative error of 0.079%. When it was tested for other orbits, circular, hyperbolic, etc., it retained all of them according to theoretical predictions. The model was then extended to incorporate the effects of the Moon by launching the satellite at quarter, half and three-quarter distance from Earth to Moon. A circular orbit was used to model the effects of the Moon. The acceleration results of the model were compared with theoretical predictions. The corresponding errors in the acceleration for the three positions of launch were 0.019% and 0.20%. This showed that this model is applicable for predicting perturbated satellite orbit and it can be applied with any extra force to describe perturbated orbit of the satellite. It can also be used to model the trajectory of projectile motion, of which the exact solution is incapable of generating. Since this model gives the speed of the satellite at any instant, it can be applied when the orbit needs to be changed as it can be used to compute the required new speed.


Introduction
The forces acting on an artificial satellite are responsible for its motion as well as its deviation from the desired orbit. This deviation is referred to as perturbation. Perturbation has detrimental effects. It alters the altitude, acceleration and orientation of the satellite [1]. The disturbance in motion and underlying causes are important to look at so that they can be contained since the satellite has to maintain its desired orbit and orientation so that it can perform its task. Since the total forces acting on the satellite are not fully comprehended [2][3][4][5][6], accurate prediction of satellite position and orbit remains a continual endeavour. But also, even if known, to precisely express all of them mathematically, and provide an analytical solution is a complex task [7,8]. In the unavailability of exact solutions, numerical solutions are developed as an alternate way of predicting the location, trajectory and the general behaviour of the satellite at a given instant.
Various approaches have been used to develop numerical models of satellite orbit. They include; Euler's method [9,10], Runge-Kutta method [10,11], power series methods [12] and Lagrangian methods [13,14] among others. The Euler, Runge-Kutta and power series methods are based on, or derived from, infinite series method. It follows that their accuracy depends on the truncation, that is, the larger the number of terms the closer the method to the exact solution. But there is a trade-off, as the number of terms increase the time taken to reach the solution elongates, making them slow. At the same time, when the terms are reduced, the method becomes faster but less accurate. Lagrangian methods, on the other hand, are direct, but difficult to solve [13,14]. Thus, there is a need for a simple, faster and accurate numerical solution to model the interaction of the satellite with other celestial bodies which is what this study sought to address.

Materials and Methods
This work was divided into two parts. The first part was derivation of a two-dimensional numerical solution for the Earth-satellite orbit. This was then compared with the exact solution. The second part involved extending the numerical model to incorporate the effect of the Moon in the satellite orbit. In both cases the bodies were assumed to be coplanar. Maple program was used to compute and plot the output results.

Derivation of the Numerical Solution of the Earth-Satellite Orbit
The Earth and satellite were modelled as a system of twopoint masses, and , where the only forces were those due to interaction potential. In this case, the satellite orbits the Earth under central force motion. The numerical model of the orbit of satellite was based on the fact that, motion in two-dimensions can be modelled as two independent motions in each of the two perpendicular directions associated with radial ( ) and tangential ( ) directions. That is to say, influence in the direction does not affect the motion in the direction and vice versa [15]. In this work, only gravitational interaction between the two bodies was considered. This is because the exact solution, of a two-body system which the model was compared with, only uses the gravitational force.

Motion in the Tangential Direction
The trajectory of the satellite was segmented into time intervals and the curve swept by the satellite over each interval was approximated as a straight line with the assumption that acceleration in each interval was constant. This was more so because for an orbit of large radius, → ∞, orbital motion taken at an infinitesimal step, defined as ∆ , can be approximated as a straight line. Hence, equations of constant acceleration motion were used over each interval. Also, the velocity at which the satellite was launched was assumed to be purely tangential.
To derive the equations of motion for the tangential direction, the component of the final tangential position ( ) of the satellite at any instant was modelled using the following equation Where was the initial tangential position. The change in the tangential direction, denoted ∆ was defined using the following equation.
This was adapted from the equation of an arc ∆ = ∆ Where is the radial distance and ∆ is the infinitesimal length of the satellite trajectory. Here, ∆ approximates the arc with a straight line. The change between the final position ( ) and the initial position ( ) is given by, In this case, denote the initial angular speed, ∆ denote the change in time and , denotes the acceleration of the satellite due to the Earth. Since there is no force in the tangential direction for a two-body system, then implying that, , = 0. Equation (4) becomes, From equations of motion, the tangential velocity, is given by; = + , ∆ ⇒ = Hence the final tangential position is expressed as, = + ! ∆"

Motion in the Radial Direction
The component of the motion of satellite in the radial direction was modelled using equations of motion since acceleration is constant.
From equations of straight-line motion, the final (radial) position ( ) of the satellite is given by; Where denotes the initial radial position of the satellite measured from the centre of the Earth, # is the initial radial velocity and , denotes the radial acceleration due to force exerted by the Earth on the satellite.
The final radial velocity is given by; Acceleration of earth ( ) towards satellite was ignored in this work because The magnitude of the force exerted by the Earth on the satellite denoted, is; The negative sign shows that the force is attractive.
Substituting (13) into (12) , becomes; Substituting (14) into (10) # becomes; Since launching velocity is purely tangential, then the initial radial velocity # = 0, which implies that This shows that the satellite's final radial velocity is proportional to radial acceleration for constant ∆ . On the other hand, the radial displacement from initial position to final position after ∆ is; Since # = 0, then This shows that is less than by ∆ , that is, The results show that the radial distance decreased as the satellite moved from the initial position to the final position after ∆ . Since 9 the satellite has accumulated a radial velocity, which was initially # = 0 as it moved from position to . This implies that the radial velocity # is no longer radial to the new position but makes an angle ∆ with the new radial direction as shown in figure 1 below. The real components of the tangential ( ) and radial (#) velocities corresponding to the new radial position are computed by decomposing the magnitudes of # and to adjust them into the radial and tangential components of the velocities.
The successive positions of the satellite are a repetition of this process. The above procedure was used to calculate the rest of the successive positions for the satellite by replacing "?" and "<" by "=" and "= − 1" respectively. Maple program was used to compute the results at each iteration.

Exact Solution of a Two-Body System
The exact solution of an orbit for the gravitational interaction is expressed as; = @ 0A BC Physically, D is the radius of the circular orbit corresponding to the given values of E, F and G. That is, The dimensionless parameter N , called the eccentricity, characterises the shape of the orbit and is expressed as; Is the total energy of the system. The constant G is defined as Where 5 is the universal gravitation constant. The angular momentum, E of the system is given by, Where # is the initial tangential velocity of the satellite, that is, the launching velocity of the satellite.

Error Analysis
The method used for error analysis was the relative error method which takes the difference between the exact value and the approximate value at each iteration and then divide it with the exact value at the respective iteration. The relative error (Q R ) at each iteration was given by; The maximum of these errors was then extracted to see the furthest value from the exact in percentage form.

Extending the Model to Incorporate Effects of the Moon on the Satellite's Orbit
The essence of this work was not only comparing the numerical solution with the exact solution for a two-body system, but also to extend and use it independently in the case of three body orbit involving the Earth, Moon and satellite. As with the two-body orbit system, there were assumptions made; 1. The Moon's orbit around the Earth was assumed to be circular. 2. The Earth, Moon and the satellite were assumed to be coplanar hence the motion of the orbiting bodies was modelled using a two-dimensional approach. The approach used in this case is the same as that used in the two-body system involving the satellite and Earth with the only difference that the radial and tangential components of motion of satellite were modified to incorporate the effects of the Moon's gravity in the orbit of the satellite.
From equations of motion, the final radial position of the satellite could be modelled as; All the components of this equation are the same as those for the Earth-satellite system with the exception of \ ] which, in this case, is the sum of the accelerations due to both the Earth's and Moon's forces.
The radial acceleration due to the Earth remained the same as that for the two-body system. In order to evaluate the acceleration of the satellite due to the Moon's force, the coordinates for the Moon and satellite were transformed from plane polar to Cartesian as shown in figure 2. This was done because it was easier to work with vector products in Cartesian coordinates than in plane polar. Where In this case, f is the angular velocity of the Moon. Since the Moon is assumed to move in a circular orbit around the Earth, with centripetal acceleration # c & and c being constant, then the angular velocity (f) of the Moon can be derived from Newton's second law of motion, that is, Multiplying by 1 c & on both sides, then The satellites position in Cartesian coordinates is given by Where the value of and are as evaluated from the twobody section.
The vector (j _ ) between the satellite and the Moon is The unit vectors in the radial and tangential directions of satellite position in Cartesian coordinates are respectively given by; = T (:; ) + e ( <= ) (37) = T (− <= ) + e (:; ) With this information, the radial and tangential acceleration due to the force of the Moon can now be derived.
The force exerted by the Moon on the satellite is decomposed into tangential ( ,) ) and radial ( ,) ) components. The acceleration in the radial direction due to the force of the Moon will be derived, and that in the tangential direction will be inferred.
The force ( ) ) exerted by the Moon on the satellite is given by; Where k is the unit vector in the direction of the Moon's force and j is the magnitude of the vector between the satellite and the Moon.
On the other hand, the radial component of force of the Moon on the Satellite is; Therefore, the component of the force in the radial direction is given by the dot product between ) and , that is; Since the vector j _ is given by, From 41, the radial acceleration can be re-written as, By inference, the tangential acceleration due to the Moon's force is given by; These two quantities ( ,) and ,) ) give the acceleration of the satellite in the radial and tangential direction, correspondingly, due to the effect of the Moon. Incorporating them into the equations of motion derived in the preceding section, the equations of motion for the satellite then become; = + ∆ Where, And, = + ∆ Since the Moon's force causes the satellite to accelerate in the tangential direction, the new ∆ has this component included, that is, The velocities in the radial (#) and tangential ( ) directions, respectively, are; In this case, = + ,) ∆ (51) # = # + m ,) + , n∆ = # + ∆ Where is the total acceleration of the satellite in the radial direction.

Numerical and Exact Solution Model Results
The model results are plotted in figure 3 below for an elliptical orbit. The results are discrete because the output was generated iteratively. Using the same initial values, a model of the satellite motion was developed using the analytic solution. The results were plotted for the interval of = p0,2qr as shown in the figure 4 below. In this case, the orbit of the satellite is represented as a continuous elliptical trajectory. This is because of the domain, that is, = p0,2qr.

Comparison of the Two Solutions
When the exact and numerical solutions were superimposed, they exactly matched each other for the entire orbit as shown in figure 5 below. This is supported by the error analysis results which showed that the maximum difference between the exact and numerical solution for the entire orbit was 0.079%. This is the maximum relative error for all the iterations for the entire orbit. This shows that the numerical model developed closely matched the exact solution. The accuracy can further be improved by reducing the time step from ∆ = 1 to a lesser value.

Modelling Various Satellite Orbits Using the Model
This model can be used in many applications. These include modelling various orbits using given initial velocities, that is hyperbolic, parabolic, circular, elliptic, and projectile motion. This model can also be used to model perturbated orbits of satellite motion due to effects of the Sun, planets and Moon. In this work, only the effects of the Moon are treated and the effects of other planets or bodies can be treated in a similar manner. The perturbated orbit due to the effect of the Moon will be analysed, numerically, later in a separate section. Another application of the model, which is treated under this section is modelling projectile motion, of which, the exact solution is incapable of generating.
The figure below shows the output which is a representation of the various satellite orbits. The initial velocities were assigned based on the range of theoretical calculations so that it can be observed if the model was able to retain the expected solutions from theory. These values are found in table 2 of the appendices. These results are consistent with theoretical predictions, that is, based on the expected trajectory given the initial velocities presented in table 2 of the appendices, the model presented exactly what was expected from theoretical predictions, that is, all the predicted orbits were retained by the numerical model. The additional strength of the model is that it is able to produce projectile motion of which exact methods predict it as an imaginary (not real) solution or trajectory.
This model can also be applied when there is a need to change the satellite orbit. For instance, when the satellite needs to be removed from orbit back to Earth, that is changing the orbit from elliptical to projectile motion. In this case, the model can be used to figure out the new speed that the satellite should possess in order to follow the desired trajectory. In this case, assuming that the perigee is at v , then when the satellite is decelerated to a certain velocity, # v,w , then the point at which it will land on Earth can be precisely located. Conversely, by identifying a suitable point for landing the satellite on Earth, usually the sea, then the required velocity for the satellite at the perigee can be calculated precisely using this model. This is a safe way of removing a satellite from orbit and only numerical models can retain this solution.

Perturbated Satellite Orbit due to the Moon's Force
In the study of the effects of the Moon on the orbit of the satellite, a circular orbit was used. This is because, the circular orbit is highly unstable and any slight disturbance to the motion lead to a change in orbit, that is, an increase or decrease in velocity leads to an elliptical or parabolic orbit, respectively. This is an advantage in this situation since the corresponding results are easily appreciable due to the effects of the Moon on the satellite.
To study the effects of the Moon, the satellite was launched at three different distances between the Earth and the Moon, that is, quarter, half and three-quarter distance. This was done to pronounce the effects of the moon as depicted by the model.

Launching Radius at a Quarter Distance Between
Earth and Moon The maximum force exerted by the Moon on the satellite in the radial direction occurs when the Earth, satellite and the Moon are collinear. There are two possibilities for this to happen. It can be when the Earth is sandwiched by the satellite and the Moon or when the satellite is sandwiched by the Earth and the Moon. The relevant scenario is when the satellite is sandwiched by the Earth and the Moon. In this case, the radial acceleration of the satellite would be at its maximum since the tangential component is zero. Unless the satellite is geosynchronous, there are infinitely many occurrences of this type of configuration. For comparison purposes the initial configuration at the radius of launch was used.
From theory, when the satellite is launched at a quarter distance from the Earth, that is, = x c (53) Where c is the distance between the Moon and Earth and is the initial distance of satellite from Earth. Also, given that, the mass of Moon and that of Earth are related by; Then, the magnitude of the acceleration of the satellite due to Earth's attraction denoted is given by show that, This means, at this distance, the satellite's acceleration due to Moon's force is ≈ 0.133% that of the Earth. As a result, the satellite remains in circular orbit as if the Moon was not there.
Since the orbit being modelled was circular, and the radius of launch was = x c , where c = 3.884 Z 10 ‚ , the angular velocity of the satellite was given by And the initial tangential velocity of the satellite was; Applying the same initial conditions as those of theoretical predictions, the acceleration of the satellite, using numerical methods, was found to be 0.001366924322. Comparing this with the theoretical predictions, the corresponding error was found to be 0.019%.
The implication was that the effect of the Moon at this radius of launch and orbit was very small and did not cause a significant variation in the orbit of the satellite for a single orbit. However, for a large number of orbits, the effects of the Moon could be noticeable since the perturbation would be incremental. The plot of the model results is in complete agreement with this as shown in figure 7.

Launching Radius at Half the Distance Between
Earth and Moon At this distance, the effects of the Moon are notable. The satellite deviates from its circular orbit and the deviation is such that the radius of orbit increases with each orbit, that is, the perturbation effects of the Moon alters the altitude of the satellites as well as its acceleration. This increase in altitude continues until such a time when the satellite is in front of both the Moon and the Earth. In this case, the radial acceleration will be inwards (towards the Earth) since the satellite, Moon and the Earth will be collinear. Hence the net force will be radially inwards towards the Earth. This is the point where the satellites auto-corrects its radius of orbit by reducing it and becoming closer to the Earth.
Theoretically, when the satellite is launched at = c, it is expected that; In this case, the ratio of the accelerations is the same as that of the masses. The acceleration of the satellite due to the force of the Moon is ≈ 1% of that of the Earth which is higher than that of the first case. That explains the small deviation of the satellite from its circular orbit as shown by the model results in Figure 8. According to the model results, the ratio of the accelerations was found to be 0.01232417950. The corresponding relative error for this computation was 0.20%. This continues to show that this model is consistent with theoretical results.

Launching Radius at Three-Quarter Distance from
Earth to Moon It is expected that the effects of the Moon's force will be highest of all the three cases.
Theoretically, The acceleration due to the Moon's force at this launching distance was ≈ 11% , which is highest of the three as expected. From the model, the ratio of the accelerations was found to be 0.1109176156 and the corresponding relative error for this evaluation was 0.20% . The relatively large contribution of the Moon in the acceleration of the satellite explains the visible large perturbation in the motion of the satellite as shown in Figure 9.

Conclusions
Based on the results of the model, it can be concluded that the model developed was accurate for prediction of perturbated orbit when the Moon was used as a perturbating body. This model can be used with any additional force to describe the resultant motion.
It can also be concluded that the model is capable of predicting projectile motion of which exact methods are not able to predict.
The model can also be applied when there is a need to change an orbit of satellite. By using this program, the new speed required for the satellite to change its orbit can be determined. learnt the value of plain explanation. I hold him in high regard.
The Department of Science, including; Chemistry, Mathematics and mostly Physics for the scientific skills they have imparted in me over the duration of my study at the University of Eswatini. These skills have been applied greatly in this work.

Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.