Journal of Astronomy and Space Sciences

The Korean Space Science Society

2020; 37(1):1-9

pISSN: 2093-5587, eISSN: 2093-1409

DOI: https://doi.org/10.5140/JASS.2020.37.1.1

Research Paper

Frozen Orbits Construction for a Lunar Solar Sail

Elamira Hend Khattab

, Mohamed Radwan

, Walid Ali Rahoma†

Department of Astronomy and Space Science, Faculty of Science, Cairo University, Cairo 12613, Egypt

^†Corresponding Author Tel: +20-2-3567-6843, E-mail: Rahoma@sci.cu.edu.eg

© The Korean Space Science Society. All rights reserved. This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: Oct 9, 2019; Revised: Nov 22, 2019; Accepted: Dec 12, 2019

Abstract

Frozen orbit is an attractive option for orbital design owing to its characteristics (its argument of pericenter and eccentricity are kept constant on an average). Solar sails are attractive solutions for massive and expensive missions. However, the solar radiation pressure effect represents an additional force on the solar sail that may greatly affect its orbital behavior in the long run. Thus, this force must be included as a perturbation force in the dynamical model for more accuracy. This study shows the calculations of initial conditions for a lunar solar sail frozen orbit. The disturbing function of the problem was developed to include the lunar gravitational field that is characterized by uneven mass distribution, third body perturbation, and the effect of solar radiation. An averaging technique was used to reduce the dynamical problem to a long period system. Lagrange planetary equations were utilized to formulate the rate of change of the argument of pericenter and eccentricity. Using the reduced system, frozen orbits for the Moon sail orbiter were constructed. The resulting frozen orbits are shown by two 3Dsurface (semimajor, eccentricity, inclination) figures. To simplify the analysis, we showed inclination–eccentricity contours for different values of semi-major axis, argument of pericenter, and values of sail lightness number.

Keywords: frozen orbit; perturbation; lunar gravity; third-body

1. INTRODUCTION

The mission design to place a satellite in an orbit around a celestial body like the Moon is a complex task; hence, it is time-consuming to choose the perfect orbital initial conditions that would meet mission objectives.

Recently, frozen orbits drew particular attention owing to their space-applications like topographic mapping and planetary environment researches and observations. In such orbits, the argument of pericenter and eccentricity are enforced to be stationary by choosing suitable initial conditions, even if the influences of different perturbations are considered.

The history of research on frozen orbits began several years ago, as detailed in a study by Coffey et al. (1994), where many researchers introduced useful contributions. Smith (1986) used Delaunay variables in the canonical formulations of partial derivatives to predict the long-term evolution of eccentricity and the argument of pericenter. Rosborough & Ocampo (1991) utilized Lagrange planetary equations to derive frozen orbits of the Earth using the gravity mode. Lara et al. (1995) used angular momentum polar component to numerically outline families of frozen orbits of the Earth. The same technique was applied to lunar frozen orbits by Elipe & Lara (2003). Aorpimai & Palmer (2003) used epicycle description to compute the Earth’s frozen orbits up to arbitrary zonal harmonics. Zhigang et al. (2014) determined the initial conditions that minimize the variation of the orbital plane for nearly circular orbits with a semi-major axis 3–10 times the Earth’s radius. Santos et al. (2013) searched for frozen orbits through some missions around planetary satellites, like Jupiter’s moons Europa and Icy, and studied their stability. Masoud et al. (2018, 2019) discussed six control strategies to construct artificial frozen orbits and sun-synchronous orbits around the Earth.

Furthermore, solar sail is a modern propulsion method in which the radiation pressure of sunlight is used to propel a spacecraft. It presents long operating lifetimes in addition to low cost of operation. Although its force is lowthrust compared to that of electric engines, it is a perfect continuous force that could be collected over time to be sufficiently large to be considered as a propulsion force. The solar sailing force exerted on an 800 × 800 m² mirror at one astronomical unit (AU) is about 5 N (Jerome 1992).

Solar sails are applicable to many space missions spread over the entire solar system. Solar sails could reach up to 0.25 AU or even closer to the sun to offer solar observation payloads. They could travel between inner and outer planets. They could also be placed at a stationary point relative to either the Earth or the sun to warn about solar storms (Zubrin et al. 2011). The orbits of Earth’s satellites can be modified by solar sailing. Also, solar sail is suitable for pole sitter and deorbiting (McInnes 2004;Ceriotti & McInnes 2011;Gong et al. 2011;Charlotte et al. 2012).

Wie (2004a, b);Wie & Murphy (2007) studied the solar sail attitude dynamics and how to control it. Cui et al. (2008) derived the dynamical model of multi-body solar sails with control boom. Regarding a low Earth orbit, Polites et al. (2008) presented a numerical simulation on attitude control for a solar sail by using magnetic torquers in addition to small reaction wheels. Liu et al. (2014a, b) presented vibration equations for a solar sail with large flexible established attitude dynamics with control boom, control vanes, and sliding masses.

Carvalho (2016) found the frozen orbits for an artificial Mercuian solar sail considering Mercury’s non-sphericity caused by the zonal terms J₂, and J₃, solar radiation pressure (SRP), and the sun as a third-body.

Carvalho et al. (2012a,b, 2013) in consequent studies found f Jupiter’s moon (Europa) by, considering Jupiter as a third-body using Lie-transform as an average method within Hamiltonian framework. Abd El-Salam (2013) considered the Poynting–Robertson drag effects of solar sail. In heliocentric orbit, he obtained constraints for solar sail optimal dynamics.

This investigation aimed to determine the frozen orbits around the Moon using continuous acceleration by solar propulsion considering the uneven mass distribution and third body attraction, with the solar sailing as a perturbation force.

The gravitational potential of the Moon is quite different from that of the Earth, so that the zonal harmonics up to J₅, first sectorial harmonic C₂₂ and first tesseral harmonic C₃₁ are of the same order (Rahoma 2014;Rahoma & Abd El- Salam 2014).

2. DYNAMICAL MODEL

The present dynamical system comprised three bodies: the central body— the Moon, the solar sail spacecraft, and the third disturbing body (TB), with masses m₀, m, and m’, respectively. The reference frame Oxyz originated at the Moon and the xy-plane coincided with the lunar equator. The x-axis was along the intersection point between the lunar equator and the TB orbital plane. To define the sail orbits around the Moon, we used the classical Keplerian elements: a, orbital semi-major axis; e, eccentricity; i, inclination; ω, argument of pericenter; Ω, longitude of ascending node; and f, the true anomaly. Further, n refers to the mean motion, where n²a³ = Gm₀ with G as the gravitational constant. The TB, m’, is assumed to move in an elliptic inclined orbit around the central body with the orbital elements (a’, e’, i’, ω’, Ω’), and its mean motion n’ is given by n’²a’³ = G(m₀ + m’), with the reference frame definition, Ω’ = 0.

Considering the perturbing accelerations, the solar sail equation of motion around the Moon is described by

r → ¨ = r → ¨ M + r → ¨ 3 b + r → ¨ S R P

(1)

where $r → ¨ M$ is the gravitational acceleration induced by the lunar gravity field, $r → ¨ 3 b$ is the acceleration due to the presence of TB, and $r → ¨ S R P$ is the acceleration due to SRP.

2.1 Lunar Gravity Potential

The lunar gravity can be written as gradient of the potential U_M:

U M = − μ r + μ r ∑ n = 2 ∞ ∑ m = 0 n (R M r) n p n m (sin φ) (C n m cos m λ + S n m sin m λ)

(2)

where the quantities C_nm and S_nm are the lunar harmonics coefficients. μ and R_M are gravitational parameter and lunar equatorial radius, respectively. The quantities (r, λ, φ) are the spherical coordinates of the sail-craft. Further, $p n m (sin φ)$ are the associated Legendre polynomials.

The Moon’s gravity field is rather asymmetric, unlike that of the Earth, where the sectorial coefficient P₂₂ and the tesseral coefficient P₃₁ are larger than the zonal coefficients P₂, P₃, P₄ and P₅; thus, P₂₂ and P₃₁ cannot be ignored.

The following equations were derived based on the assumptions of Giacaglia et al. (1970):

P 2 (sin φ) = 12 [3 sin 2 i sin 2 (f + ω) − 1], P 3 (sin φ) = 12 [5sin 3 i sin 3 (f + ω) − 3 sin i sin (f + ω)], P 4 (sin φ) = 18 [35sin 4 i sin 4 (f + ω) − 30 sin 2 i sin 2 (f + ω) + 3], P 5 (sin φ) = 18 [63sin 5 i sin 5 (f + ω) − 70 sin 3 i sin 3 (f + ω) + 15 sin i sin (f + ω)] P 22 (sin φ) cos 2 λ = 6 [ζ 2 cos 2 f + χ 2 sin 2 f + χ ζ sin 2 f − 3 (1 − sin 2 i sin 2 (f + ω))], P 31 (sin φ) cos λ = 38 (5 sin 2 i − 4) cos Ω cos (f + ω) − 38 (2 + 5 sin 2 i) cos i sin Ω sin (f + ω) − 158 sin 2 i cos Ω cos 3 (f + ω) + 158 sin 2 i cos i sin Ω sin 3 (f + ω)

where we consider only the harmonics coefficients J₂, J₃, J₄, J₅, C₂₂ and C₃₁ (Rahoma 2014;Rahoma & Abd El-Salam 2014).

Also, ζ = cosω cosΩ - cosi sinω sinΩ, χ = -sinω cosΩ - cosi cosω sinΩ.

Using Ψ = a/r, the lunar potential can be rewritten as El-Saftawy (1991)

U M = U M 2 + U M 3 + U M 4 + U M 5 + U M 22 + U M 31

(3)

where

U M 2 = J 2 μ M 4 R M 2 4 L 6 Ψ 3 [(2 − 3 S 2) + 3 S 2 cos 2 (f + ω)], U M 3 = − J 3 μ M 5 R M 3 8 L 8 Ψ 4 [(15 S 3 − 12 S) sin (f + ω) − 5 S 3 sin3 (f + ω)], U M 4 = − J 4 μ M 6 R M 4 64 L 10 Ψ 5 [(105 S 4 − 120 S 2 + 24) ​ + (− 140 S 4 + 120 S 2) cos 2 (f + ω) + 35 S 4 cos 4 (f + ω)], U M 5 = − J 5 μ M 7 R M 5 128 L 12 Ψ 6 [30 S (21 S 4 − 28 S 2 + 8) sin (f + ω) + 35 S 3 (− 9 S 2 + 8) sin 3 (f + ω) + 63 S 5 sin 5 (f + ω)], U M 22 = C 22 μ M 4 R M 2 L 6 Ψ 3 [(3 (ζ 2 − χ 2) - 32 S 2 cos 2 ω) cos 2 f + (3 χ ζ + 32 S 2 sin 2 ω) sin 2 f + (3 (ζ 2 + χ 2) + 32 S 2 − 3)], U M 31 = C 31 3 μ M 5 R M 3 L 6 Ψ 4 [(5 S 2 − 4) cos Ω cos (f + ω) − (2 + 5 S 2) C sin Ω sin (f + ω) − 5 S 2 cos Ω cos 3 (f + ω) + 5 S 2 C sin Ω sin 3 (f + ω)],

where μ_M and R_M are the lunar gravitational parameter and lunar equatorial radius, respectively, $L = μ M a$ , C = cosi and S = sini

2.2 Third-Body Attraction

The acceleration vector $r → ¨ 3 b$ due to TB (Sun) is

r → 3 b ¨ = − μ ′ (r → − r → ′ ‖ r → − r → ′ ‖ 3 + r → ′ ‖ r → ′ ‖ 3)

(4)

where $r →$ ' r′ is the solar position vector with respect to the Moon. The quantity μ′ is the solar gravitational parameter. The acceleration can be read as

U 3 b = μ ′ G (m 0 + m ′) r ′ ∑ n = 2 ∞ (r r ′) n P n (cos (ψ))

which may be written as up to the second order

U 3 b 2 = μ ′ n ′ 2 a 2 2 (a ′ r ′) 3 (r a) 2 [3 cos 2 (ψ) − 1]

(5)

where ψ defines the angle between $r →$ and $r →$ '.

Also ψ can be formulated using the orbital elements as follows:

cos ψ = ∑ i, j = − 1 1 A i j cos (f + ω + i (f ′ + ω ′) + j (Ω ′ − Ω))

with the following non-vanishing coefficients of A_ij (Rahoma, 2016).

A − 1 − 1 = 14 (1 + cos i + cos i ′ + cos i cos i ′) A − 10 = 12 sin i sin i ′ A − 11 = 14 (1 − cos i − cos i ′ + cos i cos i ′) A 1 − 1 = 14 (1 + cos i − cos i ′ − cos i cos i ′) A 10 = − 12 sin i sin i ′ A 11 = 14 (1 − cos i + cos i ′ − cos i cos i ′)

2.3 Effect of the Solar Radiation Pressure

Assuming the solar radiation acceleration ( $r → ¨ S R P$ ) and Sun–sail direction are perpendicular to the direction of motion (Fig. 1), the sail acceleration equals the Sun attraction in magnitude, but in an opposite direction.

Fig. 1. Geometry of orbits, coordinates, vectors, and angles used in the dynamical system.

Download Original Figure

2.4 Solar Sail Lightness Number

Solar sailing is a new technology applied in the propulsion of a spacecraft. The reflective surface utilizes SRP to accelerate the spacecraft, consistent with the strategy adopted for IKAROS (2010), NanoSail-D (2010), and LightSail-1 (2016). The sail acceleration is unlimited and continuous, which helps a long term mission. The sailing efficiency depends mainly on the sail lightness number and orientation.

Lightness number, β, is a dimensionless number, and it can be defined as the ratio of acceleration due to solar radiation to acceleration due to solar gravitational pull (Fu et al. 2016).

β = P m a x (ρ) A G m 0 m ′ / ρ 2 = c o n s t a n t

where ρ is the distance between the sail and the Sun, A is the area of the sail, Pmax(ρ) is the maximum of SRP. It depends on orientation and reflected properties of the sail.

Fig. 2 introduces the values of Pmax which reaches the first five planets in the solar system. On Earth, Pmax is approximately 9 μN/m².

Fig. 2. Maximum SRP for the locations of the first five planets in the Solar System. SRP, solar radiation pressure.

Download Original Figure

In other words, the lightness number reflects the ability of the sail to exploit SRP, irrespective of the sail’s size.

Table 1 shows the lightness number, area-to-mass ratio and the featured acceleration for a sail of mass 10 kg (Farres 2019).

Table 1. Relation between lightness number, area-to-mass ratio, and acceleration

Download Excel Table

The acceleration represents SRP due to the Sun and can be written as the equation given by Tresaco et al. (2016).

r → ¨ S R P = β μ ′ r → − r → ′ ‖ r → − r → ′ ‖ 3

(6)

Using Legendre polynomials to expand equation (6), it returns to the second order as follows:

U S R P = − β μ ′ [r r 2 cos ψ + 1 2 r ′ (r r ′) 2 (3 cos 2 ψ − 1)]

(7)

3. REDUCED SYSTEM

To reduce the dynamical system’s degrees of freedom and eliminate the short-period terms, we present a double- averaging to terminate the short period of the sail and that of TB. The double-averaging for a function F is defined as follows:

〈 F 〉 = 1 4 π 2 ∫ 0 2 π ∫ 0 2 π F d l d l ′

Consistent with the derivations provided by Abd El-Salam et al. (2006), the average procedure on the lunar disturbing potential (wherein the orbital elements, except the mean anomaly, are constants) returns the following:

U M * = 〈 〈 U M 〉 〉 = U M 2 * + U M 3 * + U M 4 * + U M 5 * + U M 22 * + U M 31 *

(8)

with

U M 2 * = A 12 L 3 G 3 (3 S 2 − 2), U M 3 * = − 38 A 13 L 3 G 5 (5 S 3 − 4 S) e sin ω, U M 4 * = − 3 64 A 14 L 3 G 7 (1 + 32 e 2) (35 S 4 − 40 S 2 + 8) − 5 e 2 (7 S 4 − 6 S 2) cos 2 ω], U M 5 * = − 5 128 A 15 L 3 G 9 [(126 S 5 − 168 S 3 + 485) (2 e + 32 e 3) sin ω − (632 S 5 − 28 S 3) e 3 sin 3 ω] U M 22 * = 3 A 122 2 L 3 G 3 [S 2 cos 2 Ω] U M 31 * = 3 A 131 8 L 3 G 5 e [(− 4 + 5 S 2) cos ω cos Ω − C (2 + 5 S 2) sin ω sin Ω]

where

A 12 = J 2 ε μ M 4 R M 2, A 13 = J 3 ε μ M 5 R M 3, A 14 = J 4 ε μ M 6 R M 4, A 15 = J 5 ε μ M 7 R M 5, A 122 = C 22 ε μ M 4 R M 2, A 131 = C 31 ε μ M 5 R M 3 .

Averaging over the TB period yields the following (Rahoma 2016):

U 3 b 2 * = 〈 〈 U 3 b 2 〉 〉 = μ ′ n ′ 2 a 2 ∑ j, k = − 2 2 {α j k (e, e ′) ρ j k (i, i ′) cos (j ω + k (Ω ′ − Ω))}

(9)

where the non-vanishing coefficients are as follows:

α 00 = (2 + 3 e 2) 4 (1 − e ′ 2) 3 / 2, α 01 = α 02 = 3 (2 + 3 e 2) 4 (1 − e ′ 2) 3 / 2, α 20 = α 2 − 2 = α 2 − 1 = α 21 = α 22 = 15 e 2 (2 + 3 e 2) 4 (1 − e ′ 2) 3 / 2 ρ 00 = − 1 + 32 (A − 1 − 1 2 + A − 10 2 + A − 11 2 + A 1 − 1 2 + A 102 + A 112), ρ 01 = A − 1 − 1 A − 10 + A − 10 A − 11 + A 1 − 1 A 10 + A 10 A 11, ρ 02 = A − 1 − 1 A − 11 + A 1 − 1 A 11, ρ 20 = A − 11 A 1 − 1 + A − 10 A 10 + A − 1 − 1 A 11 ρ 2 − 2 = A − 1 − 1 A 1 − 1, ρ 2 − 1 = A − 10 A 1 − 1 + A − 1 − 1 A 10, ρ 21 = A − 11 A 10 + A − 10 A 11, ρ 22 = A − 11 A 11

Averaging Eq. (7) upon the spacecraft’s eccentric anomaly yields

U S R P * = 〈 〈 U S R P 〉 〉 = − β μ ′ n ′ 2 a 2 ∑ j, k = − 2 2 {α j k (e, e ′) ρ j k (i, i ′) cos (j ω + k (Ω ′ − Ω))} = − β 〈 〈 U 3 b 2 〉 〉

(10)

4. COMPUTATION OF FROZEN ORBITS

To formulate the time derivative of pericenter argument and eccentricity, Lagrange planetary equations introduced by Brouwer & Clemence (1961) yields:

d ω d t = (− 1 − e 2 n a 2 e ∂ ∂ e + cot i n a 2 1 − e 2 ∂ ∂ i) (U M * + U 3 b 2 * + U S R P *)

(11)

d e d t = (1 − e 2 n a 2 e ∂ ∂ ω − 1 − e 2 n a 2 e ∂ ∂ i) (U M * + U 3 b 2 * + U S R P *)

(12)

Frozen orbit can be guaranteed by elimination of the rate of change of pericenter argument and eccentricity of the orbit, which are referred to as the equilibrium points in the dynamical system:

d e d t = 0

(13.1)

d ω d t = 0

(13.2)

This procedure is essentially dependent on the effect of the zonal harmonics terms.

Eq. (11) depends mainly on the orbital semi-major axis, eccentricity, and inclination. Therefore, this equation can be represented by a 3D surface (a, e, i) . Each point on that surface represents initial conditions of the frozen orbits.

5. RESULTS AND ANALYSIS

In general, assume that the reference frame had its x-axis coincide with the TB nodal line at the initial time point, i.e., Ω′ = 0, ω = π/2 or 3π/2, and Ω - Ω′= 2kπ, where k is an integer.

Consistent with the findings by Meyer at al. (1994), the required data and orbital parameters of the Moon’s referenced orbit applied in simulation for frozen orbits are mentioned in Table 2.

Table 2. Dynamical model parameters

Download Excel Table

To derive the lunar frozen orbits, we used equation (13) in equation (11) to yield a solution dependent on (a, e, i), which provided a 3D surface as shown in Figs. 3 and 4.

Fig. 3. Frozen orbit surface (a, i, e) for

ω = π 2

; β=0.2.

Download Original Figure

Fig. 4. Frozen orbit surface (a, i, e) for

ω = 3 π 2

; β=0.2.

Download Original Figure

To perform detailed analyses, we presented cross sections at different values for semi-major axis and graph for the entire potential range of the eccentricity and inclination values, i.e., the initial conditions (e, i) of the frozen orbits were depicted for different values of α and β.

Figs. 5 and 6 represent the section (e, i) at β = 0.2 for chosen semi-major (a) 3,500, 5,500, and 7,500 km; ω = 3π/2 and π/2, respectively.

Fig. 5. Section (i, e) of frozen orbits where β = 0.2, ω = π/2 and different value of a.

Download Original Figure

Fig. 6. Section (i, e) of frozen orbits where β = 0.2, ω = 3π/2 and different values of a.

Download Original Figure

Figs. 7 and 8 show the sections when β = 0.2, 0.8 for α = 2,580 km, and ω = π/2, 3π/2, respectively.

Fig. 7. Section (i, e) of frozen orbits where α = 2,580 km, β = 0.8, and

ω = π 2, 3 π 2

Download Original Figure

Fig. 8. Section (i, e) of frozen orbits where α = 2,580 km, β = 0.2, and

ω = π 2, 3 π 2

Download Original Figure

From the graphs, we can observe the following:

1- The eccentricity value of frozen orbits increases when the inclination (i) changes from 0 to π and approaches the families of critical inclinations (inclinations at which apocenter drift is zero) (Rahoma & Abd El- Salam 2014).
2- The effects due to SRP and TB increase with the increase in the orbital semi-major axis; refer Figs. 7 and 8.
3- SRP perturbation is primarily based on the sail’s areato- mass ratio, which is expressed in the dynamic model within β in addition to the orbit semi-major axis.
4- At polar orbits, the eccentricity value of the frozen orbits increases as β increases.
5- Although there are different initial conditions, we observed from the 2D Figs. 5 and 6 that the frozen orbits are nearly the same at ω = π/2 and ω = 3π/2, except for the few dashed lines added to the left-hand side of the figures when ω = π/2.
6- We noticed a significant change in Figs. 7 and 8 depending on the different values for lightness number of the sail, β.

6. CONCLUSIONS

A good understanding of non-integrable dynamics is essential for deriving the locations of the frozen orbits and designing future Moon missions; furthermore, for Moon missions, if we select one parameter of a frozen orbit, we can obtain the remaining two parameters.

The dynamics of frozen orbits of a solar sail planet’s orbiter type has been introduced (applied here to the Moon). The model includes first order of the gravitational sectorial and tesseral harmonics, which in case of bodies like the Moon are sufficiently big and cannot be neglected. These harmonics are considered to be of the same order as that of the zonal harmonics up to J₅. Additionally, the model involves TB attraction and SRP. By performing double-average techniques, the short-period terms are eliminated. The study included the third-body’s inclination and eccentricity, which gives an additional complexity to the dynamical equations. In such scenarios, the initial conditions of the frozen orbits are introduced.

A 3D surface was derived where each point represented the initial conditions for lunar uneven mass associated with the chosen frozen orbits, considering different perturbations in formulation equations (11 and 12) of the dynamic model. The sail’s area-to-mass ratio (A/m) has a strong effect on sail-craft dynamics, given the lightness number value, i.e., area-to-mass ratio for the solar sail must be well-planned for future missions. In the future, we aim to extend this study to include effects of the Moon’s gravitational higher harmonics in the dynamics. Also, the attraction of both the Sun and Earth can be enforced in calculations.

ACKNOWLEDGEMENTS

The authors would like to thank the entire JASS editorial board and referees for their insightful comments, which improved the manuscript.

References

Abd El-Salam FA, El-Tohamy IA, Ahmed MK, Rahoma WA, Rassem MA, Invariant relative orbits for satellite constellations: a second order theory, Appl. Math. Comput. 181, 6-20 (2006).

Abd El-Salam FA, Some new locally optimal control laws for sailcraft dynamics in heliocentric orbits, J. Appl. Math. 2013, 353056 (2013).

Aorpimai M, Palmer PL, Analysis of frozen conditions and optimal frozen orbit insertion, J. Guid. Control Dynam. 26, 786-793 (2003).

Brouwer D, Clemence GM, Methods of Celestial Mechanics (Academic Press, New York, NY, 1961).

Carvalho JPS, Elipe A, de Moraes RV, Prado AFBA, Low-altitude, near-polar and near-circular orbits around Europa, Adv. Space Res. 49, 994-1006, (2012a).

Carvalho JPS , Mourao DC, Elipe A, Vilhena de Moraes R, Prado AFBA, Frozen orbits around the Europa, Int. J. Bifurcat. Chaos. 22, 1250240 (2012b).

Carvalho JPS, de Moraes RV, Prado AFBA, Dynamics of artificial satellites around Europa, Math. Probl. Eng. 2013, 182079 (2013).

Carvalho JPS, Orbital evolution of a solar sail around a planet, Proc. Ser. Braz. Soc. Comp. Appl. Math. 4, 010017 (2016). https://doi.org/10.5540/03.2016.004.01.0017

Ceriotti M, McInnes CR, Generation of optimal trajectories for earth hybrid pole sitters, J. Guid. Control Dynam. 34, 847-859 (2011). https://doi.org/10.2514/1.50935

10.

Charlotte L, Camilla C, McInnes C, Solar radiation pressure augmented deorbiting from high altitude sun-synchronous orbits, Proceedings of the 4S Symposium 2012, Small Satellites Systems and Services, Portoroz, Slovenia, 4-8 Jun 2012.

11.

Coffey SL, Deprit A, Deprit E, Frozen orbits for satellites close to an Earth-like planet, Celest. Mech. Dyn. Astr. 59, 37-72 (1994).

12.

Cui HT, Luo JH, Feng JH, Cui PY, Attitude control of solar sail spacecraft with control boom, J. Astronaut. 29, 560-566 (2008).

13.

Elipe A, Lara M, Frozen orbits about the Moon, J. Guid Control Dynam. 26, 238-243 (2003).

14.

El-Saftawy MI, Motion of a lunar artificial satellite, Master Thesis, Cairo University (1991).

15.

Farres A, Heiligers J, Miguel N, Road Map to L4/L5 with a solar sail, Aerosp. Sci. Technol. 95, 105458 (2019).

16.

Fu B, Sperber E, Eke F, Solar sail technology: a state of the art review, Prog. Aerosp. Sci. 86, 1-19 (2016).

17.

Giacaglia GEO, Murphy JP, Felsentreger TL, A semi-analytic theory for the motion of a lunar satellite, Celest. Mech. 3, 3-66 (1970).

18.

Gong S, Li J, Baoyin H, Solar sail transfer trajectory from L1 point to sub-L1 point, Aerosp. Sci. Technol. 15, 544-554 (2011).

19.

Jerome W, Space Sailing (Gordon and Breach, Yverdon, Switzerland, 1992).

20.

Lara M, Deprit A, Elipe A, Numerical continuation of families of frozen orbits in the zonal problem of artificial satellite theory, Celest. Mech. Dynam. Astr. 62, 167-181 (1995).

21.

Liu J, Cui N, Shen F, Rong S, Dynamics of highly-flexible solar sail subjected to various forces, Acta Astronaut. 103, 55-72 (2014a).

22.

Liu J, Rong S, Shen F, Cui N, Dynamics and control of a flexible solar sail, Math. Probl. Eng. 2014, 868419 (2014b).

23.

Masoud A, Rahoma WA, Khattab EH, El-Salam FA, Construction of frozen orbits using fontinuous thrust control theories considering earth oblateness and solar radiation pressure perturbations, J. Astronaut. Sci. 65, 448-469 (2018).

24.

Masoud A, Rahoma WA, Khattab EH, El-Salam FA, Design of artificial sun-synchronous orbits with main zonal harmonics and solar radiation pressure using continuous low-thrust control strategies, Open Astron. J. 28, 124-130 (2019).

25.

McInnes CR, Solar Sailing: Technology, Dynamics and Mission Applications (Springer Praxis, London, UK, 2004).

26.

Meyer KW, Buglia JJ, Desai PN, Lifetimes of lunar satellite orbits, NASA Technical Paper 3394 (1994).

27.

Polites M, Kalmanson J, Mangus D, Solar sail attitude control using small reaction wheels and magnetic torquers, Proc. Inst. Mech. Eng. Part G J. Aerosp. Eng. 222, 53-62 (2008).

28.

Rahoma WA, Orbital elements evolution due to a perturbing body in an inclined elliptical orbit, J. Astron. Space Sci. 31 199-204 (2014).

29.

Rahoma WA, Investigating exoplanet orbital evolution around binary star systems with mass loss, J. Astron. Space Sci. 33, 257-264 (2016).

30.

Rahoma WA, Abd El-Salam FA, The effects of moon’s uneven mass distribution on the critical inclinations of a lunar orbiter, J. Astron. Space Sci. 31, 285-294 (2014).

31.

Rosborough GW, Ocampo CA, Influence of higher degree zonals on the frozen orbit geometry, Proceedings of the AAS/AIAA Astrodynamics Conference, Durango, CO, 19-22 Aug 1991.

32.

Santos JC, de Moraes RV, Carvalho JS, Stability of frozen orbits around Europa, in American Astronomical Society, DDA meeting #44, id.204.12, May 2013.

33.

Smith, JC, Analysis and application of frozen orbits for the TOPEX mission, in AAS/AIAA Astrodynamics Conference, Williamsburg, VA, 18-20 Aug 1986.

34.

Tresaco E, Elipe A, Carvalho, JPS, Frozen orbits for a solar sail around mercury. J. Guid. Control Dynam. 39, 1659-1666 (2016).

35.

Wie B, Murphy D, Solar-sail attitude control design for a sail flight validation mission, J. Spacecr. Rocket. 44, 809-821 (2007). https://doi.org/10.2514/1.22996

36.

Wie B, Solar sail attitude control and dynamics, part 1, J. Guid. Control Dynam. 27, 526-535 (2004a).

37.

Wie B, Solar sail attitude control and dynamics, part two, J. Guid Control Dynam. 27, 536-544 (2004b).

38.

Zhigang W, Jiang F, Li J, Artificial Martian frozen orbits and Sun-synchronous orbits using continuous low-thrust control, Astrophys. Space Sci. 352, 503-514, (2014).

39.

Zubrin R, Wagner R, Clarke AC, The Case for Mars: The Plan to Settle the Red Planet and Why We Must. (Free Press, New York, NY, 2011).