Feb 15, 2017 - Keplerian orbit is defined as an orbit in which a perturbative or ... orbit to sit above/below the standard 70 km station-keeping box o...

0 downloads 0 Views 1MB Size

http://eprints.gla.ac.uk/136992/

Deposited on: 15 February 2017

Enlighten – Research publications by members of the University of Glasgow http://eprints.gla.ac.uk33640

(Preprint) AAS 17-278

OSCULATING KEPLERIAN ELEMENTS FOR HIGHLY NONKEPLERIAN ORBITS Alessandro Peloni,* Colin R. McInnes,† and Matteo Ceriotti‡ This paper presents a mapping between the elements of highly non-Keplerian orbits and classical orbital elements. Three sets of elements are discussed and mappings are derived in closed, analytical form for both the direct and inverse problem. Advantages and drawbacks of the use of each set of elements are discussed. The spacecraft thrust-induced acceleration used to generate families of highly non-Keplerian orbits is extracted from the inverse mapping from the osculating orbital elements. The key signatures of highly non-Keplerian orbits in Keplerian elements tracking data are determined through a set of representative test cases. INTRODUCTION In classical Keplerian mechanics, an object of negligible mass orbits around an attractor due to the gravitational interaction existing between the two bodies. Celestial mechanics traditionally deals with Keplerian orbits, treating any non-Keplerian effect as a perturbation to the classical Keplerian orbit. However, a nonKeplerian orbit is defined as an orbit in which a perturbative or propulsive acceleration acts in addition to the gravitational attraction of the primary body. If the magnitude of such acceleration averaged over one orbit is comparable with the gravitational and centripetal accelerations experienced by the object, the orbit can be classified as highly non-Keplerian.1 For the sake of conciseness and because only Keplerian and highly non-Keplerian orbits are considered in this paper, all the highly non-Keplerian orbits will be simply referred to as non-Keplerian orbits (NKOs) in this paper. New families of NKOs can be generated using continuous low-thrust propulsion.1-4 These include orbits displaced above/below the orbit plane of a conventional Keplerian orbit using out-of-plane thrust, or inside/outside a Keplerian orbit using radial thrust. Potential applications of NKOs include orbits displaced above/below the geostationary ring to increase the number of available slots for communications platforms4 and on-orbit inspection by formation-flying above/below or inside/outside the orbit of a target spacecraft. 5, 6 For example, a displacement distance of 35 km north/south of the geostationary ring allows a platform on a displaced highly non-Keplerian orbit to sit above/below the standard 70 km station-keeping box of a conventional geostationary platform. Such an orbit would require a thrust of 200 mN for a 1000 kg spacecraft, which is achievable with a single QinetiQ T6+ thruster.4 To date, several investigations have been carried out to study non-Keplerian orbits.1-7 Nevertheless, all prior studies on non-Keplerian orbits have focused on orbit and mission design. A key open issue in the analysis of the families of non-Keplerian orbits is the generation of a mapping from the NKO geometry to a corresponding set of osculating orbital elements. This is of importance to support the future operational use Research Assistant, School of Engineering, University of Glasgow, Glasgow G12 8QQ – UK. Professor of Engineering Science, James Watt Chair, School of Engineering, University of Glasgow, Glasgow G12 8QQ – UK. ‡ Lecturer in Space System Engineering, School of Engineering, University of Glasgow, Glasgow G12 8QQ – UK. *

††

1

of non-Keplerian orbits. Since families of NKOs are generated using a strong, continuous perturbation, their osculating orbital elements will be time varying, although still periodic. Similarly, the inverse problem is also of importance to deduce the NKO geometry from observations of orbital elements. Once these elements are known, the thrust-induced acceleration magnitude and direction required for the NKO can be determined. In order to pursue the above questions, the key objectives of this paper are to: a) map the properties of families of highly non-Keplerian orbits onto the classical orbital elements; b) generate an inverse mapping from orbital elements to the properties of the non-Keplerian orbits; and c) determine the key signatures of non-Keplerian orbits in orbital element tracking data. These objectives represent a largely unexplored aspect of the dynamics of families of NKOs and offer a route towards their future operational use. The paper is organized as follows. In the first section, mappings for the direct and inverse problem are discussed for a model in which the NKO orbital plane is parallel to the equatorial plane. In the second section, forward mappings are discussed for a model in which no assumptions are made on the NKO orbital plane. Numerical test cases for both aforementioned models are shown and discussed in the third section. Conclusions and final remarks are drawn in the last section. VERTICAL-DISPLACEMENT MODEL Circular non-Keplerian orbits can be obtained considering the dynamics of a low-thrust propelled spacecraft in a rotating reference frame.2 The angular velocity of rotation of the reference frame , together with the out-of-plane displacement z and the radius of the orbit , are used as free parameters of the problem. Stationary solutions of the equations of motion in this rotating reference frame correspond to periodic, displaced, circular orbits with the orbital plane parallel to the Xˆ Yˆ plane when viewed from an inertial reference frame. In this paper, the Earth-Centered Inertial frame Xˆ , Yˆ , Zˆ is considered. The ac-

celeration needed to generate such stationary solutions can be derived in a closed, analytical form. 2 The thrust vector lies in the plane spanned by the radius vector and the axis of symmetry, as shown in Figure 1. Following Reference 2, it can be shown that the propulsive acceleration vector can be described, as a function of the NKO elements xNKO z, , , by its magnitude a and pitch angle , as follows. T

2 a z , , 2 2 *2 z 2 *4 2 tan z , , 1 z *

(1)

The term * r 3 in Eq. (1) is the orbital angular velocity of a Keplerian orbit of radius r 2 z 2 . Viewed from an inertial reference frame, the orbits generated by the acceleration described in Eq. (1) correspond to circular orbits displaced above the central body. The angular velocity of the rotation of the reference frame corresponds to the angular velocity of the circular orbit viewed from an inertial reference frame. Three families of NKOs can be defined based on the value of the angular velocity of the rotating reference frame.2 Type 1 NKOs are defined as those orbits with the minimum required acceleration. From Eq. (1), the requirement of minimum acceleration leads to an angular velocity of the rotating reference frame for Type 1 NKOs defined as shown in Eq. (2). A second family of NKOs is characterized by orbits synchronous with a body on a circular Keplerian orbit in the z 0 plane with the same orbit radius. Lastly, a third family of NKOs is defined by setting the orbital period to a fixed value. The angular velocities that characterize the aforementioned three families of NKOs are shown in Eq. (2).

Type 1 r 3 , Type 2 3 , Type 3 0

2

(2)

Figure 1. Displaced orbit with thrust-induced acceleration for the vertical-displacement model. In this section, the mappings between NKO elements and the osculating orbit are derived in closed, analytical form for both the direct and inverse problem. Three sets of orbital elements are chosen to describe the osculating Keplerian orbit: a) classical Keplerian elements (KEP) xKEP a, e, i, , , , b) modified T

equinoctial elements (MEE) xMEE p, f , g , h, k , L , and c) augmented integrals of motion (AIoM) T

T

x AIoM hT , eT , L (with h e 0 ). Therefore, the formulation of the spacecraft thrust-induced acceleration given in Eq. (1) is derived as a function of the osculating orbital elements. Forward Map: Non-Keplerian Orbit to Osculating Orbital Elements Starting from the NKO elements xNKO z, , , the Cartesian position vector r and velocity vector T

v are given by

r cos t sin t z , v sin t cos t 0 T

T

(3)

From the Cartesian position and velocity vectors described in Eq. (3), the forward maps are derived in the following sections for each one of the osculating orbital elements considered. Moreover, a sensitivity analysis has been performed in order to study the behavior of the mappings in presence of errors in the NKO elements. Both analytical and numerical sensitivity analyses have been carried out to confirm the results found. The analytical sensitivity analysis is centered on the computation of the Jacobian J of the mapping. NKO to Osculating KEP Map. From Eq. (3), the NKO to osculating classical Keplerian elements map can be derived following the conversion from Cartesian position and velocity to KEP, as described in Algorithm 4.2 from Reference 8. The resulting map is shown in Eq. (4).

It is worth noting that the case z 0 2 2 2 z 2

is not taken into account in the formulation

of and in Eq. (4). However, this case corresponds to an osculating circular orbit, which is not possible for the types of NKOs taken into account within this study. Moreover, it is important to underline that 0 is arbitrarily chosen in the case of a planar NKO. Lastly, Eq. (4) shows that 0 only if z 0 2 3 , which is the condition for a Keplerian orbit. The first important characteristics of the osculating Keplerian orbits that describe an NKO have already been shown in the forward map. In fact, from the formulation of the true anomaly in Eq. (4), it can be noted that a spacecraft on a NKO is always either at the apogee or perigee of the osculating Keplerian orbit.

3

a

2 z2 2 2 2 2 z 2

i arccos 2 z2

e

,

2 2 2 z 2

0 if z 0 arccos sin t sgn z if cos t 0 z 0 cos t 0 z 0 2 arccos sin t sgn z otherwise if z 0 2 3

0 t t 2 3 2

if z 0 2 3 if z 0 2 3

if z 0

z z 0

if z 0 2 2 2 z 2 z 0 2 2 2 z 2 2

2

2

2

2

2

2 z2

(4)

t if z 0 2 3 0 if 2 2 2 z 2 2 2 2 2 if z

For the sensitivity analysis, an analytical study can be carried out only on those elements characterized by a continuous function (i.e. semi-major axis, eccentricity, and inclination). Only a numerical sensitivity analysis can be carried out on the last three elements of the map. The parameters used to test the sensitivity of the mapping to uncertainties in an NKO element xi depend on the classical Keplerian element and are (Eq. (5)): a) relative error of the semi-major axis; b) absolute error of the eccentricity; and c) normalized error of the general angles i, , , .

a , x : i

J a , xi xi a

, e, xi : J e, xi xi , , xi :

mod J , xi xi ,

(5)

NKO to Osculating MEE Map. In order to have a mapping with no singularities, a map from Cartesian position and velocity to osculating modified equinoctial elements is derived. In fact, MEE have been introduced to avoid singularities occurring in Keplerian orbits in the case of planar and/or circular orbits and are defined as follows.9

p a 1 e2

f e cos

h tan i 2 cos

k tan i 2 sin L

g e sin

(6)

Starting from the expressions of the Cartesian position and velocity vectors (Eq. (3)) and the analytical mapping of the osculating classical Keplerian elements shown in Eq. (4), the mapping from NKO elements to osculating MEE is shown in Eq. (8). The full derivation of the NKO to osculating MEE map is not shown here for the sake of conciseness. It is worth noting that there are no singularities in the formulations of the out-of-plane MEE in the case of zero out-of-plane displacement orbits. It can be seen, in fact, that

lim h lim k 0 z 0

z 0

4

(7)

p

2 2 2 z2

h

sin t 2 z 2

z

cos t 2 z 2 2 2 2 z2 f cos t k z g

2 2 2 z2 sin t

(8)

L t

Moreover, it is important to note that no assumptions have been made in the derivation of this map. Indeed the functions that describe the NKO to osculating MEE map are continuous in the whole spectrum of the NKO elements. The parameters used to test the sensitivity of the mapping to uncertainties in an NKO element xi depend on the MEE and are (Eq. (9)): a) relative error of the semi-latus rectum; b) absolute error of the general element y f , g , h, k ; and c) normalized error of the true longitude L .

p , x : i

J p , xi xi p

, y , xi : J y , xi xi , L , xi :

mod J L , xi xi ,

(9)

NKO to Osculating AIoM Map. A third set of osculating orbital elements has been chosen for the map. The integrals of motion (orbital angular momentum h and eccentricity vector e ) are computed by using the expressions of Cartesian position and velocity vectors given in Eq. (3). The resulting map is as follows:

z cos t 22 1 h r v z sin t , e 2 z2 2

cos t sin t , L t z

(10)

The parameters used to test the sensitivity of the mapping to uncertainties in an NKO element xi depend on the AIoM and are (Eq. (11)): a) relative error of the magnitude of the difference between the nominal angular momentum vector and the one related to uncertainties in the NKO elements; b) absolute error of the magnitude of the difference between the nominal eccentricity vector and that related to uncertainties in the NKO elements; and c) normalized error of the true longitude L .

h, x : i

J h, xi xi h

, e , xi : J e , xi xi , L , xi :

mod J L, xi xi ,

(11)

Inverse Map: Osculating Orbital Elements to Non-Keplerian Orbit After the forward mappings from non-Keplerian orbit families to osculating orbital elements described in the previous section, three different inverse mappings from orbital elements to NKO geometry have been generated. These mappings have been derived in closed, analytical form. A sensitivity analysis has been performed to understand the impact of uncertainties in the orbital elements on the NKO properties. Only a numerical sensitivity analysis has been carried out by computing the values of the osculating elements from the perturbed input elements and then computing the errors respect to the nominal values. The parameters used to test the sensitivity of the mapping to uncertainties in the osculating orbital elements are as follows.

Absolute error of the vertical displacement. The absolute error has been chosen against the relative error because all the NKOs are feasible with the near-term technology are characterized by small vertical displacements.

5

Relative error of the orbit radius. The relative error has been chosen against the absolute error because usually NKO with small vertical displacements are characterized by large orbit radii.

Relative error of the orbital angular momentum.

Osculating KEP to NKO Map. Recalling that the radius of the osculating orbit can be expressed, in terms of KEP, as r a 1 e2 1 e cos , the position and velocity vectors in the orbital reference frame

rˆ,ˆ, hˆ are given by

rorb r cos r sin 0 , vorb r cos r sin r sin r cos 0 T

T

(12)

in which

a 1 e2 h r2 r2 2 2 a 1 e ea 1 e sin 2 r t 1 e cos 1 e cos

a 1 e2

(13)

e sin

Therefore, the position and velocity vectors shown in Eq. (12) can be written in the Earth-Centered Inertial reference frame by following the rotations R3 R1 i R3 . After algebraic manipulation, the map from osculating classical Keplerian elements to NKO elements is summarized as

a 1 e2

sin i sin 1 e cos a 1 e2 3 cos 2i 2 cos 2 sin 2 i 2 1 e cos z

a 3 1 e 2

2 1 e cos 3

(14)

1 e 2 2e cos 3 cos 2i 2 cos 2 sin 2 i

Osculating MEE to NKO Map. From the formulation of Cartesian position and velocity vectors as a function of MEE given in Betts,10 the inverse map from osculating MEE to NKO elements is derived as

2r h sin L k cos L s 2 r 1 h 2 k 2 2 cos 2 L 4hk sin 2 L s z

s r

(15)

1 f 2 g 2 2 f cos L 2 g sin L

p 1 h 2 k 2 2 2 cos 2 L 4hk sin 2 L

in which

r

p , s 1 h2 k 2 , h2 k 2 1 f cos L g sin L

6

(16)

Osculating AIoM to NKO Map. From Eq. (4), it was noted that the spacecraft is always at either the apogee or perigee of its osculating orbit. Therefore, the magnitude of the position vector can be expressed, in terms of the osculating integrals of motion, as

r

h

2

1 e2

1 e

(17)

in which the ambiguity in the sign is due to the possibility of the spacecraft being at either the apogee or perigee of the osculating orbit. Defining the four-quadrant inverse tangent as* arctan X , Y tan 1 Y X 2 sgn Y 1 sgn X , this condition can be expressed, in terms of AIoM, as

apogee if mod arctan ex , ey , 2 L perigee if arctan ex , ey L

(18)

v2 1 Recalling that e r , h r v , and after some manipulation, the Cartesian position and veloc r ity vectors can be written, in terms of AIoM, as 2 h eˆ if mod arctan ex , ey , 2 L 1 e r 2 h eˆ if arctan ex , ey L 1 e

ˆ h 1 e h eˆ v 1 e hˆ eˆ h

(19)

if mod arctan ex , ey , 2 L if arctan ex , ey L

(20)

Therefore, the inverse map from osculating AIoM to NKO elements is derived as shown in Eq. (21). From Eq. (21), it can already be noted that such map is extremely sensitive to errors in both eccentricity vector and true longitude. In fact, if the spacecraft is computed to be far from both apogee and perigee because of errors in the osculating AIoM, one is not able to convert the osculating orbital elements to NKO elements by using this map.

*

Definition from MATLAB User’s Guide, Ver. R2014b, atan2 – Symbolic four-quadrant inverse tangent.

7

2 h ez 1 e z 2 h ez 1 e

h2 1 e 2 h 1 e

if mod arctan ex , ey , 2 L if arctan ex , ey L

ex e

2

ey e

ex e

2

ey e

if mod arctan ex , ey , 2 L

2

(21)

2 1 e 2 3 h 2 2 1 e 3 h

if arctan ex , ey L

2

ex e ex e

2

ey e

if mod arctan ex , ey , 2 L

2

if arctan ex , ey L

1 ey e

2

1 2

Determination of Spacecraft Thrust-Induced Acceleration The magnitude and direction of the spacecraft thrust-induced acceleration used to generate families of NKOs have been derived by using the inverse mapping from the orbital elements discussed in the previous section as inputs for Eq. (1). This analysis links the thrust-induced acceleration to the NKO geometry. The analytical expressions for the thrust-induced acceleration have been derived for osculating KEP and MEE. Because of the very poor robustness of the osculating AIoM for the inverse mapping, these elements have not been considered in the study of the thrust-induced acceleration. Classical Keplerian elements. From Eq. (1) and using the mappings shown in Eq. (14), the magnitude of the thrust-induced acceleration and pitch angle are obtained as

a xKEP

tan xKEP

1 e cos 2a 2 1 e2

2

A B

(22)

4 1 e2 2e cos 1 1 2sin i sin 1 1 e cos

(23)

in which

A 4 1 e cos 2 sin 2 i sin 2 2 B 1 1 e cos 4 1 e 2e cos 1 1 3 cos 2i 2 cos 2 sin 2 i

2

(24)

Modified equinoctial elements. From Eq. (1) and using the mappings shown in Eq. (15), the magnitude of the thrust-induced acceleration and pitch angle are obtained as shown in Eqs. (26) – (27), in which s and are defined in Eq. (16), whereas

q 1 f cos L g sin L 2 2 2 2 1 h k 2 cos 2 L 4hk sin 2 L

8

(25)

a xMEE

q sp 2

s 2 1 f 2 g 2 2 f cos L 2 g sin L 4q k cos L h sin L 2 q 2 2

2

tan xMEE

s 2 1 f 2 g 2 2 f cos L 2 g sin L 1 2 h sin L k cos L q 2

2

2

(26)

(27)

ARBITRARY-ORIENTATION MODEL The mappings derived in the previous section are related to non-Keplerian orbits whose orbital plane is displaced vertically and, therefore, parallel to the Xˆ Yˆ plane in a Cartesian reference frame (i.e. the equatorial plane, in this work). A mapping for the direct problem is now investigated which considers an arbitrary orientation of the NKO plane.

Figure 2. Displaced orbit with thrust-induced acceleration for the arbitrary-orientation model. Two extra elements are introduced for the description of the non-Keplerian orbit in the threedimensional space, as shown in Figure 2. These are

j Note that the symbols

Inclination of the non-Keplerian orbit Right ascension of the ascending node of the non-Keplerian orbit

j,

(28)

have been used for the description of the NKO elements inclination and

right ascension of the ascending node (RAAN) because these elements should not be confused with the osculating inclination and RAAN, respectively. Therefore, the set of NKO elements is

xNKO z, , , j, . The formulation of the Cartesian position and velocity vectors shown in Eq. (3) can T

be expanded, considering inclination and RAAN, as

cos cos t cos j sin sin t z sin j sin r sin cos t cos j cos sin t z sin j cos sin j sin t z cos j

(29)

cos sin t cos j sin cos t v sin sin t cos j cos cos t sin j cos t

(30)

9

Following the same procedure used to obtain the mappings in the vertical-displacement case, the arbitrary-orientation mappings from NKO elements to osculating orbital elements are obtained. Because of the very poor robustness of the osculating AIoM for the inverse mapping, these elements have not been considered in this study. NKO to Osculating KEP Map. From Eqs. (29) – (30), the NKO to osculating classical Keplerian elements map can be derived following the conversion from Cartesian position and velocity to KEP, as described in Algorithm 4.2 from Reference 8. The resulting map is as follows.

a

2 z2 2 2 2 2 z 2

e

,

2 2 2 z 2

cos j z sin j sin t 0 , * i arccos 2 2 z 0 t t * * 2

if z 0 j 0 otherwise

if 2 2 2 z 2 if j 0 z 0 2 3 if j 0 z 0 2 3

if 2

2

z 2

2

(31)

z cos j sin j sin t 0

otherwise

t if z 0 2 3 0 if 2 2 2 z 2 2 2 2 2 if z in which

arccos N x N if N y 0 * , 2 arccos N x N otherwise

N z 2 cos 2 t sin j z cos j sin t

N x sin j cos z cos j cos sin t sin cos t , * arccos sgn 2 2 2 z 2 N y sin j sin z cos j sin sin t cos cos t ,

2

sin j cos t 2 z 2 z 2 cos 2 t sin j z cos j sin t

2

(32) Note that the mapping shown in Eq. (4) is a particular case of the mapping for an arbitrary orientation of the orbital plane. In fact, it can be easily verified that Eq. (31) becomes Eq. (4) if j 0 . NKO to Osculating MEE Map. Starting from the expressions of the Cartesian position and velocity vectors (Eqs. (29) – (30)) and the analytical mapping of the osculating classical Keplerian elements shown in Eq. (31), the arbitrary-orientation mapping from NKO elements to osculating MEE is shown in Eq. (33), in which * and are defined in Eq. (32). The full derivation of the NKO to osculating MEE map is not shown here for the sake of conciseness. Note that, as in the previous case, the mapping shown in Eq. (8) is a particular case of the arbitraryorientation mapping. In fact, it can be easily verified that Eq. (33) becomes Eq. (8) if j 0 . Because of the introduction of inclination and RAAN related to the non-Keplerian orbit, in this case, it has not been possible to simplify further the mapping and the piecewise formulation cannot be avoided.

10

p

f

2 2 2 z2 2

2

2

2

g

h

k

cos t if z 0 j 0 z * cos arccos if z cos j sin j sin t 0 * cos arccos if z cos j sin j sin t 0 2

2

sin t if z 0 j 0 z * sin arccos if z cos j sin j sin t 0 * sin arccos if z cos j sin j sin t 0 2

2

2 z 2 cos j z sin j sin t sin j cos z cos j cos sin t sin cos t 2 z 2 cos j z sin j sin t

z 2 cos 2 t sin j z cos j sin t

2

2 z 2 cos j z sin j sin t sin j sin z cos j sin sin t cos cos t 2 z 2 cos j z sin j sin t

z 2 cos 2 t sin j z cos j sin t

t if z 0 j 0 * L arccos if z cos j sin j sin t 0 * arccos if z cos j sin j sin t 0

2

(33)

NUMERICAL TEST CASES Several test cases are chosen to assess the validity of the mappings, investigate the behavior of the osculating elements and understand the impact of uncertainties on the mappings. For the sake of conciseness, only the most interesting results are shown here. In the following sections, the test cases related to both the vertical-displacement and arbitrary-orientation models are shown and discussed. Vertical-Displacement Model A Geostationary Earth Orbit (GEO) is chosen as the reference Keplerian orbit. In terms of NKO elements, a GEO is characterized by11

z0 GEO : rGEO 42,164.1696 km 3 rGEO

(34)

Two displaced GEOs are chosen as test cases in this study: a Type 1 NKO and an in-plane displaced Type 3 NKO. Station-keeping regulations currently require a station-keeping box of 0.1 0.2 deg .11 Therefore, the out-of-plane displacement for the Type 1 NKO is chosen such that the spacecraft hovers just above the GEO station-keeping box. The same displacement is considered in the in-plane case. The NKO elements of the two test cases are

z 35 km Type 1 NKO : rGEO cos z rGEO 3 GEO rGEO

Type 3 NKO rGEO , : (in-plane)

11

z0 rGEO 35 km 3 GEO rGEO

(35)

Forward Maps. Figure 3 shows the Type 1 NKO with z 0 described by Eq. (35), together with its osculating Keplerian orbits, corresponding to the instantaneous orbital elements of the spacecraft if the thrust is nulled. In this case, it is shown how the semi-major axis and eccentricity of the osculating orbits do not change, but the osculating orbit plane rotates around the vertical axis. The spacecraft, being always at the apogee of the osculating Keplerian orbits, describes the desired NKO. Therefore, the envelope resulting from the osculating Keplerian orbits after an entire orbit is well approximated by a truncated cone, the characteristics of which depend on inclination and apogee and perigee radii of the osculating orbits. Moreover, from Figure 3, another characteristic of the NKO to KEP map for a Type 1 NKO can be seen. That is, is the only osculating KEP that varies in time. All other osculating KEP are constant, as shown in Eq. (4). On the other hand, the in-plane displaced Type 3 NKO (not shown here for the sake of conciseness) is characterized by an osculating argument of perigee being the only element that varies in time. Figure 4 shows the double cone created by the osculating eccentricity unit vectors due to the rotation of the osculating Keplerian orbits. The height of the cone and its base radius depend on both the vertical displacement z and the radius of the orbit . The osculating eccentricity vector depend on the sign of the vertical displacement, whereas the osculating angular momentum does not. On the contrary, the osculating eccentricity vectors for the in-plane Type 3 NKOs lie on the Xˆ Yˆ plane, while the eccentricity vector in the GEO case is arbitrarily set to eGEO 1 0 0 . T

Figure 3. Type 1 NKO displaced GEO (solid black line) and osculating Keplerian orbits.

Figure 4. Time evolution of the osculating eccentricity unit vectors.

To investigate the sensitivity analysis of the forward maps, the relative error of the NKO elements is considered fixed and constant, such that z 103. Considering the characteristics of the orbits shown in Eq. (35), the absolute errors on the NKO elements considered for the sensitivity analysis are

z1 z z1 0.15 km Type 1 Type 3 NKO : 1 1 42.2 km , : NKO (in-plane) 1 1 0.36 deg day

z3 z z1 0.15 km 3 3 42.2 km 0.36 deg day 3 3

(36)

Note that the term z3 in Eq. (36) depends on the displacement of the Type 1 NKO. This is made on purpose in order to have a non-zero error for the out-of-plane displacement. Figure 5 shows the evolution over one orbit of the osculating in-plane MEE f for the case of Type 1 NKO. The nominal value and the values due to errors in the NKO elements are shown. It can be seen that the error on the in-plane element f due to uncertainties in the NKO elements is of the same order of mag-

12

nitude as the uncertainties themselves. The same behavior is noted in the other in-plane modified equinoctial element g . On the other hand, Figure 6 shows the evolution over one orbit of the out-of-plane MEE h for the case of Type 1 NKO. The nominal value and the values due to errors in the NKO elements are shown. Figure 6 shows that the error on the out-of-plane element h due to uncertainties in the NKO elements is orders of magnitude smaller than the value of the uncertainties. The same behavior is noted in the other out-of-plane modified equinoctial element k . Here, the semi-latus rectum shows a constant error due to uncertainties in the NKO elements of the same order of magnitude as the uncertainties themselves. Lastly, the true longitude exhibits a linearly increasing error due to a 0.1% uncertainty on the orbit angular velocity. After one orbit, the absolute error on the true longitude is therefore 0.36 deg .

Figure 5. Evolution of the f over one orbit.

Figure 6. Evolution of h over one orbit.

Inverse Maps. The same test cases used to study the forward maps are now considered for the study of the inverse maps. All the results of the mappings have been numerically verified against their nominal values. That is, the osculating orbital elements have been obtained through the forward mappings starting from the nominal NKO elements. Then, the inverse mappings are used to obtain the original NKO elements. The errors between the nominal values of the NKO elements and those obtained after the conversion gives an estimate of the accuracy of both forward and inverse mappings. This results in errors of less than 109 km for z and , and errors of less than 1012 deg day for . This validates the inverse maps. For the sensitivity analysis, the errors of the osculating KEP and MEE are considered fixed and constant, as follows.

a 103 (relative) KEP : e 103 (absolute), MEE : 3 {i ,, ,} 10

p 103 (relative) 3 { f , g , h, k } 10 (absolute) L 103

(37)

in which the errors relative to angular quantities are considered as 0.1% of the maximum possible error ( ). As discussed earlier, the AIoM to NKO elements map is extremely sensitive to errors in both eccentricity vector and true longitude. The issue is not only the piecewise formulation per se, but the fact that the formulation takes into account only two points in the orbit, which are the apogee and perigee. If an error in either eccentricity vector or true longitude moves the point from the apogee/perigee, the inverse map fails. Therefore, there is no need to perform a sensitivity analysis for the AIoM case. The sensitivity analysis study demonstrates that the planar elements (i.e. , ) of both mappings are robust to uncertainties in the osculating orbital elements. Both mappings of the out-of-plane displacement z are robust to uncertainties in the in-plane osculating elements. However, the out-of-plane displacement is very sensitive to errors in the out-of-plane osculating elements. This is understandable if the osculating

13

Keplerian orbit has a small inclination and a large semi-major axis ( a ~ rGEO in this case). For the sake of conciseness, only the most interesting plots of the sensitivity analysis study are shown here. Figure 7 shows the relative error evolution of the orbit radius over one orbit obtained from the osculating MEE with a 0.001 error in the in-plane MEE f . Both test cases are shown. It is shown that the relative error on the orbit radius is at most the same order of magnitude as the initial uncertainty itself. Figure 8 shows the evolution of the out-of-plane displacement z over one orbit obtained from the osculating MEE. Both test cases are shown. Both the nominal value and a 0.001 error in the out-of-plane MEE h are considered. It is shown how a small error in the out-of-plane element can cause a large error in the out-of-plane displacement, as discussed above.

Figure 7. Relative error in orbit radius in case of 0.001 error in f.

Figure 8. Out-of-plane displacement. Nominal case and 0.001 error in h.

Thrust-induced acceleration. The acceleration, in terms of magnitude and pitch angle, resulting from both Eqs. (22) – (23) and Eqs. (26) – (27) has been validated using two separate checks. In the first check, the values of magnitude and pitch angle are compared with those given directly by using Eq. (1). The values of both the magnitude and direction of the acceleration have been verified for all test cases considered with both sets of elements. The errors in magnitude and direction of acceleration are a 1011 mm s2 and 1012 rad , respectively. The second check has been carried out by propagating the equations of motion with the propulsive acceleration given by Eqs. (22) – (23) and Eqs. (26) – (27), respectively. The relative errors in position and velocity vectors between the initial state and the state after one orbital period are used to validate the acceleration. A variable order Adam-Bashford-Moulton PECE solver (as implemented in MATLAB ode113) with relative and absolute tolerances equal to 1010 and 1012 , respectively, has been used for the propagation of the equations of motion. The relative errors in the final position and v v 109 . velocity are r r Lastly, the accelerations required for the test cases taken into account have magnitudes a 0.19,0.56 mm s2 for Type 1 and Type 3 NKOs, respectively. Considering a 1000-kg spacecraft, such accelerations can be easily converted into maximum required thrusts Tmax 0.19,0.56 N that are feasible values for near-term low-thrust technology.12 Summary. The use of classical Keplerian elements to describe the osculating orbits is impractical for the forward map mainly because of the piecewise formulation of most of the elements. On the other hand, KEP are good candidates to be used for the inverse map because of their easy formulation and robustness to uncertainties. The use of modified equinoctial elements to describe the osculating orbits is a good choice due to the easy formulation and robustness to uncertainties. Moreover, for the same reasons, MEE are good candidates to be used for the inverse map as well. The use of augmented integrals of motion to describe the

14

osculating orbits resulted is also a good choice due to the easy formulation of the map, even if the elements are less intuitive than the other sets. On the other hand, the formulation of the inverse map by means of AIoM is extremely sensitive to errors in both eccentricity vector and true longitude to the point that the map fails if there is an error in either of these elements. In conclusion, the use of MEE guarantees an easy and robust formulation for both forward and inverse mappings in the case of the vertical-displacement model, with the NKO orbital plane parallel to the equatorial plane. Lastly, the analytical formulations of the thrust-induced acceleration starting from both KEP and MEE have been demonstrated to be valid. A summary of advantages and drawbacks of the mappings for both the direct and inverse problems is shown in Tables 1 and 2, respectively, for each set of osculating elements considered. Table 1. NKO to osculating orbital elements map. Summary of advantages and drawbacks. KEP

Advantages

Easy to understand Only one element is time dependent

Drawbacks

Piecewise formulation Very sensitive to uncertainties Issues around z 0 Singularity: 0 if z 0 Very sensitive to the switch apogee/perigee

MEE Easy formulation (no piecewise cases) Robust to uncertainties Robust to the switch apogee/perigee No singularities 5 elements are time dependent

AIoM Easy formulation No singularities

Less intuitive 7 elements instead of 6 Very sensitive to the switch apogee/perigee (eccentricity vector changes sign)

Table 2. Osculating orbital elements to NKO map. Summary of advantages and drawbacks. KEP Easy to understand Advantages Easy formulation Robust to uncertainties Out-of-plane displacement very sensitive to uncertainties on inDrawbacks clination

MEE Easy formulation Robust to uncertainties on in-plane elements Out-of-plane displacement very sensitive to uncertainties on out-of-plane elements

AIoM

Piecewise formulation Extremely sensitive to errors in eccentricity and true longitude

Arbitrary-Orientation Model Three scenarios have been tested in order to both verify the validity of the arbitrary-orientation forward mappings and study the behavior of the osculating orbital elements. The three test cases have been chosen to be feasible mission scenarios in the near future. For this reason, only small-displacement NKOs are considered which can be generated with near-term technology. The first scenario consists of a 35-km out-of-plane (Type 1) and in-plane (Type 3) displaced GEO, as described in the case of the vertical-displacement model. The second scenario is a 5-km out-of-plane (Type 1) and in-plane (Type 3) displaced global positioning system (GPS) orbit. In this scenario, the spacecraft hovers above/below or inside/outside the GPS spacecraft for visual inspection. Moreover, a distance of 5 km can be considered within the range for proximity operations. 13 The orbit considered for describing the GPS orbit is related to the PRN-06 satellite, and its orbital parameters are accurate with respect to the Yuma Almanac 2016* (week 0896, GPS Time of Applicability 405504 s). PRN-06 has been chosen because is

*

Data available online at https://celestrak.com/GPS/almanac/Yuma/2016/ [retrieved 25 October 2016].

15

the one with the lowest value of eccentricity ( e 0.22 103 ) among those reported in the Yuma Almanac and, therefore, a circular orbit is a good approximation. The third and last scenario taken into account is a 1-km out-of-plane (Type 1) and in-plane (Type 3) displaced Sun-synchronous orbit (SSO). As for the previous scenario, a 1-km range can be used for visual inspection of the SSO satellite and is within the range for proximity operations. In terms of NKO elements, the characteristics of the GEO orbit are defined in Eq. (34) together with j 0 . The GPS and SSO orbits are characterized, in terms of NKO elements, by

z0 rGPS 26,560.9478 km 3 GPS : rGPS , SSO : j 55.2885 deg 77.7881 deg

z0 rSSO 7378.16 km 3 rSSO j 99.4845 deg 0

(38)

Signatures. The aforementioned test cases have been used to seek key characteristics of signatures that can provide a simple and reliable indication that a spacecraft is being forced along a non-Keplerian orbit. Both KEP and MEE are characterized by four elements that can be used to identify a NKO. Therefore, assuming perfect knowledge of the NKO geometry, the elements that show a peculiar trend are shown in Eq. (39), in which is the semi-amplitude of the cosine wave related to the element and is a generic phase.

i t i cos t t 0 t KEP : t cos t t t 0 t 0

if z 0 j 0

f t g t if z 0 j 0 , MEE : h t if z 0 j 0 k t

if z 0 j 0

f cos t g cos t h cos t

if z 0

k cos t

if z 0

(39)

Impact of noise. The impact of noise on the NKO elements is considered here to study the signatures of NKOs in the presence of noise on orbit determination. The same scenarios used to study the signatures have been considered. However, only the Type 1 NKOs with a positive displacement are considered. Three values for the noise have been considered in this study, consisting of upper/lower boundaries and a value inside this range. For the upper boundary, the maximum range error r 24 m and range rate error r 0.16 m s , due to the instrumental errors, are considered. 14 A second set of values for the noise is considered by dividing the upper boundary by a factor of 10. It will be seen that these values lie within the range defined by upper and lower boundaries. A lower boundary on the noise can be set by considering the errors in the orbit determination of the test cases. Because these errors are the result of filtering processes (i.e. a Keplerian orbit is assumed in the first place), these values are considered here as the lower boundary, as follows.15-17

r GEO : i

0.1 m

r , GPS/SSO : v deg

103 deg 2 10

3

1m 2 mm s

(40)

A Monte Carlo analysis has been carried out to obtain the maximum values of the errors on the Keplerian inclination i and RAAN from the aforementioned errors. 500 samples of Cartesian position and velocity vectors in a normal distribution have been considered with mean values corresponding to the nominal Cartesian states. The standard deviation considered for each component of the position vector is r , assuming the range error being equal to r . The standard deviation considered for each component of the

16

velocity vector is v , assuming the range rate error being equal to v . Considering the conversion from Cartesian state to KEP, the errors given by the aforementioned values provide an estimate of i and

. The three sets of errors considered here are, therefore, the following. r 24 m Case 1 : i 0.01 deg , Case 2 : (upper bound) 180 deg if GEO 0.01 deg otherwise : GEO: (lower bound) Case 3

r 2.4 m i 0.001 deg if GEO 180 deg (41) 0.001 deg otherwise

r 0.1 m r 1 m 3 4 i 10 deg , GPS/SSO: i 10 deg 2 103 deg 10 4 deg

It is important to note that 180 deg in the cases 1-2 related to GEO. This is due to the planar nature of GEO and, therefore, the fact that RAAN is arbitrarily set to zero. Even a small error can introduce a small inclination and, therefore, a value of 0, 2 . On the other hand, the value of the error on RAAN for case 3 is significantly lower because of the filtering. The errors

r

, i , are used to estimate

the noise to be considered for the NKO elements. A Monte Carlo analysis has been carried out with 500 samples of NKO elements in a normal distribution. The mean values considered are the nominal values of NKO elements

z0 , 0 , j0 , 0 ,

whereas the standard deviations are

r sin z

0

0 , 3 r , i , ,

respectively. The value of the angular velocity is simply given by 3 , in which is the NKO orbit radius with noise. Figures 9 and 10 both show GEO and Type 1 NKOs in the presence of noise as from Eq. (41) for cases 1 and 2, respectively. Both figures have been generated using the Cartesian position vectors resulting from the osculating KEP obtained through the forward map from NKO elements with noise. The figures show that the Keplerian GEO and the NKO are clearly distinguishable for the value of noise related to Case 2. On the contrary, the same is not true for the three-dimensional orbits of both the GPS and SSO cases. That is, in the GPS and SSO cases, one cannot understand if the orbit is Keplerian or non-Keplerian. Figure 11 shows a comparison between all scenarios and noise levels. The evolution of the out-of-plane MEE k has been chosen for comparison. It can be seen that the first level of noise is too large to be able to distinguish Keplerian from non-Keplerian orbits. On the other hand, with the third level of noise, which has been considered as the lower bound, all test cases show a clear division between the two types of orbits. If the second level of noise is considered, the signatures of a spacecraft being forced along an NKO are still clear in the GPS and SSO cases. In fact, the constant trend of k that is peculiar to a Keplerian orbit is clearly different to the sinusoidal trend of k , which is characteristic of an NKO. On the contrary, the GEO case is still too noisy to be able to separate the two trends. This is due to the very large error on RAAN considered (Eq. (41)). Lastly, it is interesting to study the thrust-induced acceleration needed to follow the osculating elements in presence of noise. However, the formulation of the acceleration as a function of the osculating elements shown in Eqs. (22) – (27) has been derived for a reference frame so that viewed from an inertial frame of reference, such solutions correspond to circular orbits displaced above the central body. 2 That is, the reference frame should be rotated following the inclination and RAAN of the nominal Keplerian orbit before computing the magnitude of the acceleration. Assuming a knowledge of the nominal Keplerian orbit, the reference frame can be rotated by the RAAN and inclination following the rotations R3 R1 j . This

17

is a good assumption since the data are noisy around the nominal values. Therefore, the values of RAAN and inclination can be considered as the mean values of the data available.

Figure 9. GEO and NKO. Noise level case 1.

Figure 10. GEO and NKO. Noise level case 2.

a) GEO. Noise level case 1.

b) GEO. Noise level case 2.

c) GEO. Noise level case 3.

d) GPS. Noise level case 1.

e) GPS. Noise level case 2.

f) GPS. Noise level case 3.

g) SSO. Noise level case 1.

h) SSO. Noise level case 2.

i) SSO. Noise level case 3.

Figure 11. Evolution of k over one orbit. All mission scenarios and noise levels are shown.

18

a)

b) Figure 12. Magnitude of acceleration over one orbit for noise level case 2. a) GEO. b) GPS.

This way, the rotated Cartesian state describes a noisy orbit with approximately zero inclination. Therefore, the analytical formulation of the acceleration given either osculating KEP or MEE shown in Eqs. (22) and (26) can be used to provide a first-guess approximation of the magnitude of the acceleration needed to follow the observed NKO. Moreover, the value of the magnitude of acceleration can be a further indication of the use of a NKO. Figure 12 shows the magnitudes of acceleration computed considering the second level of noise for the GEO and GPS cases. Summary. In Table 3, the impact of noise on the mappings is summarized. Those signatures that, despite the noise, can provide a clear and reliable indication that a spacecraft is being forced along a nonKeplerian orbit are highlighted for each case study under consideration. Table 3. Summary of the impact of noise on the NKO signatures. Signatures GEO and NKO are clearly defined from the following (Cases 2-3 only): 3D orbits GEO Osculating inclination a (vertical-displacement analytical approach) GPS/SSO and NKO are clearly defined from the following (Cases 2-3 only): Elements i, , h, k

GPS SSO

a (vertical-displacement analytical approach)

Issues Case 1 is too noisy From , , and all MEE, GEO and NKO are not clearly defined

Case 1 is too noisy From , , GPS/SSO and NKO are not clearly defined The magnitude of the in-plane MEE is too small GPS/SSO and NKO are not clearly defined from the 3D orbits a can be used only if j and are known.

CONCLUSIONS This paper has presented a mapping between highly non-Keplerian orbit (NKO) geometry and classical orbital elements for both the direct and inverse problem. Three sets of elements have been discussed and mappings have been derived in closed, analytical form. Two models have been considered for the NKOs: a vertical-displacement model, in which the NKO orbital plane is parallel to the equatorial plane, and an arbitrary-orientation model, in which no assumptions are made. For the vertical-displacement model, it has been shown that the main drawback of using the classical Keplerian elements (KEP) is due to the piecewise formulation required for a forward mapping. A simpler formulation characterizes, instead, a forward mapping that uses the augmented integrals of motion. However, the main drawback of using this set of elements

19

is the extreme sensitivity to uncertainties in eccentricity and true longitude that characterizes the inverse mapping. Lastly, the modified equinoctial elements (MEE) provide both a simple formulation and have low sensitivity to uncertainties in both the forward and inverse mappings. For the arbitrary-orientation model, the forward map has been derived in closed, analytical form. It has been shown that both KEP and MEE are characterized by piecewise formulations. Three test cases have been chosen that show a broad range of applicability of the maps. Using the same test cases, noise has been added to the initial NKO elements to assess the impact on estimates of the classical orbital elements. It has been demonstrated that, despite the noise, signatures exist that can provide a clear and reliable indication that a spacecraft is being forced along an NKO. Lastly, it has been shown that a first-guess approximation of the magnitude of the acceleration needed to follow the observed NKO can be obtained. Both KEP and MEE have equal advantages and drawbacks in the case of the arbitrary-orientation forward map. Therefore, in order to clearly distinguish a Keplerian orbit from a NKO using their orbital elements, it is important to look at: a) three-dimensional orbit; b) osculating i and ; c) osculating out-of-plane MEE h and k ; and d) magnitude of the thrustinduced acceleration. ACKNOWLEDGMENTS This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-16-1-0219. REFERENCES McKay, R., Macdonald, M., Biggs, J. and McInnes, C., “Survey of Highly Non-Keplerian Orbits with Low-Thrust Propulsion”, Journal of Guidance, Control, and Dynamics, Vol. 34, No. 3, 2011, pp. 645-666. doi: 10.2514/1.52133 2 McInnes, C. R., “The Existence and Stability of Families of Displacement Two-Body Orbits”, Celestial Mechanics and Dynamical Astronomy, Vol. 67, No. 2, 1997, pp. 167-180. doi: 10.1023/a:1008280609889 3 McInnes, C. R., “Dynamics, Stability, and Control of Displaced Non-Keplerian Orbits”, Journal of Guidance, Control, and Dynamics, Vol. 21, No. 5, 1998, pp. 799-805. doi: 10.2514/2.4309 4 Heiligers, J., McInnes, C. R., Biggs, J. D. and Ceriotti, M., “Displaced geostationary orbits using hybrid low-thrust propulsion”, Acta Astronautica, Vol. 71, 2012, pp. 51-67. doi: 10.1016/j.actaastro.2011.08.012 5 Macdonald, M., McKay, R. J., Vasile, M., Bosquillon De Frescheville, F., Biggs, J. and McInnes, C., “Low-ThrustEnabled Highly-Non-Keplerian Orbits in Support of Future Mars Exploration”, Journal of Guidance, Control, and Dynamics, Vol. 34, No. 5, 2011, pp. 1396-1411. doi: 10.2514/1.52602 6 McKay, R., Macdonald, M., Bosquillon de Frescheville, F., Vasile, M., McInnes, C. R. and Biggs, J. D., “NonKeplerian Orbits Using Low Thrust, High ISP Propulsion Systems”, 60th International Astronautical Congress, Daejeon, South Korea, 2009. 7 Wang, W., Yuan, J., Mengali, G. and Quarta, A. A., “Invariant Manifold and Bounds of Relative Motion Between Heliocentric Displaced Orbits”, Journal of Guidance, Control, and Dynamics, 2016, pp. 1-13. doi: 10.2514/1.G001751 8 Curtis, H., Orbital Mechanics for Engineering Students, Second Edition, Elsevier Science, 2009, pp. 209 - 212. 9 Walker, M. J. H., Ireland, B. and Owens, J., “A set modified equinoctial orbit elements”, Celestial mechanics, Vol. 36, No. 4, 1985, pp. 409-419. doi: 10.1007/bf01227493 10 Betts, J. T., Practical Methods for Optimal Control and Estimation Using Nonlinear Programming, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2010, pp. 265-267. 11 Heiligers, J., Ceriotti, M., McInnes, C. R. and Biggs, J. D., “Displaced geostationary orbit design using hybrid sail propulsion”, Journal of Guidance, Control, and Dynamics, Vol. 34, No. 6, 2011, pp. 1852-1866. doi: 10.2514/1.53807 12 Gates, M., Stich, S., McDonald, M., Muirhead, B., Mazanek, D., Abell, P. and Lopez, P., “The Asteroid Redirect Mission and sustainable human exploration”, Acta Astronautica, Vol. 111, 2015, pp. 29-36. doi: 10.1016/j.actaastro.2015.01.025 13 Howard, R. T., Bryan, T. C., Brewster, L. L. and Lee, J. E., “Proximity operations and docking sensor development”, 2009 IEEE Aerospace conference, edited by I. A. a. E. S. Society, 1287, IEEE Aerospace and Electronic System Society, Big Sky, MT, USA, 2009, pp. 1-10. 14 Kronmiller, G. C. and Baghdady, E. J., “The Goddard Range and Range Rate Tracking System: Concept, Design and Performance”, Space Science Reviews, Vol. 5, No. 2, 1966, pp. 265-307. doi: 10.1007/bf00241057 15 Knogl, J. S., Henkel, P. and Gunther, C., “Precise positioning of a geostationary data relay using LEO satellites”, 53rd International Symposium ELMAR-2011, edited by IEEE, 12306201, Zadar, Croatia, 2011, pp. 325-328. 16 Wang, F., Zhang, X. and Huang, J., “Error analysis and accuracy assessment of GPS absolute velocity determination without SA”, Geo-spatial Information Science, Vol. 11, No. 2, 2008, pp. 133-138. doi: 10.1007/s11806-008-0038-3 17 Švehla, D. and Rothacher, M., “Kinematic positioning of LEO and GPS satellites and IGS stations on the ground”, Advances in Space Research, Vol. 36, No. 3, 2005, pp. 376-381. doi: 10.1016/j.asr.2005.04.066 1

20