Solar sail equilibria with albedo radiation pressure in the circular restricted three-body problem

Accepted Manuscript Solar Sail Equilibria with Albedo Radiation Pressure in the Circular Restricted Three-Body Problem Mariusz E. Grøtte, Marcus J. Holzinger PII: DOI: Reference:

S0273-1177(16)30646-9 http://dx.doi.org/10.1016/j.asr.2016.11.020 JASR 12977

To appear in:

Advances in Space Research

Received Date: Revised Date: Accepted Date:

27 October 2015 14 November 2016 15 November 2016

Please cite this article as: Grøtte, M.E., Holzinger, M.J., Solar Sail Equilibria with Albedo Radiation Pressure in the Circular Restricted Three-Body Problem, Advances in Space Research (2016), doi: http://dx.doi.org/10.1016/j.asr. 2016.11.020

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

Solar Sail Equilibria with Albedo Radiation Pressure in the Circular Restricted Three-Body Problem Mariusz E. Grøtte1 Georgia Institute of Technology, 270 Ferst Drive Atlanta, GA 30332-0150, USA

Marcus J. Holzinger2 Georgia Institute of Technology, 270 Ferst Drive Atlanta, GA 30332-0150, USA

Abstract Solar Radiation Pressure (SRP) and albedo effects are investigated in the circular restricted three-body problem for a system consisting of the Sun, a reflective minor body and a solar sail. As an approximation of albedo radiation pressure (ARP), the minor body is treated as Lambertian with reflected flux scattered by the bidirectional reflectance distribution function. Incorporating the ARP, which is inherently a function of SRP, into solar sail equations of motion renders additional artificial equilibrium points in a volume between the L1 and L2 points which is defined as the region of influence. Based on the model, characterization of the findings are provided that are theoretically applicable to any body with discernible albedo such as for instance Earth, Mars or an asteroid. Example results are presented for a Sun-Vesta system which show that inclusion of ARP generates artificial equilibrium points requiring solar sail designs with very low mass-to-area ratio. In general the equilibrium points are found to be unstable but asymptotic stability may be enforced with sail attitude feedback control. Keywords: Solar sail; Restricted three-body problem; Albedo; Lambertian; Region of influence

Email addresses: [email protected] (Mariusz E. Grøtte), [email protected] (Marcus J. Holzinger ) 1 Graduate Student, School of Aerospace Engineering, Georgia Institute of Technology 2 Assistant Professor, School of Aerospace Engineering, Georgia Institute of Technology

1. Introduction The concept of solar sail propulsion is applicable to a wide range of missions and it is claimed that the sailing capability offers potential low-cost missions and flexible maneuvers for exploring the solar system (Macdonald and McInnes, 2011). By utilizing the solar radiation pressure (SRP) this type of spacecraft design provides the benefits of fuel-free propulsion and ability to hover at artificial Lagrange points, or equilibrium solutions, where the sail also theoretically may be transferred to configurations such as Halo orbits by changing its attitude (Farres and Jorba, 2012; McInnes et al., 1994; Howell, 1984; Farquhar and Kamel, 1973). In the restricted three-body problem (R3BP), body-fixed station-keeping maintains the spacecraft relative position which can prove to be advantageous in acquiring close high-resolution images, continuous surveillance, dropping off payload and even sampling. In a conceptual and simplified version of a solar sail in the circular restricted three-body problem (CR3BP), when modeled as a perfectly reflecting flat plate, the locations of equilibrium points mainly depend on its chosen orientation and the impact of the SRP and gravity forces. Solar sails may alter these locations by changing its attitude, trimming the sail area or even modulating its sail surface reflectivity as discussed by Gong and Li (2015); Mu et al. (2015); Aliasi et al. (2013). Burns et al. (1979) has already investigated how SRP and strength of solar tides alter the stability characteristics of artificial equilibrium points in vicinity of the Sun, inducing controlled orbits with low effort. McInnes et al. (1994) found the equilibrium points to be Lyapunov unstable in general but establishes that the system is controllable with feedback to sail attitude thus asymptotic stability may be achieved. Although there has been extensive research on solar sail dynamics with SRP alone, there has not yet been characterized and analytically defined how albedo from the minor body may affect the equations of motion (EoM). Scientific communities have shown increasing interest in primitive bodies with aim to understand their environment, value as mining-resource and role as indicators about the creation of the Solar System (Nesvorny et al., 2003; Sonter, 1997). Several recent missions illustrate this growing effort such as for instance the Hayabusa, Rosetta and Dawn missions (Hand, 2015; Taylor et al., 2015; Russell et al., 2012; Fujiwara et al., 2006; Rayman et al., 2006). Shankaran et al. (2011) provides a detailed analysis of a photogravitational R3BP of a small particle system which can be theoretically linked to smaller minor bodies in the Solar System and their emitted radiation effects 2

on the spacecraft. For reasons to model photometric signatures, several primitive bodies are approximated as Lambertian with reflected light scattered by the bidirectional reflectance distribution function (BRDF) (Torppa et al., 2003; Kaasalainen et al., 2002; Magnusson et al., 1996). However, there are limitations to approximate the albedo for larger bodies that have varying landscapes and atmospheres such as the Earth. For such cases, improved estimation of the reflected solar flux would require advanced models based on empirical data. Control methods applicable for spacecraft hovering over Near-Earth Asteroids are suggested in Broschart and Scheeres (2005). These methods solve the difficulty of orbiting small, irregularly sized and low-gravity objects. Morrow et al. (2001) investigated the solar sail dynamics in CR3BP with Hill approximation, a mathematical tool for analyzing the EoM around primitive bodies. Related to studies about solar sail dynamics in vicinity of comets, it has been shown that the sails may hover above the comet at Sun-side equilibrium points during its perihelion passage (Scheeres and Marzari, 2002). Considering a fair amount of various types of perturbations coming from primitive bodies, it is worth establishing the mathematics for how albedo affects equations of motion (EoM) and conditions for equilibrium in the CR3BP. Since solar sails are inherently sensitive to photons it is proposed in this paper to investigate the solar sail dynamics with additional forces due to albedo of a minor body. In order to establish an analytical framework to incorporate albedo radiation pressure (ARP) into the EoM, the initial step in this paper is to review the CR3BP with SRP in a Sun-minor body system. Emphasis is then made on ARP as an additional force to SRP in the EoM and conditions for equilibria are substantiated. As an approximation, although with its limitations, the minor body is modeled as a Lambertian diffuse gray body with ARP characterized by means of optical signatures. The reflected solar flux is scattered by the BRDF which is a function of the relative distance between the viewer, i.e. solar sail, and the minor body. Since the dynamics potentially are more complex in vicinity of a reflective minor body, it is important to investigate local stability and control at the equilibrium points. Due to high number of degrees of freedom, the Routh-Hurwitz criterion method is applied for stability analysis. Controllability is established through analysis of the local linear dynamics with solar sail pitch and precession angles as control inputs, claimed to be good choices together with trimming the sail area (Hexi and McInnes, 2005). 3

The contributions in this paper are listed as follows: (a) definition and characterization of equilibrium surfaces as a function of solar sail lightness number in a Sun-minor body system with ARP added to the solar sail EoM, (b) extension to constraints for identifying regions of where equilibrium solutions may exist, (c) definition of a region of influence (RoI) of where solar sail acceleration due to ARP is possible, and (d) extended analysis of stability and controllability at the equilibrium points. This paper is organized as follows. Section 2 describes the background to the methods used in this work and provides the solar sail dynamics in CR3BP with SRP that have previously been established in McInnes (1999); McInnes et al. (1994). Section 3 describes the means necessary to approximate ARP. In section 4 an extension to the CR3BP model is employed by adding ARP to the solar sail EoM. Modified conditions and constraints for equilibrium points are shown analytically. Results are shown in section 5 as a validation of the analytically defined framework where characterization of ARP and analysis of equilibrium points in a Sun-Vesta system are presented. Given the modified EoM, section 6 provides a general analysis of the local stability and control at the equilibrium points. Conclusions are given in section 7 which summarize the findings and discuss the theoretical impacts on the solar sail design for current and future missions. Appendix A and Appendix B contain partial derivatives of SRP and ARP acceleration components used in section 6 for stability and control analysis. 2. Background In the CR3BP, both large and minor primary objects are treated as point masses denoted as m1 and m2 , respectively, and the solar sail mass m3 is assumed to be infinitesimally small. The objects revolve around a common center of mass in the reference frame B: {ˆ x, y ˆ, ˆ z} rotating with an angular velocity ω B/I with respect to the inertial frame I: {ˆi1 , ˆi2 , ˆi3 }. By definition, the angular velocity of the rotating frame ω B/I = ˆ z is equal to the mean motion. The dimensions are normalized such that the distance between the primary masses r12 , the sum of primary masses m1 +m2 and the gravitational constant G are all taken to be unity and constant. Fig. 1 shows the geometry of the CR3BP system in the rotating reference frame, where all units are nondimensionalized. The position vectors for the solar sail with respect to

4

Figure 1: Schematic geometry of the solar sail circular restricted three-body problem in the x ˆ−y ˆ rotating frame. m1 is a luminous source, taken here to be the Sun, and m2 is the minor body that partially reflects light. The relative sizes are m1 > m2 >> m3 .

m1 and m2 in the CR3BP, as seen in Fig. 1, are defined as r1 = [(x + µ) y z]T

(1a)

r2 = [(x − (1 − µ)) y z]T

(1b)

where µ = m2 /(m1 + m2 ) is the mass ratio of the system. 2.1. Equations of motion The reduced vector form for solar sail equations of motion (EoM), with r being the vector from the center of mass to the solar sail, is defined by McInnes et al. (1994) to be ¨r + 2ω × r˙ + ∇U = a1 where the three-body pseudo-potential U may be defined as   1 2 1−µ µ 2 (x + y ) + U =− + 2 r1 r2 5

(2)

(3)

Table 1: Solar sail lightness numbers and sail loadings for relevant spacecraft.

Spacecraft β GeoStorm (McInnes, 1999) 0.02 Interstellar Heliopause Probe (McInnes, 1999) 0.3 − 0.6 0.001 IKAROS (Tsuda et al., 2013) 0.0388 − 0.0445 Sunjammer (Heiligers and McInnes, 2014)

σ (g/m2 ) 75 2.5 − 5 1500 33.71 − 38.66

The acceleration vector of a perfectly reflecting solar sail a1 due to SRP from m1 is expressed as (1 − µ) hˆr1 · n ˆ i2 n ˆ (4) a1 = β 2 r1 where the h·i term is defined as the nonnegative operator i.e.  x if x ≥ 0 hxi = 0 if x < 0 Knowing that ˆ r1 is directed along the Sun-line and since the SRP force can never be directed sunwards, the solar sail acceleration is constrained to face the Sun such that |α| ∈ (0◦ , 90◦ ) meaning only one side of the sail experiences SRP (McInnes, 1999). The solar sail lightness number β is an important design parameter, which defines the ratio of SRP acceleration to gravitational acceleration: σ∗ σ

(5)

L1 2πGm1 c

(6)

β= where σ∗ =

where σ is the sail loading m3 /A3 , c is the speed of light and L1 is the luminosity of the Sun taken to be approximately 3.846 × 1026 W . It can be seen that a large value of β is induced by lower σ = m3 /A3 , i.e. ka1 k2 is proportional to the sail area and inversely proportional to the spacecraft mass. For comparison and building intuition, Table 3 shows different values of β and corresponding σ for spacecraft that have been tested and deployed in space. It is to be noted that the examples are not perfect solar sails with 100 % specular reflection as in the theory discussed in this paper. 6

Figure 2: Schematic of solar sail cone and clock angles in the sail coordinate system.

2.2. Equilibrium solutions At equilibrium the inertial solar sail acceleration and velocity terms in Eq. (2) are zero, i.e. radiation pressure and gravitational forces are balanced such that ∇U = a1 . The sail normal vector n ˆ can thus be defined as the unit vector of the pseudo-potential gradient as shown by (McInnes et al., 1994) n ˆ=

∇U k∇U k2

(7)

With ∇U = 0 the five Lagrange points L1−5 are retrieved as the system is reduced to the classical CR3BP. The sail cone (pitch) angle α and clock (precession) angle γ define the sail attitude with respect to the coordinate system R: {ˆ r1 , ˆ r1 × ˆ z, (ˆ r1 × ˆ z) × ˆ r1 } located on the solar sail center of mass as shown in Fig. 2. The pitch angle α controls the magnitude of the radiation pressure force and the precession angle γ controls the force direction. The pitch and precession angles are   kˆ r1 × ∇U k2 α = arctan (8) ˆ r1 · ∇U 7

γ = arctan



z) × ˆ r1 ) × (ˆ r1 × ∇U )k2 k((ˆ r1 × ˆ ((ˆ r1 × ˆ z) × ˆ r1 ) · (ˆ r1 × ∇U )



(9)

Additionally, as indicated on Figs. 1 and 2, an sail orientation angle is defined with respect to r2 as   r2 k2 )(ˆ r1 · ∇U ) − (ˆ r1 · ˆ r2 )(kˆ r1 × ∇U k2 ) (kˆ r1 × ˆ (10) ψ − α = arctan (ˆ r1 · ˆ r2 )(ˆ r1 · ∇U ) + (kˆ r1 × ˆ r2 k2 )(kˆ r1 × ∇U k2 ) At an equilibrium solution, when considering no other radiation pressure forces than SRP, the solar lightness number may be expressed as β=

r12 ∇U · n ˆ (1 − µ) hˆ r1 · n ˆ i2

(11)

The classical solutions, when the solar sail is not affected by SRP, correspond to the subset β = 0. With β > 0 a particular equilibrium solution on a given surface is determined by the sail pitch and precession angles in Eq. (8) r1 · n ˆ i, a boundary and Eq. (9). By evaluating the nonnegative operator hˆ S(r1 , r2 ) = 0 can be found when ˆ r1 · n ˆ = 0 i.e. when the sail is oriented at |α| = 90◦ such that β → ∞. This boundary is expressed as S(r1 , r2 ) = x(x + µ) + y 2 −

1−µ r1 · r2 −µ 3 r1 r2

(12)

and defines two topologically disconnected surfaces that separates regions of where equilibrium solutions may and may not exist as shown in previous work (McInnes et al., 1994). The two separate forbidden volumes S1 (r1 , r2 ) and S2 (r1 , r2 ) are shown in Figs. 3a and 3b. 3. Formulation of albedo effects from minor body As an approximation to determine albedo effects, the minor body is treated as a Lambertian diffuse model with reflected light characterized by the BRDF and is commonly used in the literature (Cognion, 2013; Krag, 1974). For a Lambertian BRDF applied to the CR3BP, the reflected radiation may be expressed as a function of the phase angle between the Sun-minor body vector ˆ r12 and the viewer line of sight ˆ r2 . Since the minor body absorbs incident solar flux, black body radiation in form of emitted thermal photons would also have an effect on the solar sail, independent of whether it is located on 8

(a) x ˆ−y ˆ Plane

(b) x ˆ−ˆ z Plane Figure 3: Contour plots of forbidden regions for (a) in the rotating frame and (b) normal to the rotating frame, where is the boundary S(r1 , r2 ) = 0, 4 are the Lagrange points and ◦ is the minor body.

9

the Sun-lit or dark side of the minor body. The EoM may theoretically be extended to include black body radiation but will not be covered here. The photometric properties of planets and smaller bodies indicate the amount of light that is reflected. The albedo of the Earth, or more commonly known as ’Earthshine’, is widely varying due to the reflective properties of atmosphere, sea, ice, vegetation and land. Several techniques have been developed to estimate the albedo of the Earth and it is estimated that the global yearly average reflectance number is about ρ = 0.3 (Cognion, 2013). On the other hand, asteroids are commonly treated as gray bodies i.e. reflectance is constant for all wavelengths and has low emissivity which is consistent with previous studies (Torppa et al., 2003; Kaasalainen et al., 2002; Magnusson et al., 1996). Suitable reflectance numbers for many asteroids are estimated to be approximately ρ = 0.1−0.2 (Busch et al., 2008; Fernandez et al., 2005). Astronomers have used the visual magnitude system as technique of estimating optical signatures to determine the size and heat of distant stars. This method is applied here to characterize the albedo effects but is limited to a viewer observing only photons in the visual spectrum. Assuming the minor body is a sphere and gray body, its apparent magnitude B2 is given by   D(φ; ρ, A2 ) (13) B2 = B1 − 2.5 log10 r22 where B1 is the reference apparent magnitude, taken here to be the Sun which has the visual magnitude of about −26.78 as seen from Earth with atmospheric effects being uncorrected (Krag, 1974). Since B2 ≤ B1 then the condition D(φ; ρ, A2 )/r22 ∈ (0, 1) is naturally required. The solar flux impacting the minor body obeys the inverse square law and is proportional 2 to 1/r12 . The fraction in the logarithm in Eq. (13) is the ratio of reflected flux from the minor body experienced at the observer location to solar flux. The numerator is explicitly given by D(φ; ρ, A2 ) = ρA2 p(φ)

(14)

where A2 is the cross-sectional area of the minor body, seen as a disk from the observer, and r2 is the distance between the minor body and observer. Given that the minor body is a sphere then A2 = πd22 /4 where d2 is the diameter which is nondimensionalized in the CR3BP. The Lambertian BRDF may be expressed in terms of the body-illumination angle φ where p(φ) is the diffuse phase angle function as identified by Cognion (2013); Krag (1974). 10

The function p(φ) is defined for a diffuse Lambertian sphere as p(φ) =

2 (sin φ + (π − φ) cos φ) 3π 2

(15)

where the body-illumination angle φ and phase angle Φ are φ =π − Φ Φ = arctan

(16a) 

r2 k2 kˆ r12 × ˆ ˆ r12 · ˆ r2



(16b)

Looking at Eq. (13) and since the total solar luminosity L1 increases propor2 tionally with 1/r12 , where r12 is constant and unity in CR3BP, the apparent induced luminosity of minor body may be defined as L2 (φ; ρ, A2 ) = D(φ; ρ, A2 )L1

(17)

It is clear that the apparent induced luminosity is expressed as a fraction of the solar luminosity. 4. Equations of motion added with albedo radiation pressure Based on the methodology presented by McInnes (1999); McInnes et al. (1994) as described in section 2 and having established the basic definitions and assumptions for characterizing ARP in section 3, it is now possible to analytically incorporate ARP into the solar sail EoM to determine modified conditions for equilibrium solutions. The magnitude of the radiation pressure exerted on a surface is given by P =

W c

(18)

where c is the speed of light and W represents the energy flux of the electromagnetic radiation. Assuming the origin of light can be modeled as a point source, the energy flux varies proportionally with the inverse square of the distance traveled. Energy flux from the Sun can be written in terms of the solar luminosity L1 and is scaled by the distance r1 between Sun and solar sail as L1 (19) WSRP = 4πr12 11

Figure 4: Schematic of how SRP and ARP forces are exerted on a perfectly reflecting solar sail. Subscripts i and r stand for incident and reflected, respectively.

The reflected energy flux from the minor body is a function of the inverse square of the distance r2 . It can be written in terms of the apparent luminosity L2 (φ; ρ, A2 ) as L2 (φ; ρ, A2 ) (20) WARP = 4πr22 Starting with the force model, the total force exerted on a solar sail is expressed as X F = Fg + FSRP + FARP (21)

where Fg represents gravitational forces, FSRP is the SRP force and FARP is the ARP force. Other smaller forces are assumed to be negligible. For a specularly reflecting flat plate, the momentum transfer from incident solar flux creates a reaction force of equal magnitude from the reflected photons, generating the total SRP force FSRP . Photons due to albedo of the minor body also impact and are reflected by the surface, creating an additional force FARP in the same manner. A schematic of these radiation pressure forces exerted on the solar sail are shown in Fig. 4. The SRP and ARP forces are given as FSRP =2PSRP A3 (ˆ ui · n ˆ )2 sgn (ˆ ui · n ˆ )ˆ n FARP =2PARP A3 (ˆ vi · n ˆ )2 sgn (ˆ vi · n ˆ )ˆ n 12

(22a) (22b)

where A3 (ˆ ui · n ˆ ) and A3 (ˆ vi · n ˆ ) are the projected sail areas in u ˆ i and v ˆi directions, respectively. Extensive background on the methodology for deriving Eqs. (22a) and (22b) may be found in McInnes (1999). Geppert et al. (2011); McInnes (1999) discuss also the idea of solar sails with both of its sides being reflective. For a two-side reflective ideal or non-ideal solar sail the signum function can in general be defined as follows  if x > 0  τf 0 if x = 0 sgn(x) =  −τb if x < 0. where front and back reflectivity efficiencies τf and τb account for required sign change due to the momentum directions experienced on either side of the sail with respect to n ˆ. Using Eq. (19) and Eq. (20), then the net radiation force of Eqs. (22a) and (22b) is FSRP + FARP =

L1 A 3 (ˆ ui · n ˆ )2 sgn (ˆ ui · n ˆ )ˆ n 2 2πcr1 L2 (φ; ρ, d2 )A3 (ˆ vi · n ˆ )2 sgn (ˆ vi · n ˆ )ˆ n + 2 2πcr2

(23)

Using the definition of sail loading σ then the sum of accelerations due to SRP and ARP is X FSRP + FARP (24) a= σA3 ˆ i is For a solar sail on orbit around m1 , the direction of incident radiation u defined by the unit vector ˆ r1 . Likewise for an orbit around m2 then v ˆi is defined as ˆ r2 . Given these facts, using Eq. (5) and substituting Eq. (17) into Eq. (23), then Eq. (24) renders the net acceleration in terms of the solar sail lightness number  X 1 a =β(1 − µ) 2 (ˆ r1 · n ˆ )2 sgn (ˆ r1 · n ˆ )ˆ n r1  D(φ; ρ, A2 ) 2 + (ˆ r2 · n ˆ ) sgn (ˆ r2 · n ˆ )ˆ n (25) r22 It is clear that Eq. (25) may be explicitly decomposed into a1 and a2 . The

13

ARP acceleration, second term on the right hand side, is then explicitly given by D(φ; ρ, A2 ) a2 = β(1 − µ) (ˆ r2 · n ˆ )2 sgn (ˆ r2 · n ˆ )ˆ n (26) 2 r2 The reduced vector form for solar sail EoM in CR3BP may now be written as X ¨r + 2ω × r˙ + ∇U = a (27) which is Eq. (2) added with ARP acceleration a2 on the right hand side of the equation.

4.1. Equilibrium solutions P At equilibrium solutions it is required that ∇U = a and taking the scalar product of Eq. (27) with n ˆ , then  1 r1 · n ˆ )2 sgn (ˆ r1 · n ˆ) β(1 − µ) 2 (ˆ r1  D(φ; ρ, A2 ) 2 + (ˆ r2 · n ˆ ) sgn (ˆ r2 · n ˆ ) =∇U · n ˆ (28) r22 in which four analytical forms exist due to the signum function. The four cases describing the solar sail orientations and constraints are illustrated in ˆ is perfectly Fig. 5. If the assumption is that only the front side opposite of n reflective, i.e. τf = 1 and τb =P0, then only cases 1-3 are physically valid such that the condition ∇U · n ˆ = a·n ˆ ≥ 0 is satisfied for equilibrium solutions. Eq. (28) may therefore be rewritten as    D(φ;ρ,A2 ) 1 2 2  β(1 − µ) r2 (ˆ ˆ) + ˆ) if case 1; r1 · n (ˆ r2 · n  r22  1     β(1 − µ) 12 (ˆ X r1 · n ˆ )2 if case 2; r1 a·n ˆ= (29)  D(φ;ρ,A2 ) 2  (ˆ r β(1 − µ) · n ˆ ) if case 3;  2 r22     0 if case 4.

For any specific value of β the solar sail may direct the normal vector n ˆ ◦ ◦ unconstrainedly at |α| ∈ (0 , 180 ) and achieve equilibrium, but not when against both SRP and ARP at the same time as shown by case 4 in Fig. 5. It is clear that case 2 and 3 exclusively represent equilibrium solutions with SRP and ARP only. Rearranging Eq. (28) and assuming that τf = 1 and 14

Figure 5: Schematic geometry P of the solar sail and four orientation cases where case 1: ˆ r1 · n ˆ ≥ 0, ˆ r2 · n ˆ ≥ 0 and a·n ˆ ≥ 0; case 2: ˆ rP ˆ ≥ 0, ˆ r2 · n ˆ < 0; case 3: ˆ r1 · n ˆ < 0, 1 ·n ˆ r2 · n ˆ ≥ 0; and case 4: ˆ r1 · n ˆ < 0, ˆ r2 · n ˆ < 0 and a·n ˆ < 0.

τb = 0, then the equilibrium solutions may be represented in closed form by β as β=

∇U · n ˆ

r1 · n (1 − µ)h r12 (ˆ ˆ )2 sgn (ˆ r1 · n ˆ) + 1

D(φ;ρ,A2 ) (ˆ r2 r22

·n ˆ )2 sgn (ˆ r2 · n ˆ )i

(30)

where the h·i is a nonnegative operator. Eq. (30) is referred to as contribution (a) in section 1 and serves as an extension to previous work in McInnes (1999); McInnes et al. (1994). Locations of equilibrium points are now also a function of albedo, phase angle, minor body cross-sectional area, sail orientation and distance to the minor body, which all characterize the magnitude of ARP exerted on the solar sail. Values for β, α and γ can be evaluated at any specific equilibrium point in the three-dimensional space or equivalently equilibrium surfaces may be rendered for any chosen values of β. If both sides of the solar sail are reflective then both the numerator and denominator would be unconstrained. In this case SRP dominates in almost all regions in general given that ka1 k2 >> ka2 k2 since WSRP >> WARP . 15

4.2. Constraints on the existence of equilibrium solutions Assuming τf = 1 and τb = 0 and evaluating the gradient of the pseudopotential U as well as taking the scalar P product of Eq. (30) with n ˆ , it is required that β ≥ 0 since ∇U · n ˆ = a·n ˆP≥ 0. In the case when β → ∞, the boundary can be derived in terms of P a which corresponds to the denominator in Eq. (30). With the condition a·n ˆ = 0, the boundary V (r1 , r2 ) = 0 is derived to be V (r1 , r2 ) = r22 (ˆ r1 · n ˆ )2 sgn (ˆ r1 · n ˆ ) + r12 D(φ; ρ, A2 )(ˆ r2 · n ˆ )2 sgn (ˆ r2 · n ˆ ) (31) This boundary shows where the equilibrium solutions may and may not exist when the sail is impacted by both SRP and ARP forces. The forbidden region is shown in Figs. 6a and 6b and it is clear that the region between L1 and m2 may now have equilibrium solutions with ARP added to the EoM. The result is referred to contribution (b) in section 1 and means that the solar sail may now orient itself unreservedly such that |α| ∈ (0◦ , 180◦ ) and |(ψ − α)| ∈ (0◦ , 180◦ ) to achieve equilibrium. It can be inferred that if D(φ; ρ, A2 ) = 0, i.e. when there are no albedo effects, then V (r1 , r2 ) = S(r1 , r2 ) which is the boundary for the SRP model represented by Eq. (12). 4.3. Photometric region of influence Assuming again that τf = 1 and τb = 0, and since it is now required that SRP and ARP forces balance gravity to achieve equilibrium, it is important to investigate the scale of influence ARP has on the locations of equilibrium solutions. The equation for balanced magnitudes of SRP and ARP accelerations ka1 k2 = ka2 k2 define a boundary, similar to the mathematical definition for sphere of influence (Battin, 1999). Knowing that the normal vector has the property kˆ nk2 = 1 and dividing the magnitudes of Eq. (4) and Eq. (26) by the common term β(1 − µ), then the boundary ka1 k2 − ka2 k2 = 0 may be explicitly defined as follows r1 · n ˆ )2 sgn (ˆ r1 · n ˆ ) − r12 D(φ; ρ, A2 )(ˆ r2 · n ˆ )2 sgn (ˆ r2 · n ˆ) = 0 r22 (ˆ

(32)

This equation may also be represented by the distance from solar sail to minor body r2 in terms of explicit angles rRoI ≡ r2 = r1 D(φ; ρ, A2 )1/2 16

cos(ψ − α) cos α

(33)

(a) x ˆ−y ˆ Plane

(b) x ˆ−ˆ z Plane Figure 6: Plots of RoI volume and forbidden region for (a) in the rotating frame and (b) normal to the rotating frame, where is the boundary ka1 k2 = ka2 k2 , is the boundary V (r1 , r2 ) = 0, 4 are the Lagrange points and ◦ is the minor body.

17

Figure 7: Plot of RoI volume gap at approximately φ ≈ 90◦ in the rotating frame where is the boundary ka1 k2 = ka2 k2 , and is the non-existing SRP model boundary S(r1 , r2 ) = 0 for comparison.

which is the distance from m2 to the RoI boundary. This boundary is shown in and normal to the rotating reference frame in Figs. 6a and 6b, respectively. The RoI exists interior to the boundary ka1 k2 − ka2 k2 = 0 and defines where solar sail acceleration is affected by ARP with its scale determined by the extent of albedo effects, given that D(φ; ρ, A2 ) > 0. Exterior to the RoI boundary, the locations of equilibria are only affected by SRP. The result is referred to as contribution (c) in section 1. Principally, the RoI boundary is nearly identical to S(r1 , r2 ) except for the gaps between two surfaces created by ka1 k2 − ka2 k2 = 0 at φ ≈ ±90◦ as P shown in Fig. 7. This space grows with increasing D(φ; ρ, A2 ), i.e. when a has more contribution from ARP. In general if the back-side of the sail is reflective, i.e. τb 6= 0, then the RoI is reduced to a smaller volume since SRP dominates in almost all regions except in these gaps.

18

Table 2: Examples of objects in the Solar System (rounded up to 2 significant figures). Object Mercury Venus Earth Bennu 4179 Toutatis Eros Mars Vesta Ceres Jupiter

r12 (km) 5.79 × 107 1.08 × 108 1.50 × 108 1.34 × 108 1.40 × 108 1.69 × 108 2.28 × 108 3.53 × 108 4.14 × 108 7.41 × 108

ρA2 (km2 ) 1.27 × 106 6.90 × 107 3.84 × 107 6.47 × 10−3 6.38 × 10−1 1.96 × 101 9.06 × 106 8.67 × 104 6.33 × 104 6.14 × 109

µ 1.65 × 10−7 2.45 × 10−6 3.00 × 10−6 3.90 × 10−20 2.51 × 10−17 3.36 × 10−15 3.23 × 10−7 1.30 × 10−10 4.50 × 10−10 9.54 × 10−4

(PARP /PSRP )|L1 1.74 × 10−5 4.44 × 10−5 1.13 × 10−5 4.27 × 10−6 5.25 × 10−6 4.25 × 10−6 5.11 × 10−6 3.75 × 10−6 8.70 × 10−7 1.42 × 10−6

5. Results for a solar sail in vicinity of a minor Lambertian body This section covers simulation results based on the theory established in sections 2, 3 and 4. Before analyzing results for equilibrium points in CR3BP, it is first essential to characterize the magnitude of ARP in vicinity of the minor body with Lambertian BRDF characteristics. The chosen volume to investigate is between L1 and L2 . 5.1. Characterization of albedo radiation pressure The effect ARP has on solar sail acceleration depends on Eq. (14), i.e. the minor body distance relative to the solar sail r2 , mass m2 , albedo ρ, crosssectional area A2 , body-illumination angle φ, and distance between the Sun and the minor body r12 . To illustrate the ARP effect, Table 2 shows parameters for various bodies in the Solar System and their respective magnitude of PARP /PSRP at Lagrange point L1 , roughly assuming all orbits are circular and that objects are Lambertian. Due to its dense and highly reflective atmosphere, it is clear that Venus has the greatest ARP contribution from the bodies listed. In reality planets with atmospheres or varying landscapes, such as Earth, Venus and Jupiter, have albedos that are difficult to estimate and would need a better developed practical model to characterize reflective flux rather than using the optical signature method presented in section 3. For the following results the geometrical and physical parameters of the minor body are chosen to represent Vesta as an example, a large asteroid with m2 = 2.6 × 1020 kg, ρ = 0.35, d2 = 525.4 km. Its apparent magnitude M2 as seen by a solar sail is shown in Fig. 8. Comparative plots of PSRP and PARP are given in Figs. 9a and 9b, respectively, where it is evident that 19

Figure 8: Contour plot of apparent magnitude M2 in the rotating frame where are the apparent magnitude contours, 4 are the Lagrange points and ◦ is the minor body.

PSRP > PARP . The results are not restricted to exclude other bodies in the Solar System with discernible albedo such as for instance Jupiter, Venus, Earth, Mercury and Mars listed in Table 2 given that they are approximated as Lambertian bodies. 5.2. Equilibrium solutions in Sun-Vesta system Since Vesta is an adequate example for a Lambertian gray body, the equilibrium solutions for a Sun-Vesta system are presented here. Figs. 10a and 10b show the equilibrium solutions in terms of β contours in the rotating frame and normal to the rotating frame, respectively. The highlighted contours are given in Table 3 with corresponding sail loading σ. It can be seen that additional equilibrium solutions appear interior to the RoI boundary ka1 k2 − ka2 k2 = 0 for very high values of β or equivalently the solar sail must have a very low sail loading in order to be in tension between gravity and ARP. It is also demonstrated that the β contours exterior to the RoI are identical to the model with SRP only as defined by Eq. (11) when the sail is oriented away from ARP. The shaded region of V (r1 , rP 2 ) defines a forbidden volume for equilibrium points due to the constraint a·n ˆ ≥ 0, and as β → ∞ the equilibrium points approach m2 and the boundary V (r1 , r2 ). 20

(a) PSRP

(b) PARP Figure 9: Contour plots of (a) PSRP and (b) PARP in the rotating frame, where are the radiation pressure contours, 4 are the Lagrange points and ◦ is the minor body.

21

(a) x ˆ−y ˆ plane

(b) x ˆ−ˆ z plane Figure 10: Plots of artificial equilibrium solutions for SRP and ARP by solar sail lightness number in vicinity of the minor body for (a) in the rotating frame and (b) normal to the rotating frame, where are the β contours, is the boundary V (r1 , r2 ) = 0, 4 are the Lagrange points and ◦ is the minor body.

22

Table 3: Solar sail lightness number and respective sail loading for equilibrium solutions in vicinity of the minor body due to SRP and ARP. Contour β σ (g/m2 )

1 0.0005 3075

2 0.0008 1922

3 0.005 307.5

4 0.05 30.75

5 1 1.54

6 1000 0.0016

7 1500 0.0010

8 2000 0.00079

9 5000 0.00031

Figure 11: Three-dimensional equilibrium surfaces with solar sail lightness number β = 1500 where 4 are the Lagrange points and ◦ is the minor body.

It can be seen that the previous forbidden S1 (r1 , r2 ) volume has disappeared, and the volume size of S2 (r1 , r2 ) has decreased due to the contribution from the reflected flux scattered by BRDF. Fig. 11 shows the three-dimensional surface of equilibrium solutions for a chosen value of β = 1500. The darker surface inscribed in the lighter one is created due to ARP from the minor body and is valid for |α| ∈ (90◦ , 180◦ ), thus a solar sail with a particular β may orient itself in this range to achieve tension between ARP and gravitational forces. The particular value of β is equivalent to a sail loading of σ = 0.0010 which indicates a requirement for very low mass-to-area ratio. Fig. 12 explicitly shows a zoomed-in section of the equilibrium points for a constant β = 1500 with |α| ∈ (0◦ , 180◦ ) and |γ| = 90◦ . It can be explicitly seen that altering attitude alone changes the 23

Figure 12: Various sail pitch angles α required at equilibrium points for β = 1500 and γ = 90◦ with SRP and ARP in the rotating frame where is the β contour, × are equilibrium points, 4 is L1 Lagrange point and ◦ is the minor body.

locations of equilibrium solutions and as |α| → 180◦ then these approach m2 . 6. Stability and and control at equilibrium points Motion of a spacecraft in CR3BP is highly nonlinear and complex. The behavior of the system at artificial equilibrium points of interest is analyzed by investigating local linear dynamics. This section intentionally presents the methodology and notation provided by McInnes (1999); McInnes et al. (1994) in order to understand the similarities and differences between the extended model in this paper and those presented previous studies. 6.1. Linearization The equilibrium points, or otherwise called libration points, have the coordinates (xLi , yLi , zLi ) in the rotating reference frame B with Li being the artificial equilibrium point of interest. With only first order perturbation terms considered in the Taylor Series Expansion, the perturbed dynamics

24

around an equilibrium point with the coordinates (δx, δy, δz) may be written as rLi + δr = [(xLi + δx) (yLi + δy) (zLi + δz)]T (34) The local solar sail EoM with respect to the equilibrium points at rLi (denoted with 0∗0 ) are expressed in scalar form as ∗ ∗ + a∗1,xx + a∗2,xx )δx + (−Uxy + a∗1,xy + a∗2,xy )δy δ¨ x − 2δ y˙ =(−Uxx ∗ + (−Uxz + a∗1,xz + a∗2,xz )δz

(35a)

∗ ∗ δ y¨ + 2δ x˙ =(−Uyx + a∗1,yx + a∗2,yx )δx + (−Uyy + a∗1,yy + a∗2,yy )δy ∗ + (−Uyz + a∗1,yz + a∗2,yz )δz

δ¨ z

∗ =(−Uzx

+

+ a∗1,zx + a∗2,zx )δx + ∗ (−Uzz + a∗1,zz + a∗2,zz )δz

(35b) ∗ (−Uzy

+

a∗1,zy

+

a∗2,zy )δy (35c)

∗ where Ujk = ∂ 2 U/∂j∂k and Ujk = Ujk |Li , (j, k) ∈ (x, y, z). The partial derivatives of the acceleration terms are al,jk = ∂al,j /∂k which are, with respect to position j, the derivatives of the k component of the acceleration vector al due to radiation pressure from body ml , l ∈ (1, 2), and furthermore a∗l,jk = al,jk |Li . The partial derivatives of the accelerations a1 and a2 with respect to position are given in Appendix A. The perturbed dynamics in state space form are

x ˜˙ = A˜ x

(36)

where the state vector is x ˜ = [δx δy δz δ x˙ δ y˙ δ z] ˙ T and the state matrix A ∈ R6×6 is   I3×3 03×3 (37) A= T N where   −Uxx + a1,xx + a2,xx −Uxy + a1,xy + a2,xy −Uxz + a1,xz + a2,xz T =  −Uyx + a1,yx + a2,yx −Uyy + a1,yy + a2,yy −Uyz + a1,yz + a2,yz  (38) −Uzx + a1,zx + a2,zx −Uzy + a1,zy + a2,zy −Uzz + a1,zz + a2,zz



0  N = −2 0

25

2 0 0

 0 0 0

(39)

6.2. Stability of the local linear system A nontrivial solution is rendered by taking det(A) = 0 and gives the characteristic polynomial P(λ) =

6 X

qj λ6−j

(40)

j=0

The coefficients of the polynomial P(λ) are the following ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ (T22 T33 − T23 T32 ) − T12 (T33 T21 − T23 T31 ) q6 = T11 ∗ ∗ ∗ ∗ ∗ − T13 (T22 T31 − T21 T32 ) ∗ ∗ ∗ ∗ ∗ ∗ ∗ q5 = 2T33 (T21 T12 ) + 2(T32 T13 − T23 T31 ) ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ q4 = T11 T22 + T11 T33 + T23 T32 − T13 T31 − T12 T21 + 4T33 ∗ ∗ − T12 ) q3 = 2(T21 ∗ ∗ ∗ q2 = T11 + T22 + T33 +4 q1 = 0 q0 = 1

(41a) (41b) (41c) (41d) (41e) (41f) (41g)

To investigate whether the system is asymptotically stable according to the Routh-Hurwitz criterion, then for a nth-degree polynomial P(λ), all coefficients qi must exist (qi 6= 0) and be positive qi > 0. If there is any sign change in the Routh Array then it means the system is unstable (Franklin et al., 2014). Looking at Eqs. (41a)-(41g) implies that at least one eigenvalue will not lie in open left half-plane of the root locus diagram since q1 = 0. Thus the system is not asymptotically stable. To investigate Lyapunov stability √ by substituting for purely imaginary eigenvalues, such that λ = ικ (ι = −1), the characteristic polynomial becomes as shown in McInnes et al. (1994) P(ικ) = −κ6 + q2 κ4 − ιq3 κ3 − q4 κ2 + ιq5 κ + q6

(42)

For the condition P(λ) = 0 to be satisfied then both real and purely imaginary parts of the polynomial must be zero, thus κ6 + q2 κ4 − q4 κ2 + q6 = 0 ικ(q5 − κ2 q3 ) = 0

26

(43a) (43b)

2 Six solutions appear from this q set of equations with κi > 0, i = (1, ..., 6). With κ1 = 0 and κ2,3 = ± qq53 , the solution κ is not a consistent solution.

Consistency in Eq. (47) may be achieved by satisfying q3 = q5 = 0 where eigenvalues in Eq. (43a) are represented in conjugate pairs which may or may not have real solutions. By necessity, in order to have Lyapunov stability then ∗ ∗ q3 = 0 ⇒ (T21 − T12 ) = 0. Since potential is conservative Uxy − Uyx = 0 such that q3 = 0 ⇒ (a1,yx + a2,yx − a1,xy − a2,xy ) = 0 and with q5 = 0, it is also required that (a1,zx +a2,zx −a1,xz −a2,xz ) = 0 and (a1,yz +a2,yz −a1,zy −a2,zy ) = 0. This implies that β = 0 or ⇒ ∇ × a1 = −∇ × a2

(44)

Similar to the CR3BP with SRP acceleration a1 as indicated by McInnes et al. (1994), this implies both a1 and a2 must either be conservative and one of the accelerations must be balanced by opposite curl. Requirements for Lyapunov stability are therefore that SRP and ARP accelerations are zero such that β = 0, or the radiation pressure accelerations are conservative when n ˆkˆ r1 and n ˆkˆ r2 , i.e. |α| = 0◦ and |ψ| = 180◦ or |α| = 180◦ and |ψ| = 0◦ . These findings for instability and necessary conditions for Lyapunov stability are part of contribution (d) in section 1. 6.3. Controllability of the local linear system Since the equilibrium points are in general unstable, the sail controllability of position and velocity will be investigated using the solar sail orientation angles as control inputs which are defined as u∗ + δu = [(αLi + δα) (γLi + δγ)]T

(45)

where αLi and γLi are the nominal sail angles at the equilibrium point of interest. Thus the state space form becomes x ˜˙ = A˜ x + Bδu where A is defined in Eq. (37) and the input matrix B ∈ R6×2 is   03×2 B= ˜ a∗

27

(46)

(47)

where

 ∗ a1,xα + a∗2,xα ∗ ˜ a = a∗1,yα + a∗2,yα a∗1,zα + a∗2,zα

 a∗1,xγ + a∗2,xγ a∗1,yγ + a∗2,yγ  a∗1,zγ + a∗2,zγ

(48)

where al,jm = ∂al,j /∂m and is, with respect to orientation m ∈ (α, γ), the derivative of the j component of the acceleration vector al , l ∈ (1, 2), and furthermore a∗l,jm = al,jm |Li . The partial derivatives of a1 and a2 with respect to orientation are given in Appendix B. The controllability matrix Γ ∈ R6×12 is expressed as following   Γ = B AB A2 B A3 B A4 B A5 B (49)

where B is given by Eq. (47) and remaining columns are  ∗  ˜ a AB = N˜ a∗   N˜ a∗ 2 A B= (T + N2 )˜ a∗   a∗ (T + N2 )˜ 3 A B= (2NT − 4N)˜ a∗   (2NT − 4N)˜ a∗ 4 A B= (T2 − 4N2 − N2 T)˜ a∗  2  a∗ (T − 4N2 − N2 T)˜ 5 A B= (16N − NT2 )˜ a∗

(50a) (50b) (50c) (50d) (50e)

It is found that rank(Γ) = 6 and the system is in general controllable. The system is uncontrollable for two special cases: when the solar sail is oriented ˆ ⊥ ˆ r2 simultaneously and secondly when a1,jm = such that n ˆ ⊥ ˆ r1 and n −a2,jm as these conditions would render B = 06×2 . First case is naturally similar to the CR3BP with SRP acceleration a1 alone, i.e. the system is uncontrollable if n ˆ ⊥ˆ r1 or ˆ r1 · n ˆ = 0 as shown in McInnes et al. (1994). For |α| = 90◦ and |ψ| = 180◦ , then the solar sail is positioned at the classical collinear Lagrange points, L1−3 , and no SRP and ARP effects are induced on the solar sail. Since the system generally is controllable, then a feedback gain K may be constructed by pole placement or optimal control such that the closed-loop system Acl = (A + BK) is asymptotically stable. The findings about controllability are referred to contribution (d) stated in section 1. 28

7. Conclusions For a solar sail in a Sun-minor body system it is demonstrated that adding albedo radiation pressure (ARP) to solar radiation pressure (SRP) in the circular restricted three-body problem renders additional artificial equilibrium points that are possible for designs with low sail loading. New equilibrium points are found to exist between the Lagrange point L1 and the minor body, in a previously determined forbidden volume when SRP is only accounted for. The extensions to the equations of motion theoretically motivate novel opportunities for solar sails given that a reflective minor body has high enough albedo. Hovering locations for closer observations are possible in vicinity of both L1 and L2 points which may be considered in mission operations for flexible observations, continuous surveillance, sampling and dropping off payload. In case of a real flying solar sail it would be worth investigating how a solar sail may expand its sail area necessary to capture the albedo effects to achieve tension between ARP and gravity. The question of practical feasibility lies in a solar sail design with low sail loading, i.e. high area-to-mass ratio would be necessary to exploit the effects of ARP but may be challenging in terms of structural loads. Limitations also lie in the optical signature model used to approximate ARP as it is restricted to albedo photons from Lambertian bodies being only considered in the visual spectrum, hence alternative methods should also be investigated further. Linearization about the equilibrium points has to take ARP forces into account in order to have a well-defined framework for solar sail control in vicinity of a reflective minor body. Analyzing the local linear dynamics, it is demonstrated that the equilibrium points are in general unstable. In a special case, by necessity, Lyapunov stability may be achieved when the solar sail is oriented with its normal vector parallel to SRP and ARP. The system is found to be controllable in general for attitude angles as control inputs, thus asymptotic stability may be enforced. The theoretical findings in this paper are beneficial for analysis of potential solar sail station-keeping missions to planets, asteroids and comets, given that the bodies have appreciable albedo and the sail has very high area-to-mass ratio. Related further studies may involve black body thermal radiation, periodic orbits, transfers between equilibrium points, oblateness effects and elliptical models for a system consisting of the Sun, a reflective minor body and a solar sail. 29

References M. Macdonald, C. R. McInnes, Solar Sail Science Mission Applications and Advancement, Adv. Space Res. 48 (11) (2011) 1702–1716, doi: \let\@tempa\bibinfo@[email protected]/j.asr.2011.03.018. A. Farres, A. Jorba, Station Keeping of a Solar Sail around a Halo Orbit, Acta Atronaut. 94 (1) (2012) 527–539, doi:\let\@tempa\bibinfo@X@ doi10.1016/j.actaastro.2012.07.002. C. R. McInnes, A. J. C. McDonald, J. F. L. Simmons, Solar Sail Parking in Restricted Three-Body Systems, J. Guid. Control Dyn. 17 (2) (1994) 399–406, doi:\let\@tempa\bibinfo@[email protected]/3.21211. K. C. Howell, Three-Dimensional, Periodic, ’Halo’ Orbits, Celestial Mech. 32 (1) (1984) 53–71, doi:\let\@tempa\bibinfo@[email protected]/ BF01358403. R. W. Farquhar, A. A. Kamel, Quasi-Periodic Orbits About the Translunar Libration Point, Celestial Mech. 7 (4) (1973) 458–473, doi:\let\@tempa\ bibinfo@[email protected]/BF01227511. S. Gong, J. Li, Equilibria near asteroids for solar sail with reflection control devices, Astrophys. Space Sci. 355 (2) (2015) 213–223, doi:\let\@tempa\ bibinfo@[email protected]/s10509-014-2165-7. J. Mu, S. Gong, J. Li, Coupled Control of Reflectivity Modulated Solar Sail for GeoSail Formation Flying, J. Guid. Control Dyn. 38 (4) (2015) 740–751, doi:\let\@tempa\bibinfo@[email protected]/1.G000117. G. Aliasi, G. Mengali, A. A. Quarta, Artificial Lagrange Points for Solar Sail with Electrochromic Material Panels, J. Guid. Control Dyn. 36 (5) (2013) 1544–1550, doi:\let\@tempa\bibinfo@[email protected]/1.58167. J. A. Burns, P. L. Lamy, S. Soter, Radiation Forces on Small Particles in the Solar System, Icarus 40 (1) (1979) 1–48, doi:\let\@tempa\bibinfo@ [email protected]/0019-1035(79)90050-2. D. Nesvorny, W. F. Bottke, H. F. Levison, et al., Recent Origin of the Solar System Dust Bands, Astrophys. J. 591 (1) (2003) 486–497, doi: \let\@tempa\bibinfo@[email protected]/374807. 30

M. J. Sonter, The Technical and Economic Feasibility of Mining the Near-Earth Asteroids, Acta Atronaut. 41 (4-10) (1997) 637–647, doi: \let\@tempa\bibinfo@[email protected]/S0094-5765(98)00087-3. E. Hand, Dawn probe to look for a habitable ocean on Ceres, Science 347 (6224) (2015) 813–814, doi:\let\@tempa\bibinfo@[email protected]/ science.347.6224.813. M. G. G. T. Taylor, C. Alexander, N. Altobelli, et al., Rosetta begins its comet tale, Science 347 (6220) (2015) 387–387, doi:\let\@tempa\bibinfo@ [email protected]/science.aaa4542. C. T. Russell, C. A. Raymond, A. Coradini, et al., Dawn at Vesta: Testing the Protoplanetary Paradigm, Science 336 (6082) (2012) 684–686, doi:\let\ @tempa\bibinfo@[email protected]/science.1219381. A. Fujiwara, J. Kawaguch, D. K. Yeomans, et al., The Rubble-Pile Asteroid Itokawa as Observed by Hayabusa, Science 312 (5778) (2006) 1330–1334, doi:\let\@tempa\bibinfo@[email protected]/science.1125841. M. D. Rayman, T. C. Fraschetti, C. A. Raymond, et al., Dawn: A mission in development for exploration of main belt asteroids Vesta and Ceres, Acta Atronaut. 58 (11) (2006) 605–616, doi:\let\@tempa\bibinfo@ [email protected]/j.actaastro.2006.01.014. Shankaran, J. P. Sharma, B. Ishwar, Out-of-plane equilibrium points and stability in the generalised photogravitational restricted three body problem, Astrophys. Space Sci. 332 (1) (2011) 115–119, doi:\let\@tempa\bibinfo@ [email protected]/s10509-010-0494-8. J. Torppa, M. Kaasalainen, T. Michalowski, et al., Shapes and rotational properties of thirty asteroids from photometric data, Icarus 164 (2) (2003) 346–383, doi:\let\@tempa\bibinfo@[email protected]/S0019-1035(03) 00146-5. M. Kaasalainen, J. Torppa, J. Piironen, Models of Twenty Asteroids from Photometric Data, Icarus 159 (2) (2002) 369–395, doi:\let\@tempa\ bibinfo@[email protected]/icar.2002.6907.

31

P. Magnusson, M. Dahlgren, M. A. Barucci, et al., Photometric Observations and Modeling of Asteroid 1620 Geographos, Icarus 123 (1) (1996) 227–244, doi:\let\@tempa\bibinfo@[email protected]/icar.1996.0151. S. B. Broschart, D. J. Scheeres, Control of Hovering Spacecraft Near Small Bodies: Application to Asteroid 25143 Itokawa, J. Guid. Control Dyn. 28 (2) (2005) 343–354, doi:\let\@tempa\bibinfo@[email protected]/1.3890. E. Morrow, D. J. Scheeres, D. Lubin, Solar Sail Orbit Operations at Asteroids, J. Spacecr. Rockets 38 (2) (2001) 279–286, doi:\let\@tempa\bibinfo@ [email protected]/2.3682. D. J. Scheeres, F. Marzari, Spacecraft Dynamics in the Vicinity of a Comet, J. Astronaut. Sci. 50(1) (2002) 35–52. B. Hexi, C. R. McInnes, Solar Sail Orbits at Artificial Sun-Earth Libration Points, J. Guid. Control Dyn. 28 (6) (2005) 1328–1331, doi:\let\@tempa\ bibinfo@[email protected]/1.14598. C. R. McInnes, Solar Sailing: Technology, Dynamics and Mission Applications, Springer Praxis, 1999. Y. Tsuda, O. Mori, R. Funase, et al., Achievement of IKAROS Japanese Deep Space Solar Sail Demonstration Mission, Acta Atronaut. 82 (2) (2013) 183–188, doi:\let\@tempa\bibinfo@[email protected]/j. actaastro.2012.03.032. J. Heiligers, C. R. McInnes, Novel Solar Sail Mission Concepts for Space Weather Forecasting, in: 24th AAS/AIAA Space Flight Mechanics Meeting 2014, AAS 14-239, New Mexico, 2014. R. Cognion, Observations and modeling of GEO satellites at large phase angles, in: Advanced Maui Optical and Space Surveillance Technologies Conference, Wailea, Maui, Hawaii, 2013. W. E. Krag, Visible Magnitude of Typical Satellites in Synchronous Orbits, Tech. Rep. ESD-TR-74-278, Massachusetts Institute of Technology, 1974. M. W. Busch, L. A. M. Benner, S. J. Ostro, Physical properties of a NearEarth Asteroid (33342) 1998 WT24, Icarus 195 (2) (2008) 614–621, doi: \let\@tempa\bibinfo@[email protected]/j.icarus.2008.01.020. 32

Y. R. Fernandez, D. C. Jewitt, S. S. Sheppard, Albedos of Asteroids in Comet-Like Orbits, Astron. J. 130 (1) (2005) 308–318, doi:\let\@tempa\ bibinfo@[email protected]/430802. U. Geppert, B. Biering, F. Lura, et al., The 3-step DLR-ESA Gossamer road to solar sailing, Adv. Space Res. 48 (11) (2011) 1695–1701, doi:\let\ @tempa\bibinfo@[email protected]/j.asr.2010.09.016. R. H. Battin, An Introduction to the Mathematics and Methods of Astrodynamics (Revised Edition), American Institute of Aeronautics and Astronautics, 1999. G. F. Franklin, J. D. Powell, A. Emami-Naeini, Feedback Control of Dynamic Systems, Prentice Hall, Upper Saddle River, NJ, 7th edn., 2014.

33

Appendix A. Partial derivatives of the solar sail accelerations relative to position The scalar equations for the partial derivatives of the solar sail accelerations a1 and a2 with respect to position given in section 6, assumed that solar sail cone and clock angles (α and γ) are independent of position, are expressed here.  (1 − µ) 2(x + µ)nx cos α(y 2 + z 2 ) 2 a1,xx =β cos α − + r12 r12 r13  sin α sin γy(x + µ) sin α cos γz((x + µ)4 − y 2 z 2 − y 4 ) − + (A.1) kr1 × ˆ zk32 k(r1 × ˆ z) × r1 k32  (1 − µ) 2ynx cos α(x + µ)y sin α sin γ(x + µ)2 2 a1,xy =β cos α − − + r12 r12 r13 kr1 × ˆ zk32  sin α cos γ(x + µ)yz(2(x + µ)2 + 2y 2 + z 2 ) + (A.2) k(r1 × ˆ z) × r1 k32  (1 − µ) 2znx cos α(x + µ)z 2 a1,xz =β cos α − 2 − 2 r1 r1 r13  sin α cos γ(x + µ)((x + µ)2 + y 2 )2 − (A.3) k(r1 × ˆ z) × r1 k32  (1 − µ) 2(x + µ)ny cos α(x + µ)y sin α sin γy 2 2 cos α − − − a1,yx =β r12 r12 r13 kr1 × ˆ zk32  sin α cos γ(x + µ)yz(2(x + µ)2 + 2y 2 + z 2 ) (A.4) + k(r1 × ˆ z) × r1 k32  (1 − µ) 2yny cos α((x + µ)2 + z 2 ) 2 a1,yy =β cos α − + r12 r12 r13  sin α sin γy 2 sin α cos γz((x + µ)2 (z 2 + (x + µ)2 ) − y 4 ) + − (A.5) kr1 × ˆ zk32 k(r1 × ˆ z) × r1 k32

34

a1,yz

a1,zx

a1,zy

a1,zz

a2,xx

a2,xy

a2,xz

 (1 − µ) 2zny cos αzy 2 =β cos α − 2 − 2 r1 r1 r13  sin α cos γy((x + µ)2 + y 2 )2 − (A.6) k(r1 × ˆ z) × r1 k32  (1 − µ) 2(x + µ)nz cos α(x + µ)z =β cos2 α − − 2 r1 r12 r13  sin α cos γ(x + µ)z 2 ((x + µ)2 + y 2 ) + (A.7) k(r1 × ˆ z) × r1 k32  (1 − µ) 2ynz cos αyz 2 =β cos α − 2 − 2 r1 r1 r13  sin α cos γyz 2 ((x + µ)2 + y 2 ) + (A.8) k(r1 × ˆ z) × r1 k32  2znz cos α((x + µ)2 + y 2 ) (1 − µ) 2 =β cos α − 2 + r12 r1 r13  sin α cos γz((x + µ)2 + y 2 )2 − (A.9) k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ)r12 2(x − 1 + µ)nx 2 =β cos (ψ − α) − 2 r2 r22 cos α(y 2 + z 2 ) sin α sin γy(x + µ) + − r13 kr1 × ˆ zk32  sin α cos γz((x + µ)4 − y 2 z 2 − y 4 ) + (A.10) k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ) 2ynx cos α(x + µ)y 2 =β cos (ψ − α) − − r22 r22 r13 sin α sin γ(x + µ)2 + kr1 × ˆ zk32  sin α cos γ(x + µ)yz(2(x + µ)2 + 2y 2 + z 2 ) + (A.11) k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ) 2znx cos α(x + µ)z 2 =β cos (ψ − α) − 2 − 2 r2 r2 r13  sin α cos γ(x + µ)((x + µ)2 + y 2 )2 − (A.12) k(r1 × ˆ z) × r1 k32

35

a2,yx

a2,yy

a2,yz

a2,zx

a2,zy

a2,zz

 (1 − µ)ρA2 p(φ) 2(x − 1 + µ)ny 2 =β cos (ψ − α) − 2 r2 r22 cos α(x + µ)y sin α sin γy 2 − − r13 kr1 × ˆ zk32  sin α cos γ(x + µ)yz(2(x + µ)2 + 2y 2 + z 2 ) + k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ) 2yny 2 =β cos (ψ − α) − 2 2 r2 r2 2 2 cos α((x + µ) + z ) sin α sin γy 2 + + r13 kr1 × ˆ zk32  sin α cos γz((x + µ)2 (z 2 + (x + µ)2 ) − y 4 ) − k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ) 2zny cos αzy 2 =β cos (ψ − α) − 2 − 4 r2 r2 r13  sin α cos γy((x + µ)2 + y 2 )2 − k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ) 2(x − 1 + µ)nz 2 =β cos (ψ − α) − 2 r2 r22  cos α(x + µ)z sin α cos γ(x + µ)z 2 ((x + µ)2 + y 2 ) − + r13 k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ) 2ynz cos αyz =β cos2 (ψ − α) − 2 − 2 r2 r2 r13  sin α cos γyz 2 ((x + µ)2 + y 2 ) + k(r1 × ˆ z) × r1 k32  (1 − µ)ρA2 p(φ) 2znz 2 =β cos (ψ − α) − 2 2 r2 r2  2 2 cos α((x + µ) + y ) sin α cos γz((x + µ)2 + y 2 )2 + − r13 k(r1 × ˆ z) × r1 k32

36

(A.13)

(A.14)

(A.15)

(A.16)

(A.17)

(A.18)

Appendix B. Partial derivatives of the solar sail accelerations relative to orientation The scalar equations for the partial derivatives of the solar sail accelerations a1 and a2 with respect to orientation given in section 6 are expressed here.  (1 − µ) 3 cos α sin α(x + µ) (1 − 3 sin2 α) sin γy a1,xα =β cos α − + r12 r1 kr1 × ˆ zk2  2 (1 − 3 sin α) cos γ(x + µ)z (B.1) − k(r1 × ˆ z) × r1 k2  (1 − µ) 3 cos α sin αy (1 − 3 sin2 α) sin γ(x + µ) a1,yα = β cos α − − r12 r1 kr1 × ˆ zk2  2 (1 − 3 sin α) cos γyz − (B.2) k(r1 × ˆ z) × r1 k2  (1 − µ) 3 cos α sin αz a1,zα = β cos α − 2 r1 r1  2 2 (1 − 3 sin α) cos γ(y + (x + µ)2 ) − (B.3) k(r1 × ˆ z) × r1 k2   sin α cos γy sin α sin γ(x + µ)z (1 − µ) 2 cos α + (B.4) a1,xγ = β r12 kr1 × ˆ zk2 k(r1 × ˆ z) × r1 k2   (1 − µ) sin α cos γ(x + µ) sin α sin γyz 2 a1,yγ = β cos α − + (B.5) r12 kr1 × ˆ zk2 k(r1 × ˆ z) × r1 k2   (1 − µ) sin α sin γ(y 2 + (x + µ)2 ) 2 cos α − (B.6) a1,zγ = β r12 k(r1 × ˆ z) × r1 k2 (1 − µ)ρA2 p(φ) cos(ψ − α) a2,xα =β r22  (2 sin(ψ − α) cos α − cos(ψ − α) sin α) (x + µ) r1 (2 sin(ψ − α) sin α + cos(ψ − α) cos α) sin γy + kr1 × ˆ zk2  (2 sin(ψ − α) sin α + cos(ψ − α) cos α) cos γ(x + µ)z − (B.7) k(r1 × ˆ z) × r1 k2

37

(1 − µ)ρA2 p(φ) a2,yα =β cos(ψ − α) r22  (2 sin(ψ − α) cos α − cos(ψ − α) sin α) y r1 (2 sin(ψ − α) sin α + cos(ψ − α) cos α) sin γ(x + µ) − kr1 × ˆ zk2  (2 sin(ψ − α) sin α + cos(ψ − α) cos α) yz − (B.8) k(r1 × ˆ z) × r1 k2 (1 − µ)ρA2 p(φ) cos(ψ − α) a2,zα =β r22  (2 sin(ψ − α) cos α − cos(ψ − α) sin α) z r1  (2 sin(ψ − α) sin α + cos(ψ − α) cos α) cos γ(y 2 + (x + µ)2 ) − k(r1 × ˆ z) × r1 k2 (B.9)  (1 − µ)ρA2 p(φ) sin α cos γy cos2 (ψ − α) a2,xγ =β 2 r2 kr1 × ˆ zk2  sin α sin γ(x + µ)z + (B.10) k(r1 × ˆ z) × r1 k2  β(1 − µ)ρA2 p(φ) sin α cos γ(x + µ) 2 a2,yγ =β cos (ψ − α) − 2 r2 kr1 × ˆ zk2  sin α sin γyz + (B.11) k(r1 × ˆ z) × r1 k2   (1 − µ)ρp(φ)A2 sin α sin γ(y 2 + (x + µ)2 ) 2 a2,zγ =β cos (ψ − α) − r22 k(r1 × ˆ z) × r1 k2 (B.12)

38