Re-Evaluation Solution Methods for Kepler's Equation of an Elliptical Orbit

An evaluation was achieved by designing a matlab program to solve Kepler’s equation of an elliptical orbit for methods (Newton-Raphson, Danby, Halley and Mikkola). This involves calculating the Eccentric anomaly (E) from mean anomaly (M=0°-360°) for each step and for different values of eccentricities (e=0.1, 0.3, 0.5, 0.7 and 0.9). The results of E were demonstrated that Newton’sRaphson Danby’s, Halley’s can be used for e between (0-1). Mikkola’s method can be used for e between (0-0.6).The term that added to Danby’s method to obtain the solution of Kepler’s equation is not influence too much on the value of E. The most appropriate initial Gauss value was also determined to be (En=M), this initial value gave a good result for (E) for these methods regardless the value of e to increasing the accuracy of E. After that the orbital elements converting into state vectors within one orbital period within time 50 second, the results demonstrated that all these four methods can be used in semi-circular orbit, but in case of elliptical orbit Danby’s and Halley’s method use only for e ≤ 0.7, Mikkola’s method for e ≤ 0.01 while Newton-Raphson uses for e < 1, which considers more applicable than others to use in semi-circular and elliptical orbit. The results gave a good agreement as compared with the state vectors of Cartosat-2B satellite that available on Two Line Element (TLE).


Introduction
In celestial mechanics, one is concerned with the motions of celestial bodies under the effect of mutual mass attraction. The simplest form is the motion of two-body problem. In this case, there are two molds that must be taken in the consideration; the first is the bodies are spherically symmetric in order to treat the bodies according to their masses that determined at their centers. The second is assuming there are no outward or inward forces acting on the system other than the gravitational forces which act along the line linking the two bodies [1,2]. The development of the two-body problem, which can describe the motion of a satellite the the Earth, originates with Newton's Law of universal gravitation [3,4].
For artificial satellites, the mass of the satellite usually is ignoring, because it is small as compared with the mass of the central body. The two point masses rotate with respect to internal reference frame. An inertial reference frame moving at constant velocity and without rotation with respect to an inertial frame in which the universe seems spherically symmetric. Generally, it is indicated to the location and the orientation of the satellite coordinate axes [3,5]. A body border is defined by an origin at a specified point in the body and three Cartesian axes. It is used to arrange in a line the various components that will generally shift because of the large forces experienced during launch. Every power is made to border these motions, but they cannot always be ignored. Whether they are ignored or not depends on the pointing accuracy required of the satellite [6]. To describe an orbit, six parameters are required, which are: eccentricity (e) describes the shape, semi-major axis (a) describes the size, longitude of ascending node (Ω), argument of perigee (ω), inclination (i) describes the orientation of an orbit and the true anomaly (f), which gives the position of the satellite in its orbit at a particular time [7,8].

The Methods for Solving Kepler's Equation:
Kepler's equation is used for describing the position and velocity for a different time that affected by the gravity force and determining the relationship between time and angle that places the body in an elliptical orbit as shown below [9]: Where: (E) is the Eccentric anomaly, (M) is the Mean anomaly at a given instant of time and ( ) is the eccentricity measured in degrees equal to . In General e and M are given, and then the equation must be solved directly for E. In some cases, E and e are determined and M is unknown; in this case, the equation above cannot be used directly [9]. Mean anomaly at the time (t) can be calculated by [6]: ( ) Where: (t p ) is the time for an object to pass through the perigee, (n) is the mean motion calculated by [6]: μ Earth =GM Earth =398600 Km 3 / Sec 2 . The time for the satellite to go once around its orbit is called the period, and calculated by third Kepler's law [7]: or A number of researchers have advanced iterative and non-iterative numerical methods. There are many methods for assessing E, which is depending on its unknown value, a formula which gives approximation results, as listed below [9,10]:  Newton's-Raphson method: The idea of this method is to approximate the nonlinear function by the first two terms in a Taylor series expansion. It is defined as [9]: Where: , . To solve the above equation the initial value of Eccentric anomaly (E n equal to M), then the repetition is continued until the absolute value of (-) be smaller than the accuracy, this accuracy is chosen by the user to correspond the desired precision in calculation [10]. We can continue with term .
Where: is the fourth derivative. So eq.(8) can be written as the following:  Mikkola's method: presented a formula for initial Guess as below [10]: Where: , √ , , .
 Halley's method: the formula for initial Guess as below [10]:

The State Vectors Calculation:
Exact orbit determination is a major requirement in the case of several satellite missions, which is specified by the satellite's state at any time, can be completely described by six independent parameters in the equatorial plane, the position (X, Y, Z) and the velocity (v X , v Y , v Z ) [11,12,13]. The Cartesian coordinates (x, y, and z) for satellite in its orbit, are given by [ [14]: √ To transformation of the position and velocity components from the object plane to the equatorial plane could be achieved by using a Gaussian vector (conversion matrix), which is given by [14]: Where: R -1 is the inverse of Gauss matrix, which contents Eular angles (i, ω and Ω), as below [14]:

[ ]
Where: Eq.(18) and Eq.(19) can be written as the following: Where: X, Y, Z, X , Y, and Z represent the state vectors (position and velocity components) of the satellite in the equatorial plane at time (t).

Results and Discussion
Kepler's equation was solved by using some different initial Guess values, which depend on M and e. A program was designed by using Matlab software. These initial Guess values for E n were considered as a base value. As a result, those refined initial Guess values were used another time to solve Kepler's equation to obtain the more accurate value of E and to select the more appropriate initial value, which was calculated for each step from M. The using program for the four methods will stop according to the input accuracy 10 -10 , which was selected by the user. The relationship between M and E was plotted for several values of e, as Figures-(1-5). Figure-(1) represents the values of E at the beginning of orbital period (perigee M=0°), middle (apogee M=180°) and at the end of period (perigee M=360°), the value of E and M are close to each other, and as e converges to zero the behavior of E for the four methods approaches to linear relationship. Figures-(2 and 3) indicate to as e starts to increase, these methods begin to shift from each other with an approximately the same values of E. This value decreases at the beginning of orbital period, then be approximately the same at the middle, then increases until complete 360°. Figs (4 and 5) represent the value of E for e ≥ 0.7; these values start to be different from each other. Danby's, Mikkola's and Halley's method will not complete 360° like Newton's-Raphson. The term that added to Danby's method to obtain the solution of Kepler's equation is not influence too much on the value of E.
The most appropriate initial Gauss value was also tested and selected to be (E n =M), this initial value gave a good result for E for our using methods regardless the value of eccentricity, as Figures-(6  and 7).
The satellite orbit determination involves a set of method, which measure the satellite motion in terms of position and velocity without perturbations effect; therefore these four methods for calculating the values E were used to calculate the state vectors of the satellite (X,Y,Z, R, X , Y , Z , ). The input parameters were selected from Cartosat-2B satellite that available on TLE (i=97.9448°, Ω=207.1202°, e=0.0016257, ω=44.4835°, M= 315.7690°, h p = 622 km), and the state vectors were (X= -6231.7560 km,Y= -3189.4018 km, Z=14.8069 km, R=7000.5204 km, X = -0.4537 km/s, Y = 0.9400 km/s, Z = 7.4776 km/s, V= 7.5500 km/s). Our results by the four methods of solution gave the same values (X= -6234.3849 km,Y= -3190.7472km, Z= 14.8132 km, R= 7003.4736 km, X = -0.4536 km/s, Y = 0.9398 km/s, Z = 7.4760 km/s, = 7.5485 km/s). This refers to the value of e was too small and is not affect the values of position and velocity components this happens in semi-circular orbit as figs (8 and 9). Fig (8) represents the shortest distance of the satellite from perigee at the beginning of the orbital period and become larger at apogee, the end of the orbital period. Fig (9) represents the largest velocity of the satellite from perigee and became smaller at apogee. These two figures represent the distances and velocities that are identical for too small e (semi-circular orbit). As e begins to increase the values of the state vectors were different as figs (10)(11)(12)(13). The value of M was calculated from mean motion that changes until the orbital period of the satellite completed; the step of time was selected by the user to be about 50 seconds. Figs (10 and 11) represent the values of distances and velocities components are still identical to each other by Danby's, Newton's Raphson and Halley's method, but the value of E by Mikkola's started to be different (fluctuated). As e more increases the same behavior was occurred for e between (0.01-0.6) and Mikkola's was more fluctuated. For e ≥ 0.8 Newton's-Raphson starts to shift from Danby's and Halley's method, as Figures (12 and 13). That means the solution for Kepler's equation of an elliptical orbit by Newton's-Raphson is the best solution, but in case of semi-circular orbit all these methods can be used.

Conclusions
1. For a small value of e the E is close to M and as e increases the value of E will shift from M and the difference between E and M increases with e. 2. Newton's-Raphson, Danby's and Halley's method gave a good result of E for e between (0-1). Mikkola's method is used for e between (0-0.6). 3. The term that added to Danby's method to solve Kepler's equation is not influence too much on the value of E. 4. Multi initial Gauss values were tested and the most appropriate E n selected to be (E n =M), this initial value gave a good result for E for our using methods regardless of the value of eccentricity. 5. In case of calculating the state vectors for semi-circular orbit all these methods can be used, but for elliptical orbit, Newton's-Raphson more appropriate than others and still applicable for e between (0-1). Danby's and Halley's use for e ≥ 0.7. Mikkola's is used for e ≤ 0.01.