Dynamics and Control of Typical Orbits around Saturn

: This paper investigates the dynamics of some typical orbits around Saturn, including sun-synchronous orbits, repeating ground track orbits, frozen orbits, and stationary orbits, and corresponding control methods mainly based on the mean element theory. The leading terms of Saturn’s aspheric gravitational ﬁeld, J 2 and J 4 terms, are used when designing the orbits around Saturn. Two control methods of sun-synchronous orbits, including initial inclination-biased method and periodic inclination-biased method, are used to damp the local time drift at the descending node, which is caused by solar gravitation and atmospheric drag. The compensation of semimajor axis and maneuver period to maintain the recursive feature of repeating ground orbits are calculated. While only J 2 and J 3 terms are taken into account, we examine the argument that the perigee of frozen orbits around Saturn should be 270 deg to promise meaningful eccentricity. The perturbations of inclination and eccentricity of stationary orbits due to solar gravitation and solar radiation pressure are presented. Meanwhile, the preliminary control strategies of inclination perturbation and eccentricity perturbation are naturally introduced.


Introduction
Interplanetary exploration is one of the most important methods to seek the origin and evolution of the universe, suitable living planets for human beings, and the existence of other intelligent life. Over the last few decades, with the development of science and technology, we have launched numerous spaceships to explore the planets of our solar system. Saturn, since its discovery, is absolutely one of the most attractive planets in the solar system. The tremendous thin ring consisting of many ringlets has become the famous characteristic of Saturn. Saturn is the second biggest planet in the solar system with a volume about 755 times larger than Earth. It is made predominantly of hydrogen and helium, which causes its density to be smaller than water. Saturn rotates so fast that its rotation period is only 10.656 h [1]. However, its orbital period is tremendously longer than its rotation period, at 29.4 Earth years [1]. The gravitational field of the Saturn system is extremely complicated, as Saturn is the planet with the largest amount of moons in the solar system. There are 53 known moons and 29 moons waiting to be formally confirmed [2]. The largest moon in the Saturn system, named Titan, was found to have a nitrogen-rich atmosphere similar to that of ancient Earth.
Therefore, the exploration of the Saturnian system could promote the process of searching for habitable planets and researching the evolution of planets in the solar system. However, few spacecraft had investigated Saturn. Pioneer11 is the first spacecraft to reach the Saturn system. It reached Saturn in 1979, and discovered a narrow ring, which is named the F-ring, outside the A-ring for the first time [3]. It also first discovered the existence of the magnetosphere of Saturn. After two years, Voyager 2, the third spacecraft to visit Saturn, gained a glimpse of Saturn on its way to Uranus. Using photopolarimeter, we have chosen a reasonable repetition parameter, we could calculate the decay rate of orbital semimajor axis with an initial condition of semimajor axis and eccentricity. To maintain the recursive feature, we give a solution of semimajor axis compensation and control period. When we use J 2 and J 4 terms to design frozen orbits around Saturn, we naturally see a novel result that the argument of perigee of frozen orbits around Saturn must be 270 • to make eccentricity meaningful. At last, we give some strategies to offset the influence on inclination caused by the solar gravitation and the influence on eccentricity caused by solar radiation pressure.
In this paper, we use some traditional symbols to denote the basic parameters of Saturn and the orbital elements of spacecraft motion. a is the semi-major axis of orbit. e is the orbital eccentricity. i is the orbital inclination. J 2 is the second-order zonal harmonic. J 3 is the third-order zonal harmonic. J 4 is the fourth-order zonal harmonic. J 22 is the second-degree and -order tesseral harmonic. M is the mean anomaly. f is the true anomaly. n is the mean angular velocity. n s is the mean angular velocity of Saturn around the Sun. p is the semiparameter. R e is the reference equatorial radius of Saturn. r is the position of the spacecraft. Ω is the right ascension of the ascending node. λ is the latitude of body-fixed coordinate system. µ is the Saturn's gravitational constant. ϕ is the longitude of body-fixed coordinate system. ω is the argument of perigee. ω s is the rotational angular velocity of Saturn.
In this paper, we mainly use the aspheric gravitational field zonal harmonic coefficients to design different typical orbits. Some basic parameters of Saturn are listed in Table 1 [1, 7,30]. Similar to Jupiter, the magnitude of J 2 of Saturn is not predominant among the zonal harmonic coefficients [28,29]. Therefore, we choose the dynamic model involved J 2 and J 4 terms. The gravitational potential function could be approximated as [31]:

Sun-Synchronous Orbits and Orbital Maintenance
The first order secular term for the precession rate of ascending node is [15,16,32]: Here, the second order secular perturbation considered J 4 for the precession rate is [16] Sun-synchronous orbits require the precession rate of ascending node to be equal to the Saturn rotational angular speed around the Sun. Hence, the average nodal precession rate is Here n s is the mean angular velocity around the Sun. Thus, we gained the equation from Equation (3): where For sun-synchronous orbits, we can derive the inequality systems (5) to calculate the inclination from theoretical analysis. The following inequalities are the requirements to find the roots of˙Ω = n s .
The first inequality in (5) prevents the spacecraft from crashing into the outer atmosphere of Saturn. The second inequality is the discriminant of Equation (4). Then we obtain the plot of the discriminant, Figure 1. We learn from Figure 1 that the discriminant is almost always positive for any chosen a and e. It means Equation (4) only has one real root, and there is only one sun-synchronous orbit for any chosen semimajor axis and eccentricity. Here, we choose the orbital elements a = 62,268 km and e = 0.01, substitute them into Equation (4), and solve the equation. Then the inclination i = 90.0483 deg can be obtained naturally.
Here we continue to analyze the perturbations of spacecraft in circular sun-synchronous orbits. The spacecraft located in circular sun-synchronous orbits around Saturn would be perturbed by the Sun. The derivative of inclination is presented by using the Lagrange equations as [33,34]: 3n 2 s cos(ω + f ) cos ξ n (cos β s sin i sin Ω − sin β s cos i s sin i cos Ω + sinβ s sin i s cos i), where F n is the normal perturbation force caused by the solar gravitation, n s is the angular speed of Saturn around the Sun, β s is the ecliptic longitude of the Sun, i s is the obliquity of the ecliptic of Saturn, and ξ is the angle between the unit vector from center of mass of Saturn pointing the spacecraft and the unit vector from center of mass of Saturn pointing at the Sun. Then we average the secular derivative of inclination in a period of circular sun-synchronous orbits [35]. We could see the secular inclination perturbation as: It is important to emphasize that we only consider the secular effect of perturbations.
We can substitute di dt for di dt for convenience. (2β s − 2Ω) in Equation (6) is associated with the local time drift at the descending node.
The semimajor axis decay due to atmospheric drag leads to a shorter period and damage of the recursive feature. Here we simply assume the atmosphere of Saturn is stationary. According to the preliminary study of Saturn atmosphere model [8], we could estimate the acceleration of semi-major axis caused by atmospheric drag. We describe the variation of semimajor axis due to atmospheric drag of circular orbits by [29,36] where C d is the drag coefficient, S is the projected area, ρ is the neutral atmosphere density, and m is the mass of spacecraft. When we analyze the inclination perturbation, we can calculate the local time drift at the descending node instead. Meanwhile, it is a fact that the local time drift at the descending node is also the local time drift at ascending node. Thus, we could continue to analyze the evolution of the ascending node. According to Equation (3), the first-order approximation of˙Ω by using first-order Taylor expansion is The unit of spacecraft running time is second. The local time drift is counted in minutes. When we only investigate the local time drift caused by solar gravitation, Equation (9) could be rewritten as There are four initial variables, ∆Ω(t 0 ), ∆a(t 0 ), ∆i(t 0 ), and −(β s − Ω) in Equation (10). In this paper, we assume ∆Ω(t 0 ) = 0 and only consider the effect of other three variables of a given sun-synchronous orbit around Saturn.
The initial orbital elements of the circular sun-synchronous orbits given in this where the local time drift will reach the extremum after 10 Earth years. Therefore, in this paper, we fix the value −(β s − Ω) = −135 deg, which means the local time at the descending node is 15:00.
When we fix the initial inclination deviation and choose different initial semimajor axis deviation, comparing Figure 2b with Figure 2d, the extremum of local time drift of these two figures almost have the same distribution. This means the first-order approximation of local time drift is sensitive to the initial inclination deviation and insensitive to the initial semimajor axis deviation. Thus, we fix the semimajor axis deviation in this paper ∆a 0 = 1 km.
When we consider the solar gravitation perturbation and the atmospheric drag at the same time, we should use Equation (9) to calculate the local time drift at the descending node. The fact we should notice is that the changing trend of da dt and di dt is quadratic function of time t [29]. Therefore, we try to simulate the effect of atmospheric drag by using the inclinaiton drift caused by solar gravity. In this way, we only need to choose inclination prebiased methods to offset the local time drift caused by two different kinds of force. The explicit calculation of local time drift: We find two figures to test the effect of simulating method used in Equation (11). In each figure, there are two curves representing the local time drift in 5 years in different calculation methods. The curves which have considered atmospheric drag use Equation (9). The curves simulating the atmospheric drag use Equation (11). Comparing Figure 3a with Figure 3b, we learn that the circular sun-synchronous orbits with the smaller initial inclination deviation has the better simulating effect of atmospheric drag.
Then we use two methods in this paper to offset the local time drift at the descending node.
The first control strategy is using initial inclination biased method. This control strategy suits the spacecraft whose entire lifetime is not too long. We take an initial inclination bias to make the spacecraft located in a quasi-sun-synchronous orbit. Then the inclination of the quasi-sun-synchronous orbit will oscillate around its normal inclination in its lifetime but still maintain the feature of sun-synchronous orbits. When t e = − ∆i(t 0 ) di dt , the local time drift ∆T in Equation (11) sees its extremum If the end of lifetime of the spacecraft is t f , then the limit of local time drift at descending node is |∆T(t f )| ≤ ∆T(t e ). Let ∆T(t f ) = −∆T(t e ), then we find the initial inclination bias: From Figure 3c,d, we learn that the prebiased initial strategy damps the local time drift at the descending node in the designed 5 Earth years' lifetime of spacecraft.
The second control strategy is designed for the spacecraft which need to execute a long-time mission. This method needs periodic inclination bias to damp the local time drift caused by solar gravitation. The extremum of local time drift Equation (11) in a period and the value of periodic inclination bias can be derived naturally: where the sign of ∆i(t c 0 ) is opposite to the sign of di dt . From Figure 3c,d, compared with Figure 4a,b, the periodic inclination biased method has the better effect of local time drift limitation. Meanwhile, comparing Figure 4a with Figure 4b, it is obviously also sensitive to the initial inclination deviation. Therefore, regardless of the methods, it is essential to limit the initial inclination deviation ∆i 0 .

Repeating Ground Track Orbits and Orbital Maintenance
The repeating ground track is significant for remote sensing satellites. The condition for spacecraft achieving repeating ground track is R∆λ = 2πN, where R and N are relatively positive prime numbers. The interval of the adjacent ground track on the equator is [28,29] which ω S is the angular rotational velocity of Saturn. The nodal period of the motion of spacecraft is According to the mean element theory [16], the average perturbing rates of M and ω is [37]˙M = n +Ṁ 1 +Ṁ 2 ,ω =ω 1 +ω 2 , In this paper, we use the ground track repetition parameter Q = R/N to describe different repeating ground track orbits. Here, we should pay attention to the meaningful range of repetition parameter Q [29].
To find a suitable repetition parameter, we should firstly find the lower bound of Q by using the discriminant of f (sin i): where e 2 ) + n + Q(n s − ω s ).
To promise f (sin i) has four real roots, the discriminant should be non-negative, i.e., Solving the inequality, we find the lower bound of Q: To estimate the maximum of repetition parameter, we use the rotational period of Saturn T s and the period of spacecraft in a low orbit T l . Then the upper bound of Q can be presented as: Using the data from the planet model of Saturn, we find the meaningful range of repetition parameter for engineering application: 2.4701 ≤ Q ≤ 2.5191.
For different chosen values of repetition parameter Q, we see a function f (sin i) that describes the relation between the inclination i and the semimajor axis a of repeating ground track orbits. The relation between i and a are shown in Figure 5 for three given meaningful values of Q. Using Equation (21), when initial semimajor axis and eccentricity are given, the corresponding inclination is the root of f (sin i). In Figure 6, we present the repeating ground track orbit around Saturn when initial orbital elements and repetition parameter Q are given.  An inevitable perturbation which can cause the semimajor decay in repeating ground track orbits is atmospheric drag. The semimajor axis decay due to atmospheric drag leads to a shorter period and damage of the recursive feature.
The difference between actual angular velocity and nominal angular velocity can be presented as [29]: where a 0 is the initial semimajor axis,ā is the nominal semimajor axis, andn is nominal angular speed. The increase in angular speed results in the longitude drift of ground track. Obviously, we can learn from Equation (7) thatȧ < 0. Therefore, if a 0 ≤ā, the longitude will drift eastward monotonously. To maintain the feature of repeating ground track orbits, we should take a maneuver to make initial semimajor axis a 0 >ā. Then, in the first coming time interval t s = − ∆ȧ a , the longitude will drift westward. When d∆λ dt = − 3π a (∆a +ȧt) = 0, the longitude of the ground track would reach the western boundary. In the second time interval t s , the longitude will drift eastward and finally return back to the initial longitude. Then we can find the longitude drift compared to the nominal repeating ground track orbits: where ω s is the angular rotation speed of Saturnian.
In the whole control period, when we fix the limitation of longitude drift ∆λ L , then we can find the control period ∆T C and compensation ∆a of semimajor axis [29]: Thus, in each period, when the longitude of ground track drift to the eastern boundary, we have to take a maneuver ∆a to compensate the decay of semimajor axis. Here we assumed projected area of spacecraft S = 20 m 2 , C d = 2.1, ∆λ L = 10 km, and mass m = 3000 kg. The density of simplified atmospheric model, in the range of semi-major axis a ∈ [62,268, 62,468] km, is ρ ∈ [4.7 × 10 −12 , 3.7 × 10 −12 ] kg/m 3 [8]. Then we find the Figures 7 and 8 which illustrate the the control compensation of semimajor axis and control period of repeating ground track orbits. From Figure 7, we know the semimajor compensation ∆a for a ∈ [62,268, 62,468] km varies from 5200 m to 5800 km. The corresponding compensation period ∆T C in Figure 8 varies from 18 h to 16 h. Compared to Jupiter [29], Saturnian atmospheric drag is about 2 orders of magnitude greater than Jovian atmospheric drag when the spacecraft is at the same altitude. As a result, it is essential to execute the semimajor axis compensation for repeating ground track orbits around Saturn.

Frozen Orbits
Before analyzing frozen orbits, we first investigate the orbits at the critical inclination. The traditional design of the orbits at critical inclination around Earth just considers the J 2 term by which the first order secular variation of eccentricity and argument of perigee are [15,16]:ė Different from the methods which only considered J 2 , we still need to calculate the perturbation caused by J 2 and J 4 terms. Though the magnitude of J 2 term exceeds J 3 term 1000 times, J 4 term still plays an important role in aspheric perturbation. Thus, the conditions for the orbits at the critical inclination should be combined ω 1 with ω 2 : However, the critical inclination cannot satisfy most missions for remote sensing satellites. For the purpose of missions, we must consider the frozen orbits around Saturn. Frozen orbits are the orbits which own small constant eccentricity and constant argument of perigee while arbitrary inclination and orbiting altitude can be chosen. The frozen orbits we designed which have small eccentricity must consider first order long-period terms [23]. The main parts of long-period terms are [16] Then the complete conditions for frozen orbits are combined by Equations (24) and (25) [23].
Using ω = 90 • or ω = 270 • to solve Equation (26), respectively, we learn that only when ω = 270 • can the eccentricity e be positive. In Figure 9, the values of eccentricity are presented for different combinations of semimajor axis and inclianition.

Stationary Orbits and Orbital Maintenance
Stationary orbits are more complex than those we have considered before. Therefore, we investigate stationary orbits in a special condition. When we investigate the stationary orbits around terrestrial planets, we must take the J 22 term into our gravitational field model and analyze the longitude drift caused by aspheric perturbation. Here we only analyze aspheric perturbation caused by the J 2 , J 4 terms, as Saturn is exactly a gaseous planet which means the J 22 term could be neglected.
The spherical coordinates O − r, λ, ϕ, where O is the centroid of Saturn, r is the distant from the spacecraft instant position to O, λ, and ϕ represent the longitude and latitude of the spacecraft separately. The equations of motion of stationary orbits in the spherical coordinates can be written as [28] 2r 4 (3 sin 2 ϕ − 1) + 5µJ 4 R 4 e 8r 6 (35 sin 4 ϕ − 30 sin 2 ϕ + 3), d dt (r 2λ cos 2 ϕ) =0, Spacecraft staying in stationary orbits around Saturn must satisfy these conditions [28]: Then we gain the equation of stationary orbits at ϕ = 0 as Solving Equation (28), we find the radius of stationary orbits around Saturn is r s = 112,506.0294 km.
After learning the structure of Saturn's rings [12], we know that spacecraft in the stationary orbit around Saturn may have a collision with the B-ring.
Using nonsingular elements {a, e x = e cos(Ω + ω), e y = e sin(Ω + ω), i x = sin i sin Ω, i y = sin i cos Ω, λ = Ω + ω + M}, we first calculate the secular solar gravitation perturbation of i x , i y and e x , e y . When ϕ = 0, the normal perturbation force caused by the solar gravitation can be written as [33]: F n = 3rn 2 s (sin α sin β s cos i s + cos α cos β s ) sin β s sin i s , where α is the longitude of spacecraft in the stationary orbit, n s is the angular speed of Saturn around the Sun, β s is the ecliptic longitude of the Sun, i s is the obliquity of the ecliptic of Saturn. The derivative of i x , i y perturbed by normal solar gravitation in a Saturnian rotation period are [33] di x dt = 3n 2 s n [ 1 4 sin 2 β s sin 2i s + ( 1 8 sin 2β s cos i s )i x , Here, we apply the Equations (30) and (31) to analyze the perturbation of the stationary orbit caused by solar gravitation, then we find the average rate of inclination vector perturbed by solar gravity in a sidereal orbit period: After solving this first order linear ordinary differential equation system Equation (32), we find the i x (t), i y (t) for any time t with initial condition i x (t 0 ), i y (t 0 ).
Before we analyze the perturbation of inclination vector, we firstly set the maximum of inclination perturbation |∆i max | = 0.2865 • . Then we could draw the Figure 10a,b of i x , i y separately for each initial condition value i x 0 = i x (t 0 ), i y 0 = i y (t 0 ) in the limited circle. Once we see Figure 10a,b, what we know immediately is that the inclination perturbation is so small that the magnitude is smaller than 1 × 10 −3 in 5 Earth years. Furthermore, comparing Figure 10a with Figure 10b, we learn evidently that the magnitude of ∆i y is two order magnitude smaller than ∆i x in a given time of 5 Earth years. This means when we plot the perturbation of inclination vector of stationary orbit around Saturn over a long time, we could hardly see a sight difference in i y . Though the inclination perturbation caused by solar gravity could hardly influence the spacecraft in stationary orbit around Saturn, we still analyze the control strategy of i x , i y in a preliminary way. When we know how the solar gravity perturbation acts, we could calculate ∆i for any initial (i x 0 , i y 0 ) to stay in the limited circle. The area of the cool tone circle overlap limited circle in Figure 11 indicates the initial point (i x 0 , i y 0 ) in this area are still remained in the limited circle after 100 Earth days. By observing the perturbation feature of stationary orbit around Saturn, we know that minimum maneuver − → ∆i to keep i x , i y staying in the limited circle is from the boundary of limited circle pointing towards the center of the cool tone circle. The value of − → ∆i of (i x 0 , i y 0 ) in warm tone area in Figure 11 indicates the minimum inclination maneuver of (i x 0 , i y 0 ) to keep (i x (t), i y (t)) remained in the limited circle after 100 Earth days. Figure 11. ∆i for control period ∆T C = 100 Earth days with initial inclination vector point (i x 0 , i y 0 ) in the inclination limited circle whose radius r = 4.58 × 10 −5 deg perturbed by solar gravity.
After analyzing the effect caused by solar gravity perturbation, we now turn to the solar radiation pressure which may lead to eccentricity drift. First, the mean solar irradiance around Saturn is P = 15.04 W/m 2 while the distance from the Sun to Saturn is 1.426 × 10 12 m [38]. Then we find the intensity of solar radiation pressure is p = P c = 5.0168 × 10 −8 N/m 2 , where c = 299,792,458 m/s is the speed of light. Using other supposed parameters: the reflection parameter K = 1, area of spacecraft vertical to the Sun A = 20 m 2 , the mass of spacecraft m = 3000 kg, then the solar radiation pressure of spacecraft on stationary orbits around Saturn is [34]: where S is the unit vector from the centroid of Saturn pointing at the Sun. Using the Lagrange perturbation equations, we could find the average time derivatives of the eccentricity vector e x = e cos(ω + Ω), e y = e sin(ω + Ω) in Equation (35).
where F r and F t are the radial component and tangential component of solar radiation pressure F s . When we analyze the eccentricity perturbation of spacecraft in stationary orbits over a long time period, we learn that the motion of the Sun in a time period could be described as ∆l s = l s (t) − l s (t 0 ) = n s (t − t 0 ). Then we find the time integral of the eccentricity vector for any given t.
We learn from Equation (37) that the curve of the eccentricity vector is an ellipse due to the existence of cos i s , and initial condition (e x 0 , e y 0 ) and l s 0 influence the position of this eccentricity ellipse. From Figure 12, we learn that l s 0 , which means the angular between initial solar radiation direction and x axis, can influence the perturbing direction. The centre of eccentricity limitation circle is (e x 0 , e y 0 ).

Figure 12.
A total of 72 curves of f (e x , e y ) perturbed by solar radiation pressure for l s0 = 5n deg, the integral n ∈ [0, 72], in a Saturnian year with initial eccentricity vector (e x 0 , e y 0 ) = (0, 0) in the eccentricity perturbation limited circle whose radius r = 3 × 10 −7 .
Here, we fix the (e x 0 , e y 0 ) = (0, 0). When (e x 0 , e y 0 ) = (0, 0), we can first take a maneuver to make the point (e x 0 , e y 0 ) = (0, 0) in the curve of eccentricity perturbation. Then, we can design a maintenance strategy. When taking an eccentricity control, we used to set ∆V t > 0, then we could suppose B = − ∆V t T , and find the maneuver function [39]: Combined with solar radiation pressure perturbation Equation (35), the rate of eccentricity vector (e x , e y ) after a maneuver could be rewritten as We still integrate Equation (38), then we could find the eccentricity vector at any time t after we take an eccentricity control maneuver.
To find a better control effect, we take 6 times eccentricity control for different initial solar radiation direction. Compared Figure 12 with Figure 13, we could know that e x , e y are kept staying in the eccentricity limited circle after a Saturn's year for initial (e x 0 , e y 0 ) = (0, 0) and any given solar radiation direction l s 0 . Figure 13. After taking 6 times eccentricity control when B = −5 × 10 −11 m/s 2 , 72 curves of f (e x , e y ) perturbed by solar radiation pressure for l s 0 = 5n deg, integral n ∈ [0, 72] in a Saturnian year with initial eccentricity vector point (e x 0 , e y 0 ) = (0, 0) in the eccentricity perturbation limited circle whose radius r = 3 × 10 −7 .

Conclusions
In our paper, we analyze the sun-synchronous orbits, repeating ground track orbits, frozen orbits, and stationary orbits around Saturn and corresponding control strategies based on the mean element theory while zonal harmonic coefficients J 2 and J 4 of Saturnian gravitational field are considered.
For sun-synchronous orbits, we analyze the existence of orbits, find the relation between inclination with eccentricity and semi major axis, and calculate the local time drift caused by solar gravitation and atmospheric drag. After that, we take two inclinationbiased control strategies to damp the local time drift. The initial inclination-biased method suits the short-time spacecraft. The periodic inclination biased method has a better effect on long-period missions.
For repeating ground track orbits, we find the meaningful range of repetition parameter Q. Then we find the relation between the inclination and semimajor axis. After that, we calculate the compensation for semimajor axis and maneuver period.
For frozen orbits, we learn that only when ω = 270 deg can we be sure of the eccentricity positive for any given inclination and semimajor axis.
For stationary orbits, we first calculate the radius using the conditions of equilibrium point. Then, we analyze the perturbations caused by solar gravitation and solar radiation pressure. Finally, we take corresponding maneuver strategies to control the inclination and eccentricity.