On the design of a space telescope orbit around the Sun–Venus L2 point

Available online at www.sciencedirect.com

ScienceDirect Advances in Space Research xxx (2020) xxx–xxx www.elsevier.com/locate/asr

On the design of a space telescope orbit around the Sun–Venus L2 point Maksim Shirobokov ⇑, Sergey Trofimov, Mikhail Ovchinnikov Keldysh Institute of Applied Mathematics, 4 Miusskaya Pl., Moscow 125047, Russia Received 7 July 2019; received in revised form 25 October 2019; accepted 19 December 2019

Abstract In this study, the choice of a nominal space telescope orbit in the vicinity of the Sun–Venus L2 libration point is discussed from the viewpoint of illumination conditions and station-keeping costs. Such a location for the space telescope is especially appealing for the observation of potentially hazardous asteroids approaching the Earth from the daytime side of the sky. In contrast to the case of a telescope placed near the Sun–Earth L1 point, a significantly longer warning time can be achieved. Moreover, a Lissajous-type quasi-periodic orbit can be selected so that a predefined percentage of the orbit is shadowed. When a nominal orbit is allowed to be permanently sunlit, a large halo orbit is preferred due to the lower station-keeping cost. For two sets of unstable halo and Lissajous orbits, the stationkeeping cost is evaluated by conducting Monte Carlo simulations in the ephemeris model of motion. Two modifications of the target point station-keeping technique are examined: the X-axis control and the 3-axis control. The simulation scenario includes orbit insertion and navigation errors, impulse execution errors, and constraints on the minimum imparted Dv. Dependence of the station-keeping cost on orbit insertion and navigation errors is analyzed. Some of the large halo orbits appear to be linearly stable. The corresponding central manifold location in the phase space is determined, which makes a simple targeting strategy possible. Ó 2019 COSPAR. Published by Elsevier Ltd. All rights reserved.

Keywords: Station-keeping; Sun–Venus system; Space observatory; Libration point; Halo orbit; Lissajous orbit

1. Introduction After 40 years of practical use, libration point orbits (LPOs) have become reliable and common places for space observatories. These orbits possess useful geometric and dynamic properties allowing observations of the Sun, space weather, and distant stars in different ranges of wavelength. Table 1 contains past and present missions involved LPOs. Since all those LPOs are unstable, station-keeping is required to maintain the spacecraft trajectory in the vicinity of a nominal orbit. Annual station-keeping costs are presented in the last column of Table 1. Bibliographical references for almost all the presented data can be found in (Shirobokov et al., 2017). Station-keeping information

⇑ Corresponding author.

about the Deep Space Climate Observatory (DSCOVR) mission and Queqiao spacecraft are taken from Roberts et al. (2016) and Gao et al. (2019), respectively. There are a lot of publications devoted to LPO stationkeeping; for their classification see the survey (Shirobokov et al., 2017). In (Shirobokov et al., 2017), station-keeping methods are classified into the two groups: the methods that directly exploit the dynamical effects of the threebody problem (Floquet modes, Lyapunov exponents, etc.) and the methods that are based on optimal control theory (nonlinear programming techniques, targeting approaches, linear quadratic regulator, etc.). In practice, targeting approaches are frequently used since they are reliable and simple in implementation and incorporation into the real mission design. Station-keeping costs depend on many factors, such as the type and size of the nominal orbit and the three-body system, the navigational and thrust exe-

E-mail address: [email protected] (M. Shirobokov). https://doi.org/10.1016/j.asr.2019.12.022 0273-1177/Ó 2019 COSPAR. Published by Elsevier Ltd. All rights reserved.

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

2

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

Table 1 Missions and station-keeping characteristics. Mission

Agency

Launch date

Libration point(s)

Type of orbit

AY ; AZ , 1000 km

Dv m/s/year

ISEE-3 Wind SOHO ACE WMAP Genesis

NASA/ESA NASA ESA/NASA NASA NASA NASA

12 01 02 25 30 08

ARTEMIS

NASA

17 Feb 2007

Herschel Planck Chang’e 2 Gaia Chang’e 5-T1 DSCOVR LISA Pathfinder Queqiao Spektr-RG

ESA/NASA ESA CNSA ESA CNSA NASA/NOAA ESA CNSA Roscosmos

14 14 01 19 23 11 03 20 13

SE L1 SE L1 SE L1 SE L1 SE L2 SE L1 EM L2 EM L1 EM L1 SE L2 SE L2 SE L2 SE L2 EM L2 SE L1 SE L1 EM L2 SE L2

Halo Quasi-halo Halo Lissajous Lissajous Quasi-halo Quasi-halo Quasi-halo Quasi-halo Halo Lissajous Lissajous Lissajous Lissajous Lissajous Quasi-halo Halo Lissajous

667, 120 640, 170 667, 120 264, 157 264, 264 800, 450 63.52, 35.20 (P1) 58.82, 2.39 (P1) 67.71, 4.68 (P2) 700, 400 300, 300 918, 400 350, 100 40, 35 264, 157 800, 600 35, 16 800, 600

8.5 1.0 2.4 1.0 1.2 9.0 7.39 (P1) 5.28 (P1) 5.09 (P2) 1.0 1.0 (N/A) 2.0 (N/A) 2 1.8 18.0 (N/A)

Aug 1978 Nov 1994 Dec 1995 Aug 1997 Jun 2001 Aug 2001

May 2009 May 2009 Oct 2010 Dec 2013 Oct 2014 Feb 2015 Dec 2015 May 2018 Jul 2019

cution errors, the minimal possible thrust value that the thrusters can execute, the presence of constraints on the thrust vector direction. The majority of the investigations about station-keeping around libration point orbits concern the Sun–Earth and the Earth–Moon systems. It is expected since there are many real missions in these three-body systems and because of a lower access cost when compared to other three-body systems. However, ideas of exploiting LPOs in the Sun–Venus system appeared in recent years. Libration point orbits around the Sun–Venus L2 point could serve as a good location for telescopes that observe hazardous asteroids whose aphelion is closer than the Earth’s perihelion (Dunham and Genova, 2010; Bodrova and Raikunov, 2017). These objects are called the Interior Earth Objects (IEOs) and are difficult in monitoring from the Earth, the near-Earth orbits, or the libration point orbits around the Sun–Earth L1 or L2 due to the small solar elongation angle (Jedicke et al., 2003; Greenstreet et al., 2012; Dunham and Genova, 2010). Moreover, Sun-Venus L2 LPOs provide a promising solution for determination of physical and chemical properties and general geological composition of the IEOs from the estimation of amount of reflected radiation flux (Bodrova and Raikunov, 2017). An observatory placed around the Sun–Venus L2 can extend our knowledge about dangerous asteroids and can be used in pair with an observatory placed around the Sun–Earth L1 =L2 points. Other studies propose using Venus-like orbits to cover the space regions of interest (Gold, 2001; Jedicke et al., 2003; Council, 2010; Akhmetshin, 2013; Lu et al., 2013; Dunham et al., 2013). The advantage of LPOs around the Sun–Venus L2 point over heliocentric Venus-like orbits is an additional option to choose a level of illumination (amount of light received from the Sun taking into account the eclipses by Venus) the spacecraft is exposed to. Indeed, the vicinity of the L2 point

contains the penumbra zones of different illumination levels, which can be used to help cooling the spacecraft and prevent overheating. In comparison to the orbits around Venus, there is also no necessity to rotate the orbit for keeping the favorable shadow conditions. In this paper, we estimate station-keeping costs of potentially useful libration point orbits around the Sun– Venus L2 point, taking into account the possible application in asteroid detection missions. On the one hand, requirements on low temperature regimes at which infrared optics work (Bodrova and Raikunov, 2017) can be expressed as a requirement on the minimum length of time intervals when a spacecraft is in the shadow of Venus. On the other hand, sufficient time is needed to recharge the onboard batteries. These conflicting restrictions limit the choice of a feasible operational LPO. We analyze different classes of libration point orbits (halo and Lissajous) and describe the most promising options. Moreover, a range of linearly stable halo orbits is considered, and navigational conditions that lead to the bounded motion are analyzed. The structure of the paper is the following. First, the dynamics in the vicinity of the libration point are described, and the well-known notions of the Lissajous and halo orbits are recalled within the framework of the circular restricted three-body problem. Second, the full ephemeris model augmented with the cannonball solar radiation pressure model is introduced, and the equations of motion are presented. It is also discussed which perturbations are taken into account. Then, the formulas for calculating the illumination level are given, and the orbits considered in this paper are listed. Finally, the stationkeeping simulation process is described in detail, including the optimization procedures. The two target point schemes of station-keeping are considered: the X-axis control and the 3-axis control. The results of the Monte Carlo simula-

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

3

Fig. 1. Rotating frame in the circular restricted three-body problem.

tions are provided for both Lissajous and halo orbits; sensitivity to the navigation errors is discussed. Bounded motion is then studied for the linearly stable part of the halo family. Finally, the conclusion section contains the summary of the research.

and U X ; U Y and U Z are the partial derivatives of U with respect to the position variables. The distances from the spacecraft to m1 and m2 are given by the equalities 2

D2S ¼ ðX þ lÞ þ Y 2 þ Z 2 2

D2V ¼ ðX  1 þ lÞ þ Y 2 þ Z 2 2. Dynamics in the vicinity of a collinear libration point 2.1. Equations of motion In the circular restricted three-body problem (CR3BP), two masses m1 and m2 6 m1 move in circular orbits about their barycenter C, and a spacecraft of negligible mass moves under the gravitational attraction of m1 and m2 . In the CR3BP model, the system of two masses m1 and m2 defines the force field uniquely. If m1 is the Sun and m2 is the Venus, then we obtain the Sun–Venus system. The equations of motion are usually written in the standard rotating coordinate frame (Fig. 1) with the origin at C; the X-axis connects the masses m1 and m2 towards m2 , the Z-axis is directed along the angular velocity of the orbital motion of m2 around m1 , and the Y-axis completes the right-handed system. It is also convenient to use a dimensionless system of units in which 1) masses are normalized so that m1 ¼ 1  l and m2 ¼ l where l ¼ m2 =ðm1 þ m2 Þ is the mass parameter, 2) the angular velocity of the rotating frame is normalized to one, and 3) the distance between m1 and m2 is normalized to one. Then, m1 and m2 are at fixed positions along the X-axis at ½l; 0; 0 and ½1  l; 0; 0, respectively. For the Sun-Venus system, the mass parameter l, units of distance LU, velocity VU, and time TU are listed in Table 2. The equations of motion can be expressed as follows (Szebehely, 1967): X€  2Y_ ¼ U X Y€ þ 2X_ ¼ U Y Z€ ¼ U Z

ð1Þ

The system defined by Eq. (1) has five equilibrium points; three of them lie on the X-axis and are called collinear libration points. Denoted by L1 ; L2 , and L3 , these points are proved to be unstable. Their dimensionless Xcoordinates for the Sun–Venus system are given in Table 2. For telescope mission purposes, the spacecraft should be located behind Venus (as seen from the Sun), i.e. in the vicinity of the L2 point. So, in the remainder of this paper, we pay our attention only to the L2 point. Its distance to Venus is given in Table 2. 2.2. Design of Lissajous and halo orbits The type of the L2 point is saddle  center  center, and the vicinity of L2 is filled with unstable periodic and quasiperiodic orbits Go´mez et al. (2001). For the purpose of this work, we consider only the Lissajous orbits and the halo orbits most commonly encountered in applications. Lissajous orbits are quasi-periodic trajectories, (see Fig. 2). Hereinafter, AX ; AY , and AZ stand for the X-, Y-, and Z-amplitudes of an orbit, respectively. Lissajous orbits are uniquely defined by two parameters: the Y- and Zamplitudes, for example. When looking from Venus or the Sun, these orbits densely fill a rectangular with the center at the X-axis. Construction of Lissajous orbits can be done in the following way: 1. Calculate the Lindstedt–Poincare´ coefficients1 of high degree according to the algorithm described in (Go´mez et al., 2001).

where X2 þ Y2 1  l l þ U ðX ; Y ; ZÞ ¼ þ DS DV 2

1

All the coefficients of the Lindstedt–Poincare´ series expansion up to the 15th order can be found on https://keldysh.ru/microsatellites/eng/software.html.

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

4

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

Table 2 Physical parameters and units in the Sun–Venus system adopted in this paper. Parameter

Name

Value

l LU VU TU X L1 X L2 X L3 cV

Mass parameter Distance unit Velocity unit Time unit X-coordinate of L1 X-coordinate of L2 X-coordinate of L3 Distance from Venus to L2

cS

Distance from Sun to L2

RS RV A=m

Radius of the Sun Radius of Venus Telescope area-to-mass ratio

2:447832293839728  106 1:082094745373792  108 km 3:502056725326771  101 km/s 3:576254104494544  101 days 9:906822995246675  101 1:009371016496719 1:000001019930122 9:373464329012759  103 = ¼ 1:014297649637338  106 km 1:009373464329013= ¼ 1:092237721870165  108 km 695; 700 km 6; 051:8 km 0.01 m2/kg

Fig. 2. Lissajous orbit with AX ¼ 3; 400 km, AY ¼ 10; 842 km, AZ ¼ 1; 900 km; Venus is at the origin of the frame. The arrow indicates the direction of the spacecraft’s motion.

2. Choose two parameters that uniquely identify the orbit: the Y- and Z- amplitudes, for example; calculate a few phase states (patch points) in the approximate orbit. 3. Run a multiple shooting scheme to eliminate all the residues at the patch points. As a result, a Lissajous orbit with the specified amplitudes is obtained provided the order of the Lindstedt–Poincare´ approximation is high enough (15th order is enough for our purposes). This orbit can be saved as a series of the refined patched points, and the intermediate phase states can be recovered by integrating the equations of motion. Halo orbits are three-dimensional periodic orbits symmetric with respect to the XZ-plane. They bifurcate from planar periodic orbits (also known as planar Lyapunov

orbits) at the specific energy level depending on l. In the Sun–Venus system, bifurcation occurs when AY ¼ 455; 672 km. Fig. 3 shows the example of a halo orbit with the Y-amplitude 489; 290 km and the Z-amplitude 200; 000 km. Since halo orbits are non-symmetrical with respect to the XY-plane, the Z-amplitude of a halo orbit refers to the maximum jZj along the orbit. Halo orbits form a one-parameter family, beginning on the XY-plane and ending with Venus-impacting orbits. For sufficiently small halo orbits (AZ < 1; 247; 643 km for the Sun-Venus L2 case), one can choose AZ as the parameter; in general, the maximum X-coordinate can be chosen. Fig. 4 shows the dependence of AZ when moving along the family. Halo orbits can be constructed by the following scheme:

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

5

Fig. 3. Halo orbit with AY ¼ 489; 290 km, AZ ¼ 200; 000 km; Venus is at the origin of the frame. Arrows indicate the direction of the spacecraft’s motion.

Fig. 4. Z-amplitude of halo orbits as a function of the maximum X-coordinate (family parameter).

1. Calculate the Lindstedt-Poincare´ coefficients2 of high degree according to the algorithm described in (Go´mez et al., 2001). 2. Choose a parameter that uniquely identify the orbit and construct the Lindstedt–Poincare´ approximation of the orbit. 3. Apply the differential correction method (Howell, 1984) and refine the approximate Lindstedt–Poincare´ approximation.

As a result, a halo orbit with the specified parameter is obtained provided the order of the approximation is high enough (in this work, the 9th order is used). It is also convenient to prepare a table of preliminarily constructed orbits in order to design new orbits by applying a linear interpolation (estimating the phase state in the orbit) and then refining it by the differential correction technique.

2.3. Stable halo orbits 2

All the coefficients of the Lindstedt–Poincare´ series expansion up to the 9th order can be found on https://keldysh.ru/microsatellites/eng/software. html.

Stability of periodic orbits is determined by the eigenvalues of the corresponding monodromy matrix. According to

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

6

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

Fig. 5. Maximum absolute value among all eigenvalues when AZ increases; the linearly stable region corresponds to AZ in range from 1; 239; 100 km to 1; 246; 900 km; the circle and the blue dots represent respectively the unit circle in the complex plane and the eigenvalues: the left diagram is for AZ ¼ 1; 237; 000 km, the bottom diagram is for AZ ¼ 1; 242; 000 km, and the right diagram is for AZ ¼ 1; 247; 000 km.

the stable manifold theorem for periodic orbits (Perko, 2001), if all the eigenvalues lay strictly inside the unit circle on the complex plane (jki j < 1), then the orbit is stable; if there is an eigenvalue outside of the unit circle (jki j > 1), then the orbit is unstable. In the case when all the eigenvalues lay on the boundary of the unit circle (jki j ¼ 1), the orbit is called linearly stable. To calculate the monodromy matrix associated with the point X0 ¼ ½X 0 ; Y 0 ; Z 0 ; X_ 0 ; Y_ 0 ; X_ 0  of the periodic orbit XðtÞ, one needs to solve the variational equations

the solution and the existence of the first integral in the CR3BP model (Jacobi integral), there are at least two unity eigenvalues (Meyer et al., 2008). For sufficiently small halo orbits, the monodromy matrix has also two real non-unity eigenvalues and two complex conjugate eigenvalues lying on the unit circle, denoting the center manifolds, so we have

_ t0 Þ ¼ AðtÞUðt; t0 Þ Uðt; _ XðtÞ ¼ FðXÞ

where k5 and k6 are complex conjugates. When moving along the family, the eigenvalues shift on the complex plane: k1 ! 1; k2 ! 1; k5 and k6 move along the unit circle, and k3 ¼ k4 ¼ 1 stay at 1. At some moment, all the eigenvalues appear to be on the unit circle, and the orbit becomes linearly stable. The region of stability ends when one of the eigenvalues moves outside the unit circle. The region of linearly stable orbits is shown in Fig. 4: when the maximum X-coordinate decreases, AZ increases. Shortly after the linearly stable region, AZ decreases as the maximum X-coordinate decreases. Fig. 5 shows the evolution of max jki j when AZ increases. Stable halo orbits

ð2Þ ð3Þ

with the initial conditions Uðt0 ; t0 Þ ¼ I66 ; Xðt0 Þ ¼ X0 . Here _ AðtÞ is the Jacobian of the phase state X ¼ ½X ; Y ; Z; X_ ; Y_ ; Z; _ _ _ _ _ FðXÞ ¼ ½X ; Y ; Z; 2Y þ U X ; 2X þ U Y ; U Z : 3 2 0 0 0 1 0 0 6 0 0 0 0 1 07 7 6 7 6 6 0 0 0 0 0 17 7 AðtÞ ¼ 6 7 6U 6 XX U XY U XZ 0 2 0 7 7 6 4 U XY U YY U YZ 2 0 0 5 U XZ

U YZ

U ZZ

0

0 0

X¼XðtÞ

and Uðt; t0 Þ is the state transition matrix. Integrating Eqs. (2) and (3) over the orbital period P, one can compute the monodromy matrix Uðt0 þ P; t0 Þ. It is well known that the eigenvalues of the monodromy matrix do not depend on the reference point X0 and serve as the characteristics of the orbit. Because of the symplectic nature of the state transition matrix, the eigenvalues of the monodromy matrix occur in reciprocal pairs. Due to the periodicity of

k1 > 1; k2 ¼ k1 1 < 1; k3 ¼ k4 ¼ 1 k5 ¼ k6 ; jk5 j ¼ jk6 j ¼ 1

i

correspond to AZ within the range from 1; 239; 100 km to 1; 246; 900 km. The black circle and the blue dots represent respectively the unit circle and the eigenvalues when the Z amplitude is 1,237,000 km (the left diagram), 1,242,000 km (the bottom diagram), and 1,247,000 km (the right diagram). 3. Full model of motion All the results of this work are obtained in the high-fidelity model based on the Jet Propulsion Laboratory DE423 plan-

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

7

Table 3 Sources of perturbations and their characteristics. Source of perturbation

Gravitational parameter (non-dimensional) 1

9:999975521677061  10 1:660132731535470  107 2:447832293839728  106 3:003482264085227  106 3:227148138320988  107 9:547895611422649  104 2:858849702892477  104 4:366233355518391  105 5:151371103455786  105 N/A

Sun Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune SRP

etary ephemerides3 and augmented with the cannonball solar radiation pressure (SRP) model. The SRP force is calculated by taking into account that Venus fully or partially eclipses the Sun, as visible from the spacecraft. We further discuss which orbits can be chosen in real missions with respect to the mean level of illumination along the orbit. Since the station-keeping techniques considered in this research include optimization, the Jacobian matrix of the objective function is needed. Therefore, not only the equations of motion but also the variational equations are presented. Upon designed in the CR3BP model, Lissajous and halo orbits can be adapted to the high-fidelity model with the help of the multiple shooting technique by dividing the trajectory into several arcs and minimizing the residue phase state vectors at the patch points (Ascher et al., 1995; Dei and Topputo, 2017). The optimization can be performed by the Newton method or its modifications, the Jacobi matrix can be obtained by integrating the variational equations, also listed below. Since this technique is standard, it is not described here in detail. 3.1. Equations of motion Consider an inertial coordinate system originating at Venus. Let ri ðtÞ and vi ðtÞ; i ¼ 1; . . . ; N , be correspondingly the known position and velocity of the ith celestial body with respect to Venus at time t and N be the number of celestial bodies (except Venus) that are taken into account. Let r and v be the position and velocity of the spacecraft. Then the equations of motion can be written as follows: 8 > < r_ ¼ v N   X ð4Þ lV r rri ri li jrr þ aSRP > 3 þ 3 : v_ ¼  jrj3  j jr j i¼1

i

i

where lV is the gravitational parameter of Venus, li is the gravitational parameter of the ith celestial body, and aSRP is the acceleration provided by the solar radiation pressure force acting on the spacecraft:

3 The ephemerides can be found at https://ssd.jpl.nasa.gov/?planet_ eph_export and within the planetEphemeris routine in MATLAB.

Normalized acceleration 4:65  101 1:84  108 1:00  100 1:81  107 6:77  109 8:68  107 3:37  108 6:72  1010 2:28  1010 2:60  104

aSRP ¼ qðr; rS Þ  C 

         

8:63  101 2:31  106 1:00  100 4:23  105 3:45  107 3:17  106 9:14  108 1:81  109 6:53  1010 3:23  104

A r  rS  m jr  rS j3

where C ¼ 4:56  106 N/m2 is solar radiation pressure at 1 AU, A=m is the area-to-mass ratio (0:01 m2/kg in this work), q 2 ½0; 1 is the illumination coefficient at the point r, and rS is the position of the Sun. The illumination coefficient q ¼ 0 in the full shadow of Venus, q 2 ð0; 1Þ when Venus only partially eclipses the Sun, and q ¼ 1 in the full light. The formulas for calculating q are presented in the next subsection. In this paper, the epoch t ¼ 0 corresponds to Jan 1, 2020. This is the initial moment of time from which the equations of motion are integrated and the spacecraft begins it’s motion along orbits. The equations of motion in both the CR3BP model and the high-fidelity model are integrated numerically by the classical Runge–Kutta method of the 4th order with a fixed step and absolute accuracy of 108 nondimensional units. All the orbits considered in this paper are adapted to the high-fidelity model with absolute and relative precision of 107 in phase space by the multiple-shooting scheme. Table 3 contains the non-dimensional gravitational parameters of the planets and the Sun and typical gravitational accelerations that they provide with respect to the gravitational acceleration due to Venus (i.e., typical values of jai ðtÞj=jaV ðtÞj, where ai ðtÞ is the acceleration provided by the ith body and aV ðtÞ is the acceleration caused by Venus at the same moment of time t). The acceleration caused by the SRP force is also presented and normalized by aV ðtÞ. The acceleration values are obtained for the three different Lissajous orbits with the following (AY ; AZ ): 6; 500  6; 378 km, 30; 000  95; 545 km, and 100; 000  311; 261 km. In the current research, we take into account the perturbations up to the 9th degree of the non-dimensional perturbation acceleration: SRP and the gravity of all the planets except Uranus, and Neptune are kept. 3.2. Illumination The objective of this subsection is to provide the formulas for calculating the illumination coefficient q given the positions of the spacecraft r and the Sun rS with respect to Venus. This is done by accurately calculating the inter-

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

8

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

section of two sphere caps by means of spherical geometry; the most useful formulas can be found in (Wertz, 1978)4. Denote by RS ¼ 695; 700 km the radius of the Sun and by RV ¼ 6; 051:8 km the radius of Venus. Introduce also an auxiliary parameter c ¼ RV =ðRS  RV Þ. Let DS ¼ jrS  rj be the distance from the spacecraft to the Sun, and DV ¼ jrj be the distance from the spacecraft to Venus. As seen from the spacecraft, the angular radius of the Sun is aS ¼ arcsinðRS =DS Þ, and the angular radius of Venus is aV ¼ arcsinðRV =DV Þ. The angular distance between the centers of the Sun and Venus is h ¼ arccosðn1  n2 Þ where n1 ¼ ðrS  rÞ=DS ; n2 ¼ r=DV . Depending on the relative positions of the Sun, Venus, and the spacecraft, there are four different situations (Wertz, 1978):  Total eclipse: ðDS > 1Þ and ðDS < 1 þ cÞ and ðaV  aS > hÞ In this case, q¼0  Annular eclipse: ðDS > 1 þ cÞ and ðaS  aV > hÞ In this case, q¼1

1  cos aV 1  cos aS

 Partial eclipse: ðDS > 1Þ and ðaS þ aV > hÞ and jaV  aS j < h In this case, p  a1 cos aS  a2 cos aV  a3 q¼1 pð1  cos aS Þ where cos aV  cos aS cos h sin aS sin h cos aS  cos aV cos h a2 ¼ arccos sin aV sin h cos h  cos aS cos aV a3 ¼ arccos sin aS sin aV

a1 ¼ arccos

 No eclipse (all the other cases): q¼1

1 lI ¼ Dt

Z

Dt

qðtÞ dt 0

where Dt ¼ 2 Earth years and the integral is approximated by the trapezoidal method. Since Dt is greater than the synodic period of Venus, all Sun–Venus–spacecraft geometric configurations are taken into account. This time interval also corresponds to 6–7 quasi-revolutions along the Lissajous orbits considered. The integral is calculated on a uniform grid with approximately 50 nodes per revolution. This number leads to precise enough values of the means: increasing the number of nodes to 6000 nodes corrects the values by no more than 0.005%. However, the mean illumination value depends on the starting epoch at which the mission begins. For example, if t ¼ 0 corresponds to May 1, 2020, then for some orbits the mean value is decreased by 13:5%. Nevertheless, the mean illuminations values obtained for Jan 1, 2020 as a starting epoch can be considered as the reference ones. All the Lissajous orbits considered are adapted to the high-fidelity model described above. In Table 4, some Lissajous orbits are listed together with the mean, minimum, and maximum illumination coefficients. In this work, station-keeping of the 1; 900  800 km, 6; 050  3; 400 km, and 13; 000  4; 600 km orbits is analyzed (with the mean illumination levels 0.2, 0.5, and 0.8, respectively). The L2 -family of halo orbits in the Sun–Venus system bifurcate from planar Lyapunov orbits starting from AY ¼ 455; 672 km. Hence, the prevalent part of these orbits lies in the full light. Since there is no reason to choose halo orbits with small AZ in sense of illumination, we consider the orbits with greater AZ meaning that the problem of spacecraft thermal control can be solved in some other way without using shadows by Venus. The halo orbits with higher AZ can also be comparatively simply targeted after the heliocentric transfer part of the trajectory. In this research, we study the halo orbits with AZ of 200; 000 km, 500; 000 km, and 700; 000 km. The region of linearly stable halo orbits (with large AZ ) is also of our interest, but instead of considering single orbits from this region, we study targeting the central manifold associated with these orbits. In Table 5, the parameters of the orbits chosen are summarized. 3.3. Variational equations Let the equations of motion (4) be written in the form

Figs. 6 and 7 show isolines of the mean illumination coefficient q for a variety of Lissajous orbits with the Yamplitude from 1; 600 km to 32; 000 km and the Zamplitude from 1; 000 km to 10; 000 km. Mean values are calculated according to the formula

4 A different methodology for a smoothly varying illumination coefficient is also presented in (Aziz et al., 2019).

x_ ¼ fðt; xÞ

ð5Þ

where x ¼ ½r; v and the initial value be xðt0 Þ ¼ x0 . Denote by Wðt; t0 Þ the state transition matrix @xðtÞ=@xðt0 Þ and by wðtÞ the vector @xðtÞ=@t0 . Then, the variational equations are @W @f ðt; t0 Þ ¼ ðt; xÞWðt; t0 Þ @t @x

ð6Þ

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

9

Fig. 6. Mean illumination coefficient q for different Lissajous orbits.

Fig. 7. Mean illumination coefficient q for small Lissajous orbits.

Table 4 Lissajous orbits and their illumination levels.

Table 5 Orbits studied in this research.

AY  AZ , km

Mean level

Minimum level

Maximum level

Type of the orbit

AY , km

AZ , km

1,900  800 3,200  1,500 4,800  2,000 6,050  3,400 8,000  4,000 9,500  5,000 13,000  4,600

0.20 0.30 0.40 0.50 0.60 0.70 0.80

0.11 0.12 0.12 0.14 0.14 0.15 0.14

0.30 0.42 0.57 0.71 0.85 0.95 1.0

Lissajous Lissajous Lissajous Halo Halo Halo

1,900 6,050 13,000 440,933 572,597 680,007

800 3,400 4,600 200,000 500,000 700,000

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

10

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

dw @f ðtÞ ¼ ðt; xÞwðtÞ dt @x

ð7Þ

Here xðtÞ is the solution to Eq. (4) or Eq. (5) and the initial conditions are Wðt0 ; t0 Þ ¼ I66 ; wðt0 Þ ¼ fðt0 ; x0 Þ where Imm is the m  m identity matrix and Omn is the m  n zero matrix. In practice, the variational equations are integrated together with the equations of motion. Note also that since wðtÞ ¼ Wðt; t0 Þfðt0 ; x0 Þ, the equations for xðtÞ and Wðt; t0 Þ can be integrated only (Hairer et al., 1993). The right-hand side of Eq. (6) is   O33 I33 @f ¼ @x @ v_ =@r O33 where @ v_ l I33 3l rrT ¼ V 3 þ V5 @r jrj jrj N X  li i¼1

I33 jr  ri j

3



3ðr  ri Þðr  ri ÞT jr  ri j

5

! þ

@aSRP @r

and @aSRP CA r  rS @q qCA I33 þ ¼ m jr  rS j3 @r m jr  rS j3 @r



qCA 3ðr  rS Þðr  rS Þ m jr  rS j5

T

The right-hand side of Eq. (7) contains   O31 @f ¼ @t @ v_ =@t where

! N X @ v_ I33 3ðr  ri Þðr  ri ÞT ¼  li  þ vi 5 3 @t jr  ri j jr  ri j i¼1



n X i¼1

li

I33 jri j3



3ri rTi jri j

5

! vi þ

@aSRP @t

and @aSRP CA r  rS @q qCA vS ¼ vS  3 m jr  rS j @rS m jr  rS j3 @t

þ

qCA 3ðr  rS Þðr  rS ÞT vS 5 m jr  rS j

The expressions for @q @q and @r @rS are listed in Appendix A.

4. Station-keeping simulation algorithm In this research, the target point approach is used to keep a trajectory in the vicinity of the reference orbit (Howell and Pernicka, 1993; Go´mez et al., 1998). This method minimizes the deviations from the reference orbit at specified times. The choice of this station-keeping approach among many others is motivated by simplicity of the algorithms used. This approach does not require deep insight into the complex dynamics of the three- or many body systems or knowledge of advanced optimization methods of control theory. In classical works (Howell and Pernicka, 1993; Go´mez et al., 1998), targeting is used to keep the trajectory in a vicinity (a tube) of the nominal orbit. Station-keeping impulses are performed when three conditions hold: the time since last correction is sufficiently large, the deviation from the nominal orbit is sufficiently large, and that deviation is not decreasing. In our work, we use a modification of the target point approach. Station-keeping impulses are always calculated at the target points – the ones that are used for constructing a measure of deviation of the actual spacecraft’s trajectory from the reference orbit. By doing this, we lower the degree of freedom which we have for the maneuver execution time and avoid artificial constraints such as the minimum time between maneuvers (now it is the time between adjacent target points) and the constraints concerning the minimum distance from the nominal orbit at which maneuvers are not allowed to be applied. Instead of waiting for the spacecraft to reach that critical distance, we apply maneuvers at the preliminarily calculated epochs corresponding to target points, no matter what the deviation from the nominal orbit is. Now, let us proceed to the description of the method used in this work. Two strategies are considered: 1) control impulses are always directed along the X-axis of the rotating frame (i.e. along the Sun–Venus line), and 2) control impulses can have any direction (3-axis control). The first strategy simplifies the control implementation in a real mission, though the second one provides tighter station-keeping. In the target point approach, a set of target points must be selected in the reference orbit. The number of target points and time intervals between them influence the station-keeping efficiency. The number of station-keeping maneuvers per one revolution (if a halo orbit is considered) or one quasi-period of the orbit projection on the YZ-plane (in case of a Lissajous orbit) depends on the three-body system, the orbit, constraints on the minimum time between maneuvers, and many other factors and mission constraints. In this research, we take 4 target points per (quasi-) period, though other values are also worth studying. By doing so, we suppose that the time between successive maneuvers satisfies mission requirements. Otherwise, the number of maneuvers or station-keeping scheme should be reconsidered.

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

Denote by t0 < t1 < t2 < . . . < tN p the specified times at which the target points xr ðt0 Þ; xr ðt1 Þ; xr ðt2 Þ; . . . ; xr ðtN p Þ in the reference orbit are chosen; N p is the total number of the target points. Let xðt0 Þ ¼ xr ðt0 Þ þ nins be the initial state of the spacecraft where nins has a Gaussian distribution with zero mean and the covariance matrix " # r2R I33 O33 r2V I33

O33

with some standard deviations rR and rV that define the orbit insertion errors. Further actions depend on the control case. X-axis control: find Dv such that after updating vðt0 Þ vðt0 Þ þ Dv  e the objective J ðDvÞ ¼

m X

jxðti Þ  xr ðti Þj2

i¼1

is minimized; here e is the Sun-to-Venus direction, xðti Þ is the propagated state of the spacecraft at the epoch ti , and m is the number of target points taken into account in the optimization process. In this study, we choose m ¼ 3 target points downstream though this parameter can be tuned for better algorithm performance. 3-axis control: find Dv such that after updating vðt0 Þ vðt0 Þ þ Dv the objective J ðDvÞ ¼

m X

2

jxðti Þ  xr ðti Þj

i¼1

is minimized. The both optimization problems can be effectively solved numerically by the Levenberg–Marquardt method. Usually, 2–3 iterations are needed to get a solution with a relative tolerance of 106 . After the impulse is found, the system uncertainty is simulated. First, the initial state xðt0 Þ is shifted: xðt0 Þ nnav

xðt0 Þ þ nnav " r2 I33 2 N O61 ; r O33

O33 r2v I33

vðt0 Þ

1. Take two unit vectors n1 and n2 such that Dv; n1 and n2 form a right-handed orthogonal system of vectors. 2. Take u ¼ gu , where gu has a uniform distribution over the interval ½0; 2p. 3. Calculate nr ¼ n1 cos u þ n2 sin u. 4. Take a Gaussian random variable na with zero mean and standard deviation ra . 5. Calculate a rotational matrix R that performs rotation around nr by angle a. If the resulting impulse is too small, i.e. ð1 þ nDv ÞDv < Dvmin for some Dvmin , then the impulse is skipped. After that, the trajectory is propagated until the epoch t1 where the deviation dx1 ¼ xðt1 Þ  xr ðt1 Þ is computed and stored. Then, t1 becomes the initial moment of time, xr ðt2 Þ; xr ðt3 Þ, and xr ðt4 Þ are the next target points, and the simulation algorithm repeats again: optimizing the impulse, modeling the navigation errors, applying the impulse augmented with the execution errors, and propagating the trajectory until the epoch t2 . Then, t2 becomes the initial moment of time, and the above process repeats again. This repeating procedure ends at the epoch tN p 3 (when the optimal impulse is calculated based on the deviations from the reference orbit at the epochs tN p 2 ; tN p 1 , and the final tN p ). At this step, one has a Monte Carlo simulation. Another simulation begins with t0 and provides new results, slightly different from the first simulation. Performing a series of Monte Carlo simulations, we analyze the mean station-keeping cost, its variance, and the maximum deviations from the reference orbit. In Table 6, all the parameters of the station-keeping simulation used in this study are summarized. These values approximately correspond to those obtained for the Japan Aerospace Exploration Agency’s Akatsuki mission to Venus (Ryne et al., 2011). Let each of the Monte Carlo simulations be performed over a time span ½t0 ; tN p . The computed values of impulses ðkÞ

are denoted by Dvj for j ¼ 1; . . . ; N p  3 and k ¼ 1; . . . ; K where K is the number of the Monte Carlo simulations. By

#!

where nnav has a multivariate Gaussian distribution with zero mean and the diagonal covariance matrix with the standard deviations rr for the position and rv for the velocity. After that, the impulse augmented with the execution error is applied: vðt0 Þ þ ð1 þ nDv Þ  RDv

where nDv has a Gaussian distribution with zero mean and some standard deviation rDv , and R is a random rotational matrix that models the error in the impulse direction and is calculated by the algorithm:

11

DvðkÞ ¼

N p 3  X 1  ðkÞ  Dvj  tN p  t0 j¼1

Table 6 Parameters of the station-keeping simulation. Parameter

Value

rR rV rr rv 3rDv 3ra Dvmin m K

100 km 10 cm/s 10 km 1 cm/s 0:02 1 1 mm/s 3 200

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

12

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

we denote the annual station-keeping cost obtained in the kth Monte Carlo simulation; here ðtN p  t0 Þ is expressed in Earth years. The set of DvðkÞ forms a sample, and each DvðkÞ is a sample Dv value. The sample values are used to estimate the annual station-keeping cost by the formula Dvsk ¼

p 3 K K N  X X 1X 1  ðkÞ  DvðkÞ ¼ Dvj  K k¼1 KðtN p  t0 Þ k¼1 j¼1

ðkÞ

Drmax ¼ maxjdrj j k;j

5. Results 5.1. Unstable Lissajous and halo orbits

This value is the mean annual station-keeping cost. As a measure of variability of the station-keeping cost, we use an estimator of the standard deviation (SD) rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 XK 2 rDvsk ¼ ðDvðkÞ  Dvsk Þ k¼1 K 1 The position and velocity deviations from the nominal orbit are calculated at the target points on the controlled ðkÞ ðkÞ trajectory and are denoted by drj and dvj ; as previously, j stands for the target point and k is the number of the corresponding Monte Carlo simulation. For each of the orbits, we also calculate the maximum deviation from the reference orbit:

The results of 200 Monte Carlo simulations for different Lissajous and halo orbits are presented in Tables 7 (the Xaxis control) and 8 (the 3-axis control). The results show that the mean annual station-keeping cost for the Lissajous orbits takes values 1:4–1:5 m/s/year and does not depend on the size of the orbit (at least, in the considered cases). In Fig. 8, a typical example of the histogram of the sample Dv; in this particular case, it is a histogram for the Lissajous orbit 6,050  3,400 km. In all the cases, over 98% of the calculated station-keeping cost lie within the range Dvsk 3rDvsk . From Tables 7 and 8 it also follows that the 3-axis control requires higher cost while providing 2–3 times smaller deviations from the reference orbit. Strong dependence of the cost on the orbit insertion and navigation errors is also observed, as seen from Tables 9

Table 7 Mean annual station-keeping cost, its standard deviation,and the maximum deviation from the reference orbit; the X-axis control. Type of the orbit

AY  AZ , km

Dvsk , m/s/year

rDvsk , m/s/year

Max deviation, km

Lissajous Lissajous Lissajous Halo Halo Halo

1,900  800 6,050  3,400 13,000  4,600 440,933  200,000 572,597  500,000 680,007  700,000

1.140640 1.151472 1.144088 1.038363 0.908447 0.766427

0.090263 0.100858 0.095452 0.100828 0.100944 0.077872

740.230150 627.309827 612.612098 783.863879 961.489560 1,295.723246

Fig. 8. Histogram of sample Dv for the 6,050  3,400 km Lissajous orbit; the X-axis control.

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

13

Table 8 Mean annual station-keeping cost, its standard deviation, and the maximum deviation from the reference orbit; the 3-axis control. Type of the orbit

AY  AZ , km

Dvsk , m/s/year

rDvsk , m/s/year

Max deviation, km

Lissajous Lissajous Lissajous Halo Halo Halo

1,900  800 6,050  3,400 13,000  4,600 440,933  200,000 572,597  500,000 680,007  700,000

1.213267 1.216384 1.216128 1.138523 1.004817 0.871449

0.091439 0.085488 0.086184 0.093358 0.080704 0.075924

340.186621 374.733444 383.862144 404.384913 369.591763 354.378774

Table 9 Mean annual station-keeping cost, its standard deviation,and the maximum deviation from the reference Lissajous orbit 6; 050  3; 400 km for different levels of the orbit insertion and navigational uncertainties; the X-axis control. Uncertainty level

Dvsk , m/s/year

rDvsk , m/s/year

Max deviation, km

rR ¼ 20 km, rV ¼ 2 cm/s rr ¼ 2 km, rv ¼ 0:2 cm/s

0.230678

0.019658

152.120344

rR ¼ 50 km, rV ¼ 5 cm/s rr ¼ 5 km, rv ¼ 0:5 cm/s

0.574301

0.045031

344.639002

rR ¼ 100 km, rV ¼ 10 cm/s rr ¼ 10 km, rv ¼ 1 cm/s

1.151472

0.100858

627.309827

rR ¼ 200 km, rV ¼ 20 cm/s rr ¼ 20 km, rv ¼ 2 cm/s

2.319981

0.193106

1,489.591905

rR ¼ 500 km, rV ¼ 50 cm/s rr ¼ 50 km, rv ¼ 5 cm/s

5.711677

0.456067

2,975.897634

and 10 for the Lissajous orbit 6; 050  3; 400. Multiplying rR ; rV ; rr , and rv by q times changes the station-keeping cost, its SD, and the maximum deviations from the orbit by the same q times. The same observation holds true for the considered halo orbits, too. From Tables 9 and 10 it is also seen that multiplying rR ; rV ; rr , and rv by more than 2 times loses the control: the real trajectory may appear substantially far from the reference orbit, and the mission can fail. In a real mission, the levels of uncertainty rR ; rV ; rr , and rv are recommended to be not larger than in our central case (see Table 6). From Tables 7 and 8, it follows that larger halo orbits have smaller station-keeping costs: for the halo orbit with AZ ¼ 200; 000 km, the mean annual cost is 1:04 m/s/year; for AZ ¼ 500; 000 km, it is 0:91 m/s/year, and for

AZ ¼ 700; 000 km, it is 0:77 m/s/year. This result is expected since larger halo orbits are ‘‘more stable” in a sense that their Floquet multipliers are closer to the unit circle in the complex plane. 5.2. Central manifold targeting Previously, a family of linearly stable halo orbits between Venus and L2 were localized. Since these orbits belong to the central manifold, it is reasonable to target the manifold itself when constructing a transfer from the Earth to the Sun–Venus system. In this section, approximate bounds of the central manifold and permissible position and velocity errors which lead to bounded trajectories around L2 are determined.

Table 10 Mean annual station-keeping cost, its standard deviation, and the maximum deviation from the reference Lissajous orbit 6; 050  3; 400 km for different levels of the orbit insertion and navigational uncertainties; the 3-axis control. Uncertainty level

Dvsk , m/s/year

rDvsk , m/s/year

Max deviation, km

rR ¼ 20 km, rV ¼ 2 cm/s rr ¼ 2 km, rv ¼ 0:2 cm/s

0.247905

0.019124

74.197181

rR ¼ 50 km, rV ¼ 5 cm/s rr ¼ 5 km, rv ¼ 0:5 cm/s

0.611229

0.046707

175.168375

rR ¼ 100 km, rV ¼ 10 cm/s rr ¼ 10 km, rv ¼ 1 cm/s

1.216384

0.085488

374.733444

rR ¼ 200 km, rV ¼ 20 cm/s rr ¼ 20 km, rv ¼ 2 cm/s

2.436629

0.188412

870.534192

rR ¼ 500 km, rV ¼ 50 cm/s rr ¼ 50 km, rv ¼ 5 cm/s

6.504855

5.071890

95,589.376245

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

14

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

In the CR3BP model, halo orbits cross orthogonally to the XZ-plane and are uniquely determined by the X ; Z, and Y_ coordinates of the rotating coordinate frame. For the linearly stable halo orbits, these coordinates vary along the family in the following ranges: 580. . .650 thousand km along the X-axis from Venus 1,239. . .1,246 thousand km along the Z-axis 410. . .450 m/s along the Y_ -axis It is of interest whether all the phase states on the XZplane that belong to the region indicated above lead to bounded motion. Technically, by bounded motion we mean a trajectory that stays within the range 1:0 6 X 6 1:007 of non-dimensional units during the five Earth years; this interval is selected to a certain extent arbitrarily, other criteria are also possible. The investigation of this question gives a positive answer, see Fig. 9. The red line represents a family of the linearly stable halo orbits, the color cells indicate a difference in the maximum and minimum possible Y_ that result in a bounded trajectory. It follows that there is a region of size 30,000  7,000 km (610–640 thousand km along the X-axis from Venus, 1,239–1,246 km along the Z-axis) on the XZ-plane targeting at which with a possible velocity error of 1:5 m/s gives bounded motion. Some regions—for example, 580–590 thousand km along the X axis from Venus and 1,239– 1,240 thousand km along the Z axis—correspond to bounded motion as well, but a possible error in Y_ reduces to below 0:25 m/s. In the high-fidelity model, due to gravitational perturbations of the planets and, mostly, because of the eccentricity of the Venusian orbit and the SRP force, halo orbits no longer cross orthogonally to the XZ-plane and are no longer closed. Following the same idea as for the CR3BP

Fig. 10. Difference between the maximum and the minimum values of Y_ that lead to bounded motion during the two Earth years. The X values are measured from Venus.

model of motion, we investigate boundedness of the motion by propagating trajectories from the following region: 550. . .650 thousand km along the X-axis from Venus 1,239. . .1,246 thousand km along the Z-axis 450. . .410 m/s along the Y_ -axis assuming that Y ¼ X_ ¼ Z_ ¼ 0. Although X_ – 0 and Z_ – 0 for halo orbits at the moment of intersection, we take X_ ¼ Z_ ¼ 0 for simplicity because the boundedness of a trajectory is of only interest. In the high-fidelity model, by bounded motion we mean a trajectory that stays within the range 1:0 6 X 6 1:007 of non-dimensional units during two Earth years; a stronger flight constraint results in too few bounded trajectories. The results are shown in Fig. 10. It is seen that there is a region of size 60,000  7,000 km (550–610 thousand km along the Xaxis from Venus, 1,239–1,246 km along the Z-axis) on the XZ-plane targeting at which with a possible velocity error of 0:75 m/s gives bounded motion. 6. Conclusion

Fig. 9. Difference between the maximum and the minimum values of Y_ that lead to bounded motion during the five Earth years; the red line represent the points on the XZ-plane belonging to the linearly stable halo orbits. The X values are measured from Venus.

Station-keeping costs required to stabilize the small Lissajous orbits with mean illumination level in the range from 0.1 to 1.0 do not depend on the size of the orbits and equal 1:14 0:3 m/s per Earth year. Station-keeping costs required for halo orbits are slightly less than for the Lissajous orbits and decreases as the amplitude increases for the cases considered in this work. It is also highlighted that the most part of any halo orbit lie in the full light. If mission requirements state that a partial shadow is needed most of the time (for example, to keep an admissible level of temperature onboard), the small Lissajous orbits should then be considered.

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

It is observed that the station-keeping cost depends linearly on the standard deviations that characterize the orbit insertion and navigation errors. It is recommended to have orbit insertion errors below 100 km in position and 10 cm/s in velocity, and navigation errors below 10 km in position and 1 cm/s in velocity; otherwise, deviations from the reference orbit become large and the probability of mission failure increases. The results also show that the station-keeping costs for the 3-axis control are by 0.07–0.1 m/s/year (6–10%) larger than the ones for the X-axis control. As expected, the maximum deviations from the reference orbit are lower in the first case. Existence of the linearly stable halo orbits and the associated central manifold provides an option to place the spacecraft in a bounded trajectory that does not require station-keeping if some manifold insertion constraints are satisfied. In this study, we consider a projection of the central manifold onto the plane that is perpendicular to the Venus orbit and contains the Sun–Venus line. Varying points on this plane and propagating trajectories with different orthogonal velocities, bounded trajectories are registered. The calculations show that the two-year-bounded motion can be achieved provided we target some region of size 60; 000  7; 000 km with an accuracy of 0:75 m/s in velocity. A narrower region that requires an accuracy of 1:75 m/s in velocity is also found and depicted. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

15

@aS 1 RS ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðrS  rÞT @rS 1  R2 =D2 DS S

S

@aV 1 RV ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 rT @r 2 2 DV 1  R =D V

V

@aV ¼ O31 @rS

 @h 1 @n1 @n2 ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 þ n1 @r @r @r 2 1  ðn1  n2 Þ  @h 1 @n1 @n2 ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 þ n1 @rS @rS @rS 2 1  ðn1  n2 Þ @n1 I33 ðr  rS Þðr  rS ÞT ¼ þ @r DS D3S @n2 I33 rrT ¼ þ 3 @r DV DV @n1 I33 ðr  rS Þðr  rS Þ ¼  @rS DS D3S

T

@n2 ¼ O33 @rS Derivatives with respect to the angles aS ; aV , and h are case-dependent. In case of annual eclipse @q sin aV ¼ @aV 1  cos aS @q sin aS ¼ ð1  qÞ @aS 1  cos aS @q ¼0 @h

Acknowledgments

In case of partial eclipse, the derivatives with respect to the angles are the following:

This work is fully supported by the Russian Science Foundation grant 19-11-00256.

@q a1 sin aS  sin aS ðp  a2 cos aV  a3 Þ ¼ @aS pð1  cos aS Þ2 cos aS @a1 þ pð1  cos aS Þ @aS

In case of total eclipse or no eclipse,

Appendix A. @q @q ¼ ¼0 @r @rS

þ

Otherwise, @q @q @aS @q @aV @q @h ¼ þ þ @r @aS @r @aV @r @h @r @q @q @aS @q @aV @q @h ¼ þ þ @rS @aS @rS @aV @rS @h @rS where in both annual and partial eclipse cases @aS 1 RS T ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðrS  rÞ @r 2 2 DS 1  R =D S

S

cos aV @a2 1 @a3 þ pð1  cos aS Þ @aS pð1  cos aS Þ @aS

@q a2 sin aV cos aS @a1 þ ¼ @aV pð1  cos aS Þ pð1  cos aS Þ @aV

þ

cos aV @a2 1 @a3 þ pð1  cos aS Þ @aV pð1  cos aS Þ @aV

@q cos aS @a1 cos aV @a2 1 ¼ þ þ @h pð1  cos aS Þ @h pð1  cos aS Þ @h pð1  cos aS Þ @a3  @h

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022

16

M. Shirobokov et al. / Advances in Space Research xxx (2020) xxx–xxx

@a1 cos a3 sin aV ¼ @aS j sin a1 j sin aS sin h @a1 1 sin aV ¼ j sin a1 j sin aS sin h @aV @a1 cos a2 sin aV ¼ @h j sin a1 j sin aS sin h @a2 1 sin aS ¼ @aS j sin a2 j sin aV sin h @a2 cos a3 sin aS ¼ @aV j sin a2 j sin aV sin h @a2 cos a1 sin aS ¼ @h j sin a2 j sin aV sin h @a3 cos a1 sin h ¼ @aS j sin a3 j sin aS sin aV @a3 cos a2 sin h ¼ @aV j sin a3 j sin aS sin aV @a3 1 sin h ¼ j sin a3 j sin aS sin aV @h Appendix B. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j. asr.2019.12.022. References Akhmetshin, R.Z., 2013. Space patrol: Variants of the optical barrier scheme. Cosm. Res. 51 (4), 241–253. https://doi.org/10.1134/ S0010952513040011. Ascher, U.M., Mattheij, R.M.M., Russel, R.D., 1995. Numerical solution of boundary value problems for ordinary differential equations. SIAM 170–180. https://doi.org/10.1137/1.9781611971231. Aziz, J., Scheeres, D., Parker, J., Englander, J., 2019. A smoothed eclipse model for solar electric propulsion trajectory optimization. Trans. Japan Soc. Aero. Space Sci. 17 (2), 181–188. https://doi.org/ 10.2322/tastj.17.181. Bodrova, Y.S., Raikunov, K.G., 2017. Approaches to selecting rational variants of the ballistic construction of space telescopes for operative detection and determination of physicochemical properties of asteroids unfavorable for observation from the Earth. Eng. J.: Sci. Innov. 7, 16. https://doi.org/10.18698/2308-6033-2017-7-1640. Council, N.R., 2010. Defending Planet Earth: Near-Earth-Object Surveys and Hazard Mitigation Strategies. The National Academies Press, Washington, DC. https://doi.org/10.17226/12842. Dei Tos, D.A., Topputo, F., 2017. Trajectory refinement of three-body orbits in the real solar system model. Adv. Space Res. 59 (8), 2117– 2132. https://doi.org/10.1016/j.asr.2017.01.039. Dunham, D.W., Genova, A.L., 2010. Using Venus for locating space observatories to discover potentially hazardous asteroids. Cosm. Res. 48 (5), 424–429. https://doi.org/10.1134/S0010952510050084.

Dunham, D.W., Reitsema, H.J., Lu, E.T., Arentz, R., Linfield, R., Chapman, C., Farquhar, R.W., Ledkov, A.A., Eismont, N.A., Chumachenko, E., 2013. A concept for providing warning of earth impacts by small asteroids. Sol. System Res. 47 (4), 315–324. https:// doi.org/10.1134/S0038094613040096. Gao, S., Zhou, W., Zhang, L., Liang, W., Liu, D., Zhang, H., 2019. Trajectory design and flight results for Chang’e 4-relay satellite. Scient. Sin. Tech. 49 (2), 156–165. https://doi.org/10.1360/N092018-00393. ` ., Simo´, C., Masdemont, J.J., 2001. Dynamics and Go´mez, G., Jorba, A mission design near libration points: volume 3, advanced methods for collinear points. World Sci., 33–52 https://doi.org/10.1142/4337. Go´mez, G., Llibre, J., Matı´nez, R., Simo´, C., 2001. Dynamics and mission design near libration points: volume 1, fundamentals: the case of collinear libration points. World Sci., 62–68 https://doi.org/10.1142/ 4402. Go´mez, G., Howell, K.C., Simo´, C., Masdemont, J.J., 1998. StationKeeping Strategies for Translunar Libration Point Orbits. In: AAS/ AIAA Spaceflight Mechanics Conference, Monterey, CA, United States. Paper AAS 98-168, pp. 1–20. Gold, R.E., 2001. SHIELD: a comprehensive earth-protection architecture. Adv. Space Res. 28 (8), 1149–1158. https://doi.org/10.1016/ S0273-1177(01)00496-3. Greenstreet, S., Ngo, H., Gladman, B., 2012. The orbital distribution of near-earth objects inside earth’s orbit. Icarus 217 (1), 355–366. https:// doi.org/10.1016/j.icarus.2011.11.010. Hairer, E., Nørsett, S.P., Wanner, G., 1993. Solving Ordinary Differential Equations I. Nonstiff Problems. Springer-Verlag, Berlin Heidelberg, p. 95. https://doi.org/10.1007/978-3-540-78862-1. Howell, K.C., Pernicka, H.J., 1993. Station-keeping method for libration point trajectories. J. Guid. Control Dyn. 16 (1), 151–159. https://doi. org/10.2514/3.11440. Howell, K.C., 1984. Three-dimensional, periodic, ‘halo’ orbits. Celest. Mech. 32 (1), 56–59. https://doi.org/10.1007/BF01358403. Jedicke, R., Morbidelli, A., Spahr, T., Petit, J.-M., Bottke, W.F., 2003. Earth and space-based NEO survey simulations: prospects for achieving the spaceguard goal. Icarus 161 (1), 17–33. https://doi.org/10.1016/ S0019-1035(02)00026-X. Lu, E.T., Reitsema, H., Troeltzsch, J., Hubbard, S., 2013. The B612 foundation sentinel space telescope. New Space 1 (1), 42–45. https:// doi.org/10.1089/space.2013.1500. Meyer, K.R., Hall, G.R., Offin, D., 2008. Introduction to Hamiltonian Dynamical Systems and the N-body problem. Springer Science & Business Media, p. 200. https://doi.org/10.1007/978-0-387-09724-4. Perko, L., 2001. Differential Equations and Dynamical Systems. SpringerVerlag, New York, p. 225. https://doi.org/10.1007/978-1-4613-0003-8. Roberts, C.E., Case, S., Reagoso, J., 2016. Lissajous orbit control for the deep space climate observatory Sun-Earth L1 libration point mission. Adv. Astronaut. Sci. 165, 2615–2634. Ryne, M.S., Mottinger, N.A., Broschart, S.B., You, T.-H., Higa, E., Helfrich, C., Berry, D., 2011. Navigation support at JPL for the JAXA Akatsuki (PLANET-C) Venus Orbiter Mission. Adv. Astronaut. Sci. 142, 1799–1818. Shirobokov, M., Trofimov, S., Ovchinnikov, M., 2017. Survey of stationkeeping techniques for libration point orbits. J. Guid. Control Dyn. 40 (5), 1085–1105. https://doi.org/10.2514/1.G001850. Szebehely, V., 1967. Theory of Orbits: The Restricted Problem of Three Bodies. Academic Press Inc., New York, pp. 556–558. Wertz, J.R., 1978. Spacecraft attitude determination and control. Kluwer, 727–736, p. 76. https://doi.org/10.1007/978-94-009-9907-7.

Please cite this article as: M. Shirobokov, S. Trofimov and M. Ovchinnikov, On the design of a space telescope orbit around the Sun–Venus L2 point, Advances in Space Research, https://doi.org/10.1016/j.asr.2019.12.022