A fast Chebyshev polynomial method for calculating asteroid gravitational fields using space partitioning and cosine sampling

Journal Pre-proofs A Fast Chebyshev Polynomial Method for Calculating Asteroid Gravitational Fields using Space Partitioning and Cosine Sampling Hongwei Yang, Shuang Li, Jun Sun PII: DOI: Reference:

S0273-1177(19)30792-6 https://doi.org/10.1016/j.asr.2019.11.001 JASR 14525

To appear in:

Advances in Space Research

Received Date: Revised Date: Accepted Date:

16 June 2019 30 October 2019 1 November 2019

Please cite this article as: Yang, H., Li, S., Sun, J., A Fast Chebyshev Polynomial Method for Calculating Asteroid Gravitational Fields using Space Partitioning and Cosine Sampling, Advances in Space Research (2019), doi: https://doi.org/10.1016/j.asr.2019.11.001

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 COSPAR. Published by Elsevier Ltd. All rights reserved.

A Fast Chebyshev Polynomial Method for Calculating Asteroid Gravitational Fields using Space Partitioning and Cosine Sampling Hongwei Yang1, 2, Shuang Li1, 2*, Jun Sun3 Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China 2. College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China 3. Shanghai Aerospace Control Technology Institute, Shanghai, 201109, China

1. State

This paper presents a computationally fast method for solving gravitational accelerations near irregularly-shaped asteroids. This method is based on analytical three-dimensional Chebyshev polynomial approximation of the polyhedral gravity. For the purpose of improving the approximation accuracy, space partitioning schemes based on practical flight zones is used to avoid interpolation the whole space around the target asteroid. Specifically, a minimum ellipsoid close to the asteroid surface is defined to select the space for surrounding trajectories with safe distance and a cone connected to the surface is defined to select the space for descent trajectories. Moreover, interpolation points are sampled in a cosine sampling fashion according to the Chebyshev-Gauss-Lobatto nodes and a radial adaption technique. The performance of different space partitioning schemes is analyzed. The effectiveness of the proposed method is validated through simulations of solving gravitational accelerations at the test points near different shaped asteroids 1996 HW1, 433 Eros, 25143 Itokawa and 101955 Bennu. Keywords: Asteroid; Gravitational field; Chebyshev polynomials; Space partitioning; Chebyshev-Gauss-Lobatto node

I.

Introduction

The two ongoing asteroid missions, Hayabusa2 (Tsuda et al., 2013) and OSIRIS-REx (Lauretta and OSIRIS-REx Team, 2012) have attracted great attentions around the world. The new detection results indicate the significant

*

Professor, Email: [email protected] (Corresponding author) 1

scientific value of asteroid exploration (Kitazato et al., 2019; Hamilton et al., 2019; Lauretta et al., 2019). For the purpose of obtaining high resolution data via observation or sampling in asteroid missions, close proximity operations, such as hovering (Broschart and Scheeres, 2005; Zeng et al., 2016; Yang et al., 2018; D'Ambrosio et al. 2019), orbiting (Yu and Baoyin, 2012a; Chanut et al., 2015; Surovik and Scheeres, 2015), landing (Furfaro et al., 2013; Yang et al. 2017; Bazzocchi and Emami, 2018; Pinson and Lu, 2018; Yang et al., 2019) and surface hopping (Liu et al., 2017), are essential. One of the biggest challenges for close proximity operations of asteroid exploration arises from the irregular gravitational fields of irregularly-shaped asteroids. Thus far, several effective and useful methods have been proposed for modelling gravitational fields of asteroids. Among these methods, the classical polyhedral method proposed by Werner (1994) and Werner and Scheeres (1996) is usually regarded as a high-fidelity method for describing the irregular gravitational fields. This method has been widely used in the areas of orbital dynamics analysis (Yu and Baoyin, 2012b; Yu and Baoyin, 2013; Wang et al. 2014; Jiang et al. 2018; Yang et al., 2018), orbit design (Yu and Baoyin, 2012a; Chanut et al., 2015; Surovik and Scheeres, 2015), trajectory optimization (Pinson and Lu, 2018; Yang et al., 2019) and orbit control design (Broschart and Scheeres, 2005). Although the accuracy of the polyhedral method is high, its computational efficiency is low which limits its application for onboard computation. However, because of the existence of communication delays between the ground and asteroids, autonomous guidance and maneuver planning (Huang et al. 2004; Surovik and Scheeres, 2014) which require the ability of onboard computation of gravitational accelerations become necessary for close proximity. Moreover, onboard optimal landing trajectory design based on methods such as convex optimization based methods (Yang et al. 2017; Pinson and Lu, 2018) also demands the ability of fast computation of gravitational accelerations. Although these methods have already made the efficient convex optimization be applicable for trajectory optimization with irregular gravitational fields, slow computation of gravitational accelerations in orbital propagations greatly restricts onboard implementation of these methods. The spherical harmonic method is a classical modeling method for non-spherical gravitational fields and usually has high computational efficiency. This method, however, usually has the convergence problem inside the so-called Brillouin sphere (Scheeres et al., 2000; Takahashi and Scheeres, 2013; Takahashi and Scheeres, 2014; Hu and Ji, 2017). Due to this fact, trajectory design for asteroid landing using the spherical harmonic method may be not accurate enough (Takahashi and Scheeres, 2013). Takahashi and Scheeres (2014) developed the spherical harmonic method to deal with the convergence problem inside the Brillouin sphere. However, the complexity of this method may restrict

2

the practical implementation (Hu and Ji, 2017). Besides the spherical harmonic method, another common type of modeling methods with high computational efficiency is the simplified model methods. In a simplified model method, an asteroid is modeled by a common geometry or some mass points such as a straight segment (Riaguas and Elipe, 2001), an ellipsoid (Scheeres, 1994), a cube (Liu et al. 2011), a dumbbell (Li et al. 2013; Li et al. 2017), a dipole (Zeng et al. 2015) and a dipole segment (Zeng et al. 2018). These models no longer have the convergence problem as the spherical harmonic method. Meanwhile, these models is suitable for analyze the characteristics of gravitational fields. But, the accuracy, especially in the region close to the surface, is low. An alternative approach to lower the computational effort is through simplifying polyhedral models. A representative method is the mascon model proposed by Chanut et al. (2015). Their method generates the data of the mascons using the data of an asteroid polyhedral model and gives analytical expressions for the gravitational field. Compared with the polyhedral method, a speedup can be obtained by the mascon method. Burov et al. (2019) further simplified the mascon method with the help of the k-means clustering method. In their method, only four mass points is needed for describing the gravitational fields of asteroids. However, these methods cannot accurately describe the gravitational fields near the surface. Pearl et al. (2018) proposes a hybrid gravity model method which calculates the landing region by the original polyhedral method and other regions by an approximate method. This method makes a trade-off between the accuracy and the computational speed. However, this method is proposed only for the case when a spacecraft lands on a desired region. Recently, researchers also proposed techniques to fit the polyhedral gravitational fields by new easily calculated functions and then calculate the gravitational accelerations rapidly using the corresponding new functions (Hu and Ji, 2017; Gao and Wenta, 2019). Among these techniques, the Chebyshev polynomial interpolation is a classical technique in astrodynamics. As for computation of gravitational fields, it was first used to improve the efficiency of Earth gravity calculations in last century (Smith and Lyubomirsky, 1981). Junkins et al. (2012) derived the analytical least square solution for the coefficients and also applied the Chebyshev polynomials for approximating the Earth’s gravity. For the purpose of lowering the required sample points and improving the interpolation accuracy, the Chebyshev-Gauss-Lobatto (CGL) nodes and the radial adaption technique which uses cosine sampling are employed in their study. Besides, the interpolation space is divided into two subspaces by a sphere close to the Earth’s surface when interpolations are implemented. Compared with interpolation of the Earth’s gravity, interpolation of the asteroids’ gravity is more challenge due to their irregularly shapes and the associated irregular gravitational fields. Interpolation

3

schemes which are suitable for asteroids are required. Recently, Hu and Ji (2017) first used the Chebyshev polynomials for calculating the gravity of irregularly-shaped asteroids. In their study, a space partitioning scheme which divides the space around the target asteroid into many small cells is used. The sizes of the cells become larger as the distance to the asteroid surface increases. Then, the interpolations are conducted in every cell. The computational speed can be greatly improved and the approximating accuracy can be extremely high. However, the required storage for the coefficients can be tens of megabytes due to a large number of sampling cells. In the work of Zeng et al. (2018), a socalled inner box which contains the target asteroid is proposed to partition the sampling area when solving the parameters of the simplified model. According to their results, such a space partitioning scheme is very effective for lowering the approximating error. However, the inner box scheme uses rectangular coordinates which are not suitable for radial adaptive sampling. This paper aims to propose a fast Chebyshev polynomial method with low storage coefficients for calculating the gravity of irregularly-shaped asteroids. For the purpose of accurate interpolation in a large region with low order polynomials, the Chebyshev-Gauss-Lobatto (CGL) nodes and the radial adaption technique (Junkins et al. 2012) are adopted. Importantly, two new kind of space partitioning schemes are proposed to replace the sphere based space partitioning scheme (Junkins et al. 2012) for accurate interpolations. The space partitioning schemes are proposed according to two practical types of close proximity operations. As for orbiting around an asteroid, considering the spacecraft needs to keep a distance from the surface for safety, a minimum ellipsoid close to the asteroid surface is defined and ellipsoid based space partitioning schemes are used. Because the space very close to the surface is excluded, high precision is achievable in interpolation without dividing the space into small cells. As for asteroid landing, considering the cone region constraint for the descent trajectory, a cone connected to the surface is defined to select the interpolation space. Although the target space contains points very close to the surface, it is easy to achievable high precision interpolation due to the relatively small size of the interpolation space. The rest of the paper is organized as follows. In Section II, the polyhedral method for calculating the gravitational fields of irregularly shaped asteroids is briefly introduced. In Section III, Chebyshev polynomial based threedimensional approximating function is presented and the analytical least square solution for the coefficients is introduced. Then, simplified manipulations are proposed for reducing the storage of matrices caused by Kronecker matrix product operations. Section IV describes the space partitioning schemes, variables for interpolation and an algorithm for the Chebyshev polynomial method. In Section V, the performance of difference space partitioning

4

schemes are compared and the effectiveness of the proposed method is analyzed. In Section V, the proposed method is further applied to three more irregularly-shaped asteroids for validating its effectiveness. Section VI concludes this paper.

II.

Polyhedral Method for Gravitational fields of irregularly shaped asteroids

We consider the problem of calculating gravitational fields of irregularly shaped asteroids. The polyhedral method (Werner, 1994; Werner and Scheeres, 1996) is a classical method for computation of gravitational fields of irregularly shaped asteroids. This method does not have the convergence problem inside the Brillouin sphere and is usually regarded as a high-fidelity gravity method for irregularly shaped asteroids. In this work, the polyhedral method is adopted to generate the gravitational accelerations at the sample points for determining the coefficients of the Chebyshev polynomials. For the calculation of the gravitational field near an asteroid, an asteroid body-fixed frame is built, as shown in Fig 1, where the asteroid is plotted using the shape data of 1996 HW1 (Magri et al. 2017). The origin of the frame o-xyz coincides with the mass center of the asteroid, z-axis is aligned with the maximum principle axis, x-axis is aligned with the minimum principle axis, and y-axis completes the right-handed system.

z

o

y

x

Fig 1. Schematic of asteroid body-fixed frame In the body-fixed frame, the gravitational acceleration computed by the polyhedral method is (Werner and Scheeres, 1996) g  r   U  G



Ee  re  Le  G

eedges



f  faces

5

F f  rf   f

(1)

where r = [x, y, z]T is the position vector of a field point P near the asteroid, G is the gravitational constant, σ is the bulk density and ra (a = e or f) is a vector from P to any point on an edge or face. The rest variables in Eq. (1) are defined as follows (Werner and Scheeres, 1996): F f  nˆ f nˆ f

(2)

Ee  nˆ f1 nˆ12f1  nˆ f2 nˆ 21f2

(3)

Le  ln

 f  2 arctan

r1  r2  r12 r1  r2  r12

(4)

r1  r2  r3 r1r2 r3  r1  r2  r3   r2  r1  r3   r3  r1  r2 

(5)

where Ff and Ee are dyads. The variables in Eqs. (2) - (5) are illustrated in Fig 2 where two adjacent tetrahedrons are drawn. As shown in this figure, nˆ f is a face normal vector, nˆ12f1 and nˆ 21f2 are edge normal vectors, and r1, r2 and r3 are position vectors of the vertices of a surface triangular.

3

1 f1

P

2

f2

o Fig 2. Illustration of the polyhedral method The polyhedral method can be applied to any asteroid with given polyhedron shape data. However, the computational efficiency of the polyhedral method is relatively low, especially for the asteroids with refined shape models. Therefore, this method is limited to onboard implementation. In order to improve the computational speed, we propose to use low scale Chebyshev polynomials for fitting the gravitational fields obtained by the polyhedral method and calculating the gravitational accelerations at desired points. The method using Chebyshev polynomials will be presented in next sections.

6

III.

Overview on Three Dimensional Chebyshev Polynomial Approximation

The gravitational acceleration (1) is a three dimensional function and has three independent variables. This section briefly will introduce the method proposed by Junkins et al. (2012) for approximation in three dimensions with Chebyshev polynomials. This method has been applied to approximating the Earth’s gravity and shown high accuracy. Then, an approach for reducing the storage of matrix manipulation will be presented. A. Approximating functions with three variables by Chebyshev polynomials In the method of Junkins et al. (2012), the Chebyshev polynomials of the first kind, which are orthogonal polynomials, are chosen as independent basis functions. The definition for the Chebyshev polynomials T can be obtained by the following recursive equation: T0    1, T1         Tk 1    2 Tk    Tk 1   , k  1

(6)

where ξ is a non-dimensional variable which ranges from -1 to 1. The coordinate variables of the gravitational acceleration function can be transformed to non-dimensional variables as

x

 1 2

 xmax  xmin   xmin ,

xmax  xmin

(7)

where the subscript max and min represent the maximum and minimum values, respectively. In this study, the selected coordinate variables will be described in Section IV. B. The maximum and minimum values of the coordinate variables depend on the defined coordinate system and the size of the target interpolation space. According to Hu and Ji (2017) and Junkins et al. (2012), a function f with three independent variables can be approximated by N

N

N

f  , ,     a pqr Tp  Tq   Tr  

(8)

p 0 q 0 r 0

where Nξ, Nη and Nζ are degrees of the Chebyshev polynomials, and a pqr are coefficients which needs to be determined using the values of the function at the sample points. By denoting the function vector as f at the sample points i , j ,  k  where i  0, 1,  , M  , j  0, 1,  , M  and k  0, 1,  , M  , the analytical least square solution for

the corresponding coefficient vector a becomes (Junkins et al., 2012) aC f

(9)

7

f  W 1/ 2 f

(10)

C  C  C  C

(11)

W  W  W  W

(12)

where

The operator  in Eqs. (11) and (12) represents the Kronecker matrix product. The matrices C and W have the following expressions (Junkins et al., 2012): 1 1 1 1 1 W =diag  , , ,  , ,  4 2 2 2 4  M  1 M  1 

 2 T0  0  T0 1    2  2T   2T1 1  1 0 1  C     M   2TN 1  0  2TN 1 1      2T   2TN 1  0 N 

(13)

2 T0  M  2

 





T0  M  1

M  1

1



   2T    

  2T  



1

   

M

   

 2TN 1  M  1

2TN 1  M 



2TN

2TN

M  1

M

  M   N      

(14)

or      1  C  M      

2 T0  M  2

 





T0 1 



2T1  0 

2T1 1  



T0  M  1

M  1

1



2TN 1  0 

2TN 1 1  

2TN 1  M  1

2 TN  0  2 

TN 1 

TN  M  1









   2T    

  2T  

2 T0  0  2

1



M

   

2TN 1  M  2 TN 2 

M

  M   N       

(15)

In addition, W and W have the same form as W while C and C have the same form as C . Thus, their expressions are omitted herein. B. Reducing storage of matrix manipulation Due to the use of the Kronecker matrix product operations in Eqs. (11) and (12), the size of C and W increases quickly with M. To reduce the storage of matrix manipulation, the following approach is used in our code. First, a vector composed of the non-zero elements of the matrix W is used instead of W, as

8

W   [W ,1W ,1 , W ,1W ,2 , , W ,1W , M 1 W ,2W ,1 , W ,2W ,2 , , W ,2W , M 1    W , M 1 M 1W ,1 , W , M 1 M 1W ,2 , , W , M 1 M 1W , M 1 ]T       

(16)

where

W   [W ,1W ,1 , W ,1W ,2 , , W ,2W M W ,2W ,1 , W ,2W ,2 , , W ,2W , M 1   W , M 1W ,1 , W , M 1W ,2 , , W , M 1W , M 1 ]T    

(17)

Then, the element of the vector f is computed by f  i   W  i  f  i  , i  1, 2, ,  M   1 M   1 M   1

(18)

where f  i  is the element of f which is f  [ f (1 ,1 ,  1 ), , f (1 ,1 ,  M  1 ) f (1 , M  1 ,  1 ), , f (1 , M  1 ,  M  1 )  f ( M  1 ,1 ,  1 ), , f ( M  1 ,1 ,  M  1 ) f ( M  1 , M  1 ,  1 ), , f ( M  1 , M  1 ,  M  1 )]T

(19)

Second, the elements of a are computed by the following equation to avoid using the matrix C: M  M

a  kmin : kmax    C  i  1, p  1C  j  1, q  1 C f  lmin : lmax 

(20)

p 0 q 0

where i  0, 1,  , N , j  0, 1,  , N for each i , and kmin, kmax, lmin and lmax are defined as kmin  i  N  1 N  1  j  N  1  1, kmax  i  N  1 N  1   j  1  N  1

(21)

lmin  p  M   1 M   1  q  M   1  1, lmax  p  M   1 M   1   q  1  M   1

(22)

IV.

Calculation of Gravitational Accelerations using Chebyshev Polynomials

The previous section has already presented the approximation equation for three dimensional functions using Chebyshev polynomials. In this section, space partitioning schemes will be proposed and an algorithm for calculating gravitational accelerations using Eq. (8) will be described. A. Space partitioning Before interpolation, the desired space should be selected. In this study, we use space partitioning to avoid interpolation in the whole space near the target asteroid while improving the interpolation performance. The space partitioning schemes based on the practical types of close proximity operations are used to select the target

9

interpolation space. The space which is not very close to the surface is selected using a minimum ellipsoid while the space containing the surface points are selected by cones. The detailed description is given below. (1) Space outside of a minimum ellipsoid In this case, the desired space is the space outside of a pre-defined ellipsoid. Here, an ellipsoid is used rather than a sphere because the elongated shapes of many asteroids are better approximated by ellipsoids. The ellipsoid can be used as a safe ellipsoid to keep a spacecraft away from the surface of the target asteroid. Calculating gravitational accelerations quickly in this space is necessary for onboard trajectory computation when a spacecraft orbits around the target asteroid outside of the ellipsoid. An ellipsoid can be defined by the following equation:

x2



y2



z2

  a   b   c  2

ref

2

2

1

(23)

ref

ref

where a, b, c are the axes of ellipsoids, the subscript ref is the referenced ellipsoid, and ρ is a scaling parameter which controls the size of the ellipsoid. The referenced ellipsoid is chosen as the circumscribed ellipsoid which is approximately solved by the following steps. First, searching three sets of points,

 p11 , p12  ,  p21 , p22 , p23 , p24  ,

 p31 , p32 , p33 , p34  , with the vertices of the polyhedrons. In the first set, p11 has the maximum x and p12 has the minimum x. In the second set, p21 has the maximum y and p22 has the minimum y among the points with positive x; p23 has the maximum y and p24 has the minimum y among the points with negative x. In the third set, p31 has the maximum z and p32 has the minimum z among the points with positive x; p33 has the maximum z and p34 has the minimum z among the points with negative x. Second, selecting an arbitrary point in each set and solving the three axes of ellipsoids with the three selected points. Third, the maximum a, b and c are chosen as aref, bref and cref. As shown in Fig 3, the approximately solved circumscribed ellipsoid is the inner ellipsoid with ρ =1. For safety, the selected ellipsoid can be a little larger than the circumscribed ellipsoid. For instance, the outer ellipsoid with ρ =1.1 as shown in this figure. It should be also pointed out that the gravitational accelerations inside of the circumscribed ellipsoid can be also approximated solved by Chebyshev polynomials if ρ is selected to be less than 1. The inner bound has been defined by the minimum ellipsoid. The outer bound of the interpolation space can also be described by defining an ellipsoid with an appropriate ρ. Because the interpolation space is large, we may further divide this space into several smaller sub spaces for achieving better approximating performance. Specifically, three space partitioning schemes between two defined ellipsoids will be used. These three schemes are denoted as EScheme1, 10

EScheme2 and EScheme3, respectively. EScheme1 do not divide the target space. EScheme2 divides the target space into two equal sub spaces. EScheme3 divides the target space into four equal sub spaces. Besides, a scheme without inner bound ellipsoid will also be used in comparison to show the necessary of setting the inner ellipsoid. This scheme is denoted as EScheme0.

Fig 3. Ellipsoids for selecting target interpolation space (2) Space in a cone In this case, the desired interpolation space is a cone above a surface point. A cone space is considered because it contains a surface point and satisfies descent glide slope constraints (Pinson and Lu, 2018). Calculating gravitational accelerations quickly in this kind of region can be used for onboard trajectory computation when a spacecraft needs to descent and land on the surface of the target asteroid. Besides, other implementations such as observations (Yang et al. 2018) and hopping (Liu et al. 2017) may also require calculating the gravitational accelerations in such conic regions. As shown in Fig 4, six cones on the surfaces of the asteroid are presented. As for the origins of these cones, two of them locate at the ends of the asteroid and four of them locate at the neck of the asteroid. The interpolation scheme using cones is denoted as CScheme.

11

Fig 4. Cones for selecting target interpolation space B. Variables for Interpolation As for the two types of interpolation space aforementioned, corresponding coordinate systems and variables are needed to be defined, respectively. The descriptions are given below. (1) Points on ellipsoids For this case, a coordinate system   ,  ,   is established as shown in Fig 5, where ρ reflects the magnitude of the position vector r, γ is the angle between the position vector, and the plane yz and θ is the angle between the projection of position vector in the plane yz and axis y. The range of γ is from -π/2 to π/2 and the range of θ is from –π to π. To describe a point on an ellipsoid described by Eq. (23), the position of this point can be expressed by x   aref sin 

(24)

y   bref cos  cos 

(25)

z   cref cos  sin 

(26)

With these defined coordinates, the description of the aforementioned space partitioning schemes is given in Table 1. According to Eqs. (24) - (26), θ, γ and ρ can also be solved analytically using the position vector [x, y, z]T, as follows:

z c

  a tan 2  ,

12

y  b

(27)

 y  x a tan 2  a , b cos  ，cos   0      z  a tan 2  x , ，cos   0   a c   sin 

(28)

 x  a sin   1, sin   0   y    1, sin   0 & cos   0  b cos   z  1, sin   0 & cos   0   c cos 

(29)

z

r

ρ

γ o

θ

y

x

Fig 5. Variables for Interpolation (Case 1)

Table 1 Schemes for space partitioning Subspace

EScheme1

EScheme2

EScheme3

A

-π ≤ θ ≤ π, -π/2 ≤ γ ≤ π/2

-π ≤ θ ≤ π, -π/2 ≤ γ ≤ 0

-π ≤ θ ≤ 0, -π/2 ≤ γ ≤ 0

B

-

-π ≤ θ ≤ π, 0 ≤ γ ≤ π/2

-π ≤ θ ≤ 0, 0 ≤ γ ≤ π/2

C

-

-

0 ≤ θ ≤ π, -π/2 ≤ γ ≤ 0

D

-

-

0 ≤ θ ≤ π, 0 ≤ γ ≤ π/2

In order to apply the Chebyshev polynomials for interpolation, non-dimensional variables should be used. The angles γ and θ can be transformed to the non-dimensional variables ξ and η according to Eq. (7) as

   min  ( max   min )

13

 1 2

(30)

   min  ( max   min )

 1 2

(31)

In the process of sampling, the CGL nodes are chosen as the discrete nodes for these two non-dimensional variables, which are (Junkins et al., 2012)

  j  i  cos  1  , j1  0,1, 2, M    M    j2   i  cos  M  , j2  0,1, 2, M    

(32)

According to the technique of radial adaptation (Junkins et al., 2012), the transforming equation connecting ρ and the non-dimensional  is chosen as

  1      4  

   max    max   min  cos 

(33)

where  is also sampled at the CGL nodes, which is

 j3 M  

 i  cos 

  , j3  0,1, 2, M  

(34)

In this way, the sampled nodes near the surfaces of the asteroids are denser than the sampled nodes near the ellipsoid with ρ = ρmax. Moreover, according to the work of Hu and Ji (2017), expressing the gravitational accelerations along the spherical coordinates instead of the rectangular coordinates can obtain better approximating accuracy. Therefore, we transform the coordinates of g from the coordinate system

 x, y, z

to the coordinate system   ,  ,   . The transformation

equation is given below:

 g    sin      g    0  g    cos   

cos  cos   sin  sin  cos 

cos  sin    g x  cos    g y  sin  sin    g z 

(35)

where the transformation matrix above is derived from the rotation matrix which connects the spherical coordinate system and the rectangular coordinate system in the book written by Battin (1999). (2) Points in cones In this case, the desired interpolation space is a cone with the cone vertex locating at a surface point. The variables selected for interpolation is r, δ and θ, as shown in Fig 6. The position vector can be expressed by with the aid of a

14

local coordinate system i , j , k , where k is a unit vector along the axis of the cone, i and j are orthogonal unit vectors which are perpendicular to k. In this study, we define i, j and k using the polyhedron data. For example, the cone vertex locates in a triangular with vertices numbered 1, 2 and 3. The position vectors of these vertices are r1, r2 and r3. Then, the unit vectors i, j and k are computed by

k

 r2  r1    r3  r2   r2  r1    r3  r2 

,i 

r2  r1 , j  ki r2  r1

(36)

Denote the position of the cone vertex as rc, the position vector of a sample point can be computed by

r  dr  cos  k  sin  sin  i  sin  cos  j   rc

(37)

In the numerical simulation of this study, the cone vertex will be chosen as the center of the surface triangular. Thus,

rc   r1  r2  r3  / 3 . The equations for transformation and generating discrete nodes in case have the same forms as Eqs. (30) - (34). For simplicity, these equations are omitted. Besides, the components of g are chosen as gx, gy and gz for this case unlike the previous case. It is because good fitting accuracy is easy to achieve within a small cone. Thus, related manipulations for transforming the components of g can be avoided. In addition, according to Eq. (37), dr, θ and δ can be computed using the following equations:

dr  r  rc

(38)

  a tan 2   r  rc   i ,  r  rc   j 

(39)

   r  rc   i  ,  r  rc   k  , sin   0 a tan 2  sin         r  rc   j   a tan 2  cos  ,  r  rc   k  , sin   0   

(40)

15

k

i

z

θ j

r

δ

o

y

x

Fig 6. Variables for Interpolation (Case 2) C. Algorithm of calculating gravitational accelerations With the defined variables and sample schemes, gravitational accelerations can be calculated using Chebyshev polynomials. According to the previous section, the sample nodes for both cases are  ,  and  while the selected T

T

components of g are  g  , g , g  and  g x , g y , g z  for the two cases, respectively. For simplicity, we denote both T

T

 g  , g , g  and  g x , g y , g z  as  g1 , g 2 , g3  . Replacing f in Eq. (8) by g1, g2 or g3, we can obtain a set of T

approximating equations for calculating g using Chebyshev polynomials, which are N N N  1  g1  , ,     a pqr Tp  Tq   Tr   0 0 0 p q r     N N N   2  g 2  , ,     a pqr Tp  Tq   Tr   0 0 0 p q r     N N N   g3  , ,     a 3pqr Tp  Tq   Tr    p 0 q 0 r 0

(41)

where a1pqr , a 2pqr and a 3pqr are the corresponding coefficients for g1, g2 and g3. In summary, the formulated algorithm has the following steps: Step 1: Initialization. Set the values of Mξ, Mη, Mζ, Nξ, Nη and Nζ. As for Case 1, solve aref, bref and cref and set the values of ρmax, ρmin, γmax, γmin, θmax and θmin; As for Case 2, solve the direction of the local coordinate system and set the values of rc, drmax, drmin, δmax, δmin, θmax and θmin.

16

Step 2: For each component of the gravitational acceleration, set f  , ,    gi  i  1, 2,3 and calculate the corresponding coefficients a ipqr by Eqs. (8) - (22). Step 3: Calculate the gravitational accelerations explicitly at the desired points by Eq.(41).

V.

Analysis and Validation of the Chebyshev Polynomial Method

In this section, the proposed Chebyshev polynomial method will be applied to a highly irregularly-shaped asteroid 1996 HW1 for analysis and validation. Besides, the effectiveness of this method will be discussed by comparing with other space partitioning schemes and three mascon models. The shape model of 1996 HW1 with 1392 vertices and 2780 facets (Magri et al. 2017) are used to generate the polyhedral gravitational field. The bulk density of 1996 HW1 is set to 3.56 g/cm3 (Wang et al., 2014). In the all following examples, two relative errors of the gravitational acceleration are defined for evaluating the accuracy of the used method compared with the polyhedral method. The first one is the percent error which is defined by

E

gC  g P gP

 100

(42)

where the subscripts C and P represent the Chebyshev polynomial method and the polyhedral method, respectively. The second one is the logarithmic relative error which is defined by

 g  gP E '  log10  C  gP 

  

(43)

The logarithmic relative error is used for distinguishing very small relative errors. Besides, the values of Mξ, Mη and Mζ are set to be equally, i.e. Mξ = Mη = Mζ = M. Meanwhile, the values of Nξ, Nη and Nζ are set to be the same, i.e. Nξ = Nη = Nζ = N. In simulations, the code for both the polyhedral method and the Chebyshev polynomial method is written in FORTRAN language and run in Release mode on a personal laptop with 2.8GHz CPU. A. Analysis and validation of the proposed space partitioning schemes First, the Chebyshev polynomial method is applied for calculating gravitational accelerations in the region outside of a minimum ellipsoid with the ellipsoid based space partitioning schemes. Using the method given in Section IV.

17

A, the parameters of the reference ellipsoid obtained are 1.9670, 1.1840 and 1.0503 km, respectively. ρmin is set to be 1.1 for schemes EScheme1, EScheme2 and EScheme3, which means the inner boundaries of the target regions are close to the corresponding reference ellipsoids. As for ρmax, it should be set to an appropriate value such that desired space is inside of an ellipsoid whose axes are ρmaxaref, ρmaxbref and ρmaxcref. Herein, we set ρmax to 4. For each scheme, M is set to 25 and N is varied from 5 to 15. The results of the percent error on the minimum ellipsoid and in the equator plane are shown in Fig 7 and Fig 8. The sampled points in the equator plane are not inside of the minimum ellipsoid and their positions satisfy -4 < x < 4 and -4 < y < 4. Comparing EScheme0 with EScheme1, EScheme2 and EScheme3, we can see that the approximating error of EScheme0 is much larger in each subfigure. Thus, using a minimum ellipsoid to reduce the interpolation space is useful and necessary to lower the approximating error. As for schemes using the minimum ellipsoid, according to the figures, the tendency of both the maximum percent error and the average percent error with respect to the increasing N is decreasing. Especially for the cases when N is small, the relative error decreases quickly with the increasing N. It also can be seen that the average percent error is much smaller than the maximum percent error. Comparing the results of EScheme1, EScheme2 and EScheme3, we can see that EScheme3 leads to the minimum approximating error. The reason why EScheme3 is the most accurate is that the sampled subspaces for interpolation are much smaller than those of the other two schemes. Thus, partitioning the space between two ellipsoids before interpolation is beneficial. To see the accuracy of these schemes more clearly, a zoom-in picture is presented in Fig 7 (a). From this figure, we can see that the smallest N to make the maximum percent error be less than 1% is 10 for EScheme3 while it is 12 for EScheme2. As for EScheme1, the maximum percent error is larger than 1% even when N is 15. We emphasize that EScheme3 can have much better performance compared with EScheme1 and EScheme2 for other asteroids, which will be presented in next section. The distributions of the approximating error on the minimum ellipsoid and in the equator plane using EScheme3 with N = 10 shown in Fig 9 (a) again validates the maximum percent error on each ellipsoid is smaller than 1%. According to this figure, the percent error is close to 1% only in several small areas. It is far smaller than 1% in the most areas and even smaller than 0.1% in some areas. From Fig 9 (b), it can be seen that the approximating error is near 1% only for the points close to the surface of the minimum ellipsoids and it decreases quickly when the distance to the ellipsoid surface increases. Then, we test the computational efficiency of the proposed method. The selected test points are 101 × 101 points which uniformly distributed on the minimum ellipsoids. The reported CPU time is only 3.13 × 10-2 s for computing

18

the gravitational acceleration using EScheme3 with N = 10 while it is 2.41 s when the polyhedral method is used. Therefore, the proposed method is very fast for calculating the gravitational field of the irregularly-shaped asteroid. Meanwhile, the total required storage for the interpolation coefficients is only 372 kb. 50

25

4

EScheme0 EScheme1 EScheme2 EScheme3 1%

3.5 3

40 E (%)

2.5 2

20

1

E (%)

E (%)

1.5

30

0.5 0

20

9

10

11

12

13

14

15

15

10

N

10

0

EScheme0 EScheme1 EScheme2 EScheme3 1%

5

5

10

0

15

5

10

N

15

N

(a) maximum

(b) average

Fig 7. Percent error on the ellipsoid with ρ = 1.1 50

140 EScheme0 EScheme1 EScheme2 EScheme3 1%

120

40

80

E (%)

E (%)

100

EScheme0 EScheme1 EScheme2 EScheme3 1%

60

30

20

40 10

20 0

5

10

0

15

N

5

10

N

(a) maximum

(b) average

Fig 8. Percent error in the equator plane outside of the ellipsoid with ρ = 1.1

19

15

(a) minimum ellipsoid

(b) equator plane

Fig 9. Distribution of percent error with ρ = 1.1 Next, we analysis EScheme3’s ability to approximate the gravitational field near the surface of the asteroid. ρmax is set to 1.1 and ρmin is set to 0.4. We still choose N as 10. The distribution of the approximating error in the equator plane is shown in Fig 10 (a). From this figure, one can see the percent error reaches up to about 10% near the asteroid surface. The, we increase N to 15. The corresponding result is shown in Fig 10 (b). The maximum approximating error is decreased to about 8%. Although the maximum approximating error is not small, the approximating error is smaller than 3% at most of the test points.

(a) N = 10

(b) N = 15

Fig 10. Distribution of percent error with ρmax = 1.1 and ρmin = 0.4 Last, we test the performance of CScheme to show the effectiveness of the proposed method for calculating gravitational accelerations when a spacecraft descends in a cone. In simulations, the vertices of the cones are chosen as the points at the two ends and the neck for each elongated asteroid. M is fixed to 25 and N is fixed to 15. δmin, δmax, θmin, θmax and drmin are set to 0, π/6, 0, 2π and 0. drmax is set to 1.5 km. The distribution of the logarithmic relative error 20

for the selected conic regions is shown in Fig 11. According to this figure, the maximum relative error is less than 1% which is much smaller than that in Fig 10. The averaged error is only 0.0094 %. Thus, CScheme is more suitable for computing highly accurate gravitational acceleration than EScheme3 in the descent phase.

Fig 11. Distribution of the logarithmic relative error in the selected conic regions B. Comparison with other space partitioning schemes In this part, the proposed ellipsoid based space partitioning schemes are compared to the sphere based space partitioning schemes. A sphere based scheme has been used in study of Junkins et al. (2012), in which a middle sphere with radius of 1.05 Re is used to partition the space between two spheres with radii of Re and 7 Re, where Re represents the radius of Earth. In the current study, the sphere based space partitioning schemes are special cases of the ellipsoid based space partitioning schemes. The schemes corresponding to EScheme0, EScheme1, EScheme2 and EScheme3 are denoted as SScheme0, SScheme1, SScheme2 and SScheme3, respectively. When applying these schemes, we need to set a = b = c. Then, the equation derived in Section IV can be used. In simulations, c and ρmax are set to 1.0503 km and 4; ρmin is set to 1.1 for the schemes using a minimum sphere. The obtained results of the approximating error on the minimum ellipsoid and in the equator are shown in Fig 12 and Fig 13, respectively. Compared the results of SScheme0 with SScheme1, SScheme2 and SScheme3, it can be found that setting an inner sphere above the surface of the asteroid to partition the interpolation space can also greatly decrease the approximating error like the ellipsoid based schemes. According to results of SScheme1, SScheme2 and SScheme3 in Fig 12 and Fig 13, both the maximum and average approximating errors are much larger than 1% for

21

each scheme. Comparing these results with those in Fig 7 and Fig 8, we can see the ellipsoid based space partitioning schemes lead to much better performance in approximating accuracy than the sphere based schemes. 30

20 SScheme0 SScheme1 SScheme2 SScheme3 1%

25

15

E (%)

E (%)

20

SScheme0 SScheme1 SScheme2 SScheme3 1%

15

10

10 5 5 0

5

10

0

15

5

N

10

15

N

(a) maximum

(b) average

Fig 12. Percent error on the ellipsoid with ρ = 1.1 using sphere based schemes 30

350 SScheme0 SScheme1 SScheme2 SScheme3 1%

300

25 20

200

E (%)

E (%)

250

150

15 10

100

5

50 0

SScheme0 SScheme1 SScheme2 SScheme3 1%

5

10

0

15

N

5

10

15

N

(a) maximum

(b) average

Fig 13. Percent error in the equator plane using sphere based schemes C. Comparison with mascon models In this part, the proposed method will be compared with mascon model methods. First, a recently proposed simplified model method (Burov et al. 2019) which is indeed a simplified mascon method using four mass points is chosen. Before comparison, we need to mention that the simplified mascon method is quite simple and requires extremely low computational effort because this method only needs to compute the gravitational fields of four mass points. Here, we only compare the proposed method with the four-mass-point model to show the advantage in approximating accuracy.

22

According to the method in Burov et al. (2019), the positions of the four mass points obtained from the polyhedral model are [1.021033, 0.225192, -0.011823], [-1.188637, 0.036059, -0.003048], [0.043009, -0.033385, 0.014009], [0.716618, -0.378993, 0.008434] km, respectively. The corresponding gravitational constants are 281.462, 359.762, 205.420, 183.374 m3/s2, respectively. The results of the approximating error on the minimum ellipsoid and in the equator using the four-mass-point model are shown in Fig 14. The red dash line represents the ellipsoid with ρ = 1.1. According to Fig 14 (b), the accuracy of this model can be less than 2% even 1% when the test point is not close to the asteroid surface. However, the approximating error increases quickly as the test point move close to the asteroid. According to Fig 14 (a), the maximum error on the minimum ellipsoid is about 16%. From Fig 14 (b), one can see the error outside of the minimum ellipsoid can be larger than 5% even 10%. Compared Fig 14 with Fig 9, the proposed method with EScheme3 is able to achieve much higher accuracy than the four-mass-point model.

180

16

150

14

120 90

12

(deg)

60

10

30 0

8

-30

6

-60 -90

4

-120

2

-150 -180 -90

-60

-30

0

30

60

90

0

(deg)

(a) ellipsoid with ρ = 1.1

(b) equator plane

Fig 14. Distribution of percent error using the four-mass-point model Next, the mascon model methods in Chanut et al. (2015) are chosen for comparison. Their study provides an effective approach to obtain the data from the polyhedral model. When the Mascon 1 model is used, the number of the mass points is equal to the tetrahedrons of the polyhedral model. Specifically, the Mascon 1 contains 2780 mass points for the HW1 case. The results of the approximating error on the minimum ellipsoid and in the equator using the Mascon 1 model are shown in Fig 15. According to this figure, the approximating error using this model is still much larger than that using EScheme3. Then, the Mascon 3 model which has better accuracy than the Mascon 1 model is also tested. In simulation, the points with the coordinates 2/3rvi and 4/9 rvi, where rvi is the position of the ith vertex, are used to divide each tetrahedron into three parts. 8340 mass points are used in the Mascon 3 model of the HW1 case. The corresponding results of the approximating error are shown in Fig 16. According to this figure, the accuracy 23

of the Mascon 3 model is much higher than the Mascon 1 model. Nevertheless, according to the 4% relative error on the ellipsoid, the accuracy of the Mascon 3 model is lower than that of the proposed method with EScheme3. Besides, the reported CPU time is 0.375 s when the Mascon 3 model is used for calculate the gravity at the 101 × 101 test points on the ellipsoid, which means the proposed method can achieve higher computational speed over the Mascon 3 model. 180 12

150 120

10

90

(deg)

60

8

30 0

6

-30 -60

4

-90 -120

2

-150 -180 -90

-60

-30

0

30

60

0

90

(deg)

(a) ellipsoid with ρ = 1.1

(b) equator plane

Fig 15. Distribution of percent error using the Mascon1 model

180 150

4

120

3.5

90

3

(deg)

60 30

2.5

0

2

-30 -60

1.5

-90

1

-120

0.5

-150 -180 -90

-60

-30

0

30

60

0

90

(deg)

(a) ellipsoid with ρ = 1.1

(b) equator plane

Fig 16. Distribution of percent error using the Mascon3 model

VI.

Numerical Examples and Results

In this section, the proposed method is further applied to other three asteroids for validating its effectiveness. The selected asteroids are 433 Eros, 25143 Itokawa and 101955 Bennu, which have different size and shapes compared with 1996 HW1. The related physical parameters for these asteroids are given in Table 2. The numbers of the vertices 24

and facets in the last column indicate the chosen shape models for modelling the asteroids. It can be seen that the model resolution of Bennu is relative low while that of Eros and Itokawa is much higher. Comparing the approximating performance for asteroids with different resolution can reflect the influence of the model resolution. With the shape data, the asteroids are modelled by polyhedrons, as shown in Fig 17. From the figures, it can be seen that the asteroids have different irregular shapes. Table 2 Physical parameters for the selected Asteroids Asteroid

Density (g/cm3)a

Dimension (km)

Vertices and facets b-d

433 Eros

2.67

32.78 ×14.58 ×11.96

25350 and 49152

25143 Itokawa

1.95

0.56 × 0.30 × 0.24

25350 and 49152

101955 Bennu

0.97

0.57 × 0.53 × 0.50

1348 and 2692

a: Wang et al., 2014; b: Gaskell, 2008; c: Gaskell et al., 2008; d: Nolan et al., 2019.

(a) Eros

(b) Itokawa

(c) Bennu

Fig 17. Modeling asteroids with polyhedrons A. Calculating Gravitational accelerations outside of a minimum ellipsoid In this part, the ellipsoid based space partitioning schemes are tested for the asteroid Eros, Itokawa and Bennu. The parameters of the reference ellipsoids for the four asteroids are given in Table 3. In simulations, ρmin is set to be 1.1 and ρmax is set to 5, 7 and 4 for these three asteroids, respectively.

Table 3 Parameters of the reference ellipsoids Asteroid

aref

bref

cref

433 Eros

20.4887

10.2009

7.4300

25

25143 Itokawa

0.3590

0.1785

0.2964

101955 Bennu

0.2901

0.2730

0.2643

The results of the approximating error with EScheme1, EScheme2 and EScheme3 are shown in Fig 18 - Fig 20. According to the curves in each figure, both the maximum percent error and the average percent error decrease when N increases. EScheme3 is still the best scheme as the HW1 case. From the figures, one can find the performance of EScheme1 and EScheme2 are close for the Itokawa and Bennu cases. However, EScheme3 have much better performance in accuracy compared with EScheme1 and EScheme2 for each case, especially for the Itokawa and Bennu cases. Besides, EScheme3 is able to make the maximum percent error to be less than 1% for each asteroid. The smallest N to make the maximum percent error be less than 1% is 12, 9 and 12 for the four asteroids, respectively. When N takes these values with EScheme3, the corresponding average percent error is far smaller than 1% which indicates the approximating accuracy is very high. 30

10 EScheme 1 EScheme 2 EScheme 3 1%

25

EScheme 1 EScheme 2 EScheme 3 1%

8

E (%)

E (%)

20 15

6

4 10 2

5 0

5

10

0

15

5

10

N

N

(a) maximum

(b) average

Fig 18. Percent error on the ellipsoid with ρ = 1.1 (Eros)

26

15

12

35 EScheme 1 EScheme 2 EScheme 3 1%

30

8

20

E (%)

E (%)

25

15

6 4

10

2

5 0

EScheme 1 EScheme 2 EScheme 3 1%

10

5

10

0

15

5

N

10

15

N

(a) maximum

(b) average

Fig 19. Percent error on the ellipsoid with ρ = 1.1 (Itokawa)

10

2 EScheme 1 EScheme 2 EScheme 3 1%

1.5

6

E (%)

E (%)

8

EScheme 1 EScheme 2 EScheme 3 1%

1

4 0.5

2

0

5

10

0

15

N

5

10

15

N

(a) maximum

(b) average

Fig 20. Percent error on the ellipsoid with ρ = 1.1 (Bennu)

Then, N is set to 12, 9 and 12 for the three asteroids, respectively. The corresponding distributions of the approximating error with EScheme3 are depicted in Fig 21 and Fig 22. From Fig 21, the maximum percent error for each case is smaller than 1% and the approximating error is much smaller than 1% for the most of the area. From Fig 22 (a) and (c), the approximating error is less than 1% at each test point and it becomes far smaller than 1% when a test point is far away from the surface. Besides, according to Fig 22 (b), the distribution of the percent error in the equator plane of Itokawa is different compared with the other three. It can be seen that the percent error becomes larger than 1% at the points far away from the asteroid. This phenomenon is due to the selected ρmax for Itokawa is larger than other three asteroids and also the selected N is smaller such that the sample points far away from Itokawa

27

are relatively sparse. In order to further lower the percent error, we can simply increase N to 10. The corresponding results are shown in Fig 23. From this figure, it can be seen that the maximum percent error in the equator plane becomes smaller than 1% and the percent errors on the minimum ellipsoid also decreases. Now, the accuracy of the proposed method with EScheme3 in desired regions has again been validated for three different shaped asteroids. The proposed method is also very efficient for calculating gravitational accelerations. The test points are selected as the HW1 case. The comparison of the computational speed between the proposed method and the polyhedral method is shown in Table 4. According to this table, the speedup of the Chebyshev polynomial method over the polyhedral method ranges from 44 times to 1563 times. Besides, the total required storages for the coefficients are only 372 kb and 420 kb with EScheme3 when N is 10 and 12, respectively. Due to the high computational speed and low required storage, the proposed method is promising for onboard computation of asteroid gravitational acceleration.

(a) Eros

(b) Itokawa

(c) Bennu

Fig 21. Distribution Percent error on the ellipsoid with ρ = 1.1

(a) Eros

(b) Itokawa

(c) Bennu

Fig 22. Percent error in the equator plane outside of the ellipsoid with ρ = 1.1

28

(a) ellipsoid

(b) equator plane

Fig 23. Percent error for Itokawa with ρ = 1.1

Table 4 Required CPU time for calculating gravitational accelerations Asteroid

N

Chebyshev polynomial method (s)

polyhedral method (s)

433 Eros

12

6.25 × 10-2

52.17

25143 Itokawa

10

3.13 × 10-2

48.92

101955 Bennu

12

6.25 × 10-2

2.75

In addition, it should be noted that the bulk densities of asteroids are all assumed to be constant in numerical examples above. However, the density inhomogeneity of some natural asteroids, e.g. Itokawa, is obvious (Lowry et al. 2014). The Itokawa case is reconsidered to test the effectiveness of the method when the density is inhomogeneous. According to Lowry et al. (2014) and Scheeres and Gaskell (2008), the head density is 2.85 g/cm3 while the body density is 1.75 g/cm3. For simplicity, we set σ in Eq. (1) to 2.85 g/cm3 or 1.75 g/cm3 according to whether the x coordinate of the center point of the face or edge is greater than 150 m or not. The corresponding results are shown in Fig 24. According to this figure, the approximating error is still much less than 1 % on the ellipsoid and is a little larger than 1% only at some distant small regions, which means the proposed method is still effective for achieving high accuracy.

29

(deg)

180 150

0.5

120

0.45

90

0.4

60

0.35

30

0.3

0

0.25

-30

0.2

-60

0.15

-90 -120

0.1

-150

0.05

-180 -90

-60

-30

0

30

60

90

0

(deg)

(a) ellipsoid

(b) equator plane

Fig 24. Percent error for Itokawa with Inhomogeneous density B. Calculating Gravitational accelerations in a cone In this part, we test the effectiveness of CScheme for the three asteroids. δmin, δmax, θmin, θmax and drmin are set to 0, π/6, 0, 2π and 0. drmax is set to 10 km, 0.3 km and 0.2 km for the three asteroids, respectively. The distributions of the relative errors in the conic regions are shown in Fig 25. According to these figures, the maximum relative errors appear at the points in the cones which locate at the necks of the elongated asteroids. Importantly, the relative errors are extremely small. Even for the maximum relative error, it can be much less than 1%. As for Bennu, the maximum relative error for each cones are all less than 10-4. Table 5 gives the details about the maximum, minimum and averaged relative errors. According to the approximating results shown in Fig 25, the proposed method is suitable for calculating gravitational accelerations when a spacecraft descends in a cone. In addition, the required storages for one cone with N = 15 is only 195 kb. And the speedup of the Chebyshev polynomial method over the polyhedral method is about 20 times for the asteroids with low resolution shape models and over 300 times for the asteroids with high resolution shape models according to the reported CPU time. Thus far, the effectiveness of the proposed method with CScheme is again validated for the asteroids Eros, Itokawa and Bennu.

Table 5 Percent errors in the conic regions with M = 25 and N = 15 Asteroid

Maximum error (%)

Minimum error (%)

Average error (%)

433 Eros

0.1080

1.51 ×10-5

0.0027

25143 Itokawa

0.0874

7.42

×10-5

0.0034

101955 Bennu

0.0024

5.07 ×10-6

30

2.83 ×10-4

(a) Eros

(b) Itokawa

(c) Bennu

Fig 25. Logarithmic relative error in the selected conic regions

VII.

Conclusion

This paper presents an improved Chebyshev polynomial method for solving gravitational accelerations near irregularly-shaped asteroids efficiently. Space partitioning schemes using Chebyshev-Gauss-Lobatto nodes and radial adaption are derived for interpolation outside a minimum ellipsoid and inside a cone. The coefficients of the Chebyshev polynomials are solved in a reduced storage approach using the analytical least square solution. Through simulation, the performance of different space partitioning schemes is compared. It is found that the ellipsoid based schemes are more accurate than the sphere based schemes. The partitioning the space between two ellipsoids is necessary to obtain high approximating accuracy. Besides, the interpolation scheme using a cone is better than the ellipsoid based schemes due to its higher approximating accuracy near the asteroid surface. Simulation results indicate the method is applicable to asteroids 1996 HW1, 433 Eros, 25143 Itokawa and 101955 Bennu which have different shape, size and model resolution. The maximum relative error between Chebyshev polynomial method and the polyhedral method in the test space can be less than 1% and the average relative error is far smaller than 1%. The speedup of the proposed method over the polyhedral method is up to larger than a thousand times according to the record CPU time. Meanwhile, the required storage for the coefficients in the test cases is about three or four hundred kilobytes for interpolation outside of a minimum ellipsoid and less than 200 kilobytes for interpolation in a cone. The results indicate the proposed method is promising for onboard implementation in asteroid missions.

31

Acknowledgments This work was supported by the Natural Science Foundation of Jiangsu Province (BK20180410), the Young Elite Scientists Sponsorship Program by CAST (2018QNRC001), Project funded by China Postdoctoral Science Foundation (2018M640482) and Jiangsu Postdoctoral Research Funding Program (2018K115C).

References Battin, R. H. 1999. An introduction to the mathematics and methods of astrodynamics. AIAA Education Series, AIAA, New York, pp. 82. Bazzocchi, M.C. and Emami, M.R., 2018. Formation of multiple landers for asteroid detumbling. Advances in Space Research, 62(3), 732-744. Broschart, S. B., and Scheeres, D. J. 2005. Control of hovering spacecraft near small bodies: application to asteroid 25143 Itokawa. Journal of Guidance, Control, and Dyn. 28(2), 343-354. Burov, A. A., Guerman, A. D., Nikonova, E. A., and Nikonov, V. I. 2019. Approximation for attraction field of irregular celestial bodies using four massive points. Acta Astronaut. 157, 225-232. Chanut, T. G. G., Aljbaae, S., and Carruba, V. 2015. Mascon gravitation model using a shaped polyhedral source. Monthly Notices of the Royal Astronomical Society 450(4), 3742-3749. Chanut, T. G. G., Winter, O. C., Amarante, A., and Araújo, N. C. S. 2015. 3D plausible orbital stability close to asteroid (216) Kleopatra. Monthly Notices of the Royal Astronomical Society 452(2), 1316-1327. D'Ambrosio, A., Circi, C. and Zeng, X. 2019. Solar-photon sail hovering orbits about single and binary asteroids. Advances in Space Research, 63(11), 3691-3705. Furfaro, R., Cersosimo, D., and Wibben, D. R. 2013. Asteroid precision landing via multiple sliding surfaces guidance techniques. Journal of Guidance, Control, and Dyn. 36(4), 1075-1092. Gao, A., and Liao, W. 2019. Efficient gravity field modeling method for small bodies based on Gaussian process regression. Acta Astronaut. 157, 73-91. Gaskell, R.W. 2008. Gaskell Eros Shape Model V1.0. NEAR-A-MSI-5-EROSSHAPE-V1.0. NASA Planetary Data System. https://sbn.psi.edu/pds/resource/erosshape.html [retrieved on 20 April, 2019]. Gaskell, R., Saito, J., Ishiguro, M., Kubota, T., Hashimoto, T., Hirata, N., Abe, S., Barnouin-Jha, O., and Scheeres, D., Gaskell Itokawa Shape Model V1.0. 2008. HAY-A-AMICA-5-ITOKAWASHAPE-V1.0. NASA Planetary Data System. https://sbn.psi.edu/pds/resource/itokawashape.html [retrieved on 20 April, 2019]. Huang, X., Cui, H., and Cui, P. 2004. An autonomous optical navigation and guidance for soft landing on asteroids. Acta Astronaut. 54(10), 763-771.

32

Hu, S. C., and Ji, J. H. 2017. Using Chebyshev polynomial interpolation to improve the computational efficiency of gravity models near an irregularly-shaped asteroid. Research in Astronomy and Astrophysics 17(12), 120. Hamilton, V. E., Simon, A. A., Christensen, P. R. et al., 2019. Evidence for widespread hydrated minerals on asteroid (101955) Bennu. Nature Astronomy 3(4), 332. Jiang, Y., Baoyin, H. and Li, H. 2018. Orbital stability close to asteroid 624 Hektor using the polyhedral model. Advances in Space Research, 61(5), 1371-1385. Junkins, J. L., Younes, A. B. and Bai, X. 2012. “Orthogonal Polynomial Approximation in Higher Dimensions: Applications in Astrodynamics”, In: AAS Jer-Nan Juang Astrodynamics Symposium, College Station, TX. Kitazato, K., Milliken, R. E., Iwata, T. et al., 2019. The surface composition of asteroid 162173 Ryugu from Hayabusa2 nearinfrared spectroscopy. Science 364 (6437), 272-275. Lauretta, D. S., DellaGiustina, D. N., Bennett, C. A. et al., 2019. The unexpected surface of asteroid (101955) Bennu. Nature. 568, 55-60. Lauretta, D. S. and OSIRIS-REx Team, 2012. “An overview of the OSIRIS-REx asteroid sample return mission”. In: 43rd Lunar and Planetary Institute Science Conference, Woodlands, TX. Li, X., Gao, A., and Qiao, D. 2017. Periodic orbits, manifolds and heteroclinic connections in the gravity field of a rotating homogeneous dumbbell-shaped body. Astrophysics and Space Science 362(4), 85. Li, X., Qiao, D., and Cui, P. 2013. The equilibria and periodic orbits around a dumbbell-shaped body. Astrophysics and Space Science 348(2), 417-426. Liu, X., Baoyin, H., and Ma, X. 2011. Equilibria, periodic orbits around equilibria, and heteroclinic connections in the gravity field of a rotating homogeneous cube. Astrophysics and Space Science 333(2), 409. Liu, Y., Zhu, S., Cui, P., Yu, Z., and Zong, H. 2017. Hopping trajectory optimization for surface exploration on small bodies. Advances in Space Research 60(1), 90-102. Lowry, S. C., Weissman, P. R., Duddy, S. R., et al. 2014. The internal structure of asteroid (25143) Itokawa as revealed by detection of YORP spin-up. Astronomy & Astrophysics, 562, A48. Magri, C., Howell, E.S., Nolan, M.C., Taylor, P.A., Fernandez, Y.R, Mueller, M., Vervack, R.J., Benner, L.A.M., Giorgini, J.D., Ostro, S.J., Scheeres, D.J., Hicks, M.D., Rhoades, H., Somers, J.M., Gaftonyuk, N.M., Kouprianov, V.V., Krugly, Yu.N., Molotov, I.E., Busch, M.W., Margot, J.L., Benishek, V., Protitch-Benishek, V., Galad, A., Higgins, D., Kusnirak, P., and Pray, D.P., 2017. Shape and Rotation of (8567) 1996 HW1 V1.0. EAR-A-I0037-5-SHAPE8567-V1.0. NASA Planetary Data System. https://sbn.psi.edu/pds/resource/shape8567.html [retrieved on 20 April, 2019].

33

Nolan, M.C., Magri, C., Howell, E.S., Benner, L.A.M., Giorgini, J.D., Hergenrother, C.W., Hudson, R.S., Lauretta, D.S., Margot, J.L., Ostro, S.J., and Scheeres, D.J. 2013. Asteroid (101955) Bennu Shape Model V1.0. EAR-A-I0037-5-BENNUSHAPEV1.0. NASA Planetary Data System. https://sbn.psi.edu/pds/resource/bennushape.html [retrieved on 20 April, 2019]. Pearl, J. M., Louisos, W. F., and Hitt, D. L. 2018. Hybrid Gravity Model from Asteroid Surface Topology. In: AAS/AIAA Space Flight Mechanics Meeting, Kissimmee, Florida. Pinson, R. and Lu, P. 2018. Trajectory design employing convex optimization for landing on irregularly shaped asteroids, Journal of Guidance, Control, and Dyn. 41(6), 1243-1256. Riaguas, A., Elipe, A., and López-Moratalla, T. 2001. Non-linear stability of the equilibria in the gravity field of a finite straight segment. Celest. Mech. Dyn. Astron. 81(3), 235-248. Scheeres, D. J. 1994. Dynamics about uniformly rotating triaxial ellipsoids: applications to asteroids. Icarus 110(2), 225-238. Scheeres, D. J., and Gaskell, R. W. 2008. Effect of density inhomogeneity on YORP: The case of Itokawa. Icarus, 198(1), 125129. Scheeres, D. J., Williams, B. G., and Miller, J. K. 2000. Evaluation of the dynamic environment of an asteroid: Applications to 433 Eros. Journal of Guidance, Control, and Dyn. 23(3), 466-475. Smith, R. L. and Lyubomirsky A. S. 1981. “Techniques for increasing the efficiency of Earth gravity calculations for precision orbit determination”. In: NASA. Goddard Space Flight Center Sixth Ann. Flight Mech./Estimation Theory Symp, Silver Spring, MD. Surovik, D. A., and Scheeres, D. J. 2014. “Autonomous maneuver planning at small bodies via mission objective reachability analysis”. In: AIAA/AAS Astrodynamics Specialist Conference, San Diego California. Surovik, D. A., and Scheeres, D. J. 2015. Adaptive reachability analysis to achieve mission objectives in strongly non-keplerian systems. Journal of Guidance, Control, and Dyn. 38(3), 468-477. Takahashi, Y., and Scheeres, D. J. 2014. Small body surface gravity fields via spherical harmonic expansions. Celest. Mech. Dyn. Astron. 119(2), 169-206. Takahashi, Y., Scheeres, D. J., & Werner, R. A. 2013. Surface gravity fields for asteroids and comets. Journal of Guidance, Control, and Dyn. 36(2), 362-374. Tsuda, Y., Yoshikawa, M., Abe, M., Minamino, H., and Nakazawa, S., 2013. System design of the Hayabusa 2—Asteroid sample return mission to 1999 JU3. Acta Astronaut. 91, 356-362. Wang, X., Jiang, Y., and Gong, S. 2014. Analysis of the potential field and equilibrium points of irregular-shaped minor celestial bodies. Astrophysics and Space Science 353(1), 105-121. Werner, R. A. 1994. The gravitational potential of a homogeneous polyhedron or don't cut corners. Celest. Mech. Dyn. Astron. 59(3), 253-278.

34

Werner, R. A., and Scheeres, D. J. 1996. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia. Celest. Mech. Dyn. Astron. 65(3), 313-344. Yang, H., Bai, X., and Baoyin, H. 2017. Rapid generation of time-optimal trajectories for asteroid landing via convex optimization. Journal of Guidance, Control, and Dyn. 40(3), 628-641. Yang, H., Bai, X., and Li, S. 2018. Artificial Equilibrium Points near Irregular-Shaped Asteroids with Continuous Thrust. Journal of Guidance, Control, and Dyn. 41(6), 1308-1319. Yang, H., Li, S., and Bai, X. 2019. Fast Homotopy Method for Asteroid Landing Trajectory Optimization Using Approximate Initial Costates. Journal of Guidance, Control, and Dyn. 42(3), 585-597. Yu, Y., and Baoyin, H. 2012a. Generating families of 3D periodic orbits about asteroids. Monthly Notices of the Royal Astronomical Society 427(1), 872-881. Yu, Y., and Baoyin, H. 2012b. Orbital dynamics in the vicinity of asteroid 216 Kleopatra. The Astronomical Journal 143(3), 62. Yu, Y., and Baoyin, H. 2013. Resonant orbits in the vicinity of asteroid 216 Kleopatra. Astrophysics and Space Science 343(1), 75-82. Zeng, X., Gong, S., Li, J., and Alfriend, K. T. 2016. Solar sail body-fixed hovering over elongated asteroids. Journal of Guidance, Control, and Dyn. 39 (6), 1223-1231. Zeng, X., Jiang, F., Li, J., and Baoyin, H. 2015. Study on the connection between the rotating mass dipole and natural elongated bodies. Astrophysics and Space Science 356(1), 29-42. Zeng, X., Zhang, Y., Yu, Y., and Liu, X. 2018. The dipole segment model for axisymmetrical elongated asteroids. The Astronomical Journal 155(2), 85.

35

Declaration of Interest Statement

The Authors declared that they had no conflicts of interests in their authorship and publication of this Contribution

36