An analytical description of three-dimensional heliocentric solar sail orbits

This paper introduces a simple analytical approximation to three-dimensional heliocentric solar sail orbits where the only forces considered are solar gravity and solar radiation. The approximation is based upon the previously studied hodograph transformation and provides a description of the inclination, longitude of ascending node and true latitude for a specific set of initial conditions. It is shown that the rotational symmetry of a heliocentric orbit allows this specific solution to be mapped onto a solution with arbitrary initial conditions. The approximation is then compared to the numerical results for a solar sail on an Earth escape trajectory with an area to mass ratio up to twice as high as current technology allows.


Introduction
Solar sails are an alternative form of spacecraft propulsion which use solar radiation pressure to push the spacecraft through space (McInnes 1999a). Conventional thrusters are replaced by a large, reflective membrane surface onto which the incoming solar photons transfer their momentum, accelerating the spacecraft. This provides a solar sail with an essentially unlimited supply of fuel and allows it to undertake missions that spacecraft with conventional propulsion would find impossible (McInnes 1999b).
Although the concept and mathematics behind solar sails was introduced in the 1950s (Garwin 1958), it is only in the past few years that we have finally seen the first solar sail launched. The successful demonstration of the IKAROS solar sail mission (Tsuda et al. 2011), launched in 2010 by the Japan Aerospace Exploration Agency (JAXA), has spurred on more mission concepts with both the Planetary Society (2014) and Surrey Space Centre (2014) set to launch solar sail technology demonstration and deorbiting missions in 2016.
While the majority of the objects in the solar system lie within a few degrees of the ecliptic plane, there are a significant number of objects that do not. These include Pluto, the Heliopause and minor bodies in the asteroid belt and Kuiper Belt, with many of these objects having an inclination of tens of degrees with respect to the ecliptic. Sending spacecraft to these objects is therefore challenging due to the large delta-v required to change a spacecraft's inclination. Gravitational assists are often used but this will place constraints on mission launch dates and can result in highly complex trajectories. Solar sails are able to solve this problem since the continuous source of photons allows them to change their inclination without any fuel. Solar sail use in these sorts of missions has been investigated numerically by Dachwald (2003Dachwald ( , 2004 using evolutionary neurocontrol optimization techniques. However despite this inclination change advantage, three-dimensional heliocentric orbits are still poorly understood from an analytical point of view. Early research focussed on the two-dimensional equations of motion and the solutions to these equations. Tsu (1959), London (1960) and Bacon (1959) all investigated the only known exact analytical solution; the logarithmic spiral solution. Bacon's work was actually concerned with continuous thrust spacecraft, not specifically solar sails, although the results were directly applicable. The work of Van Der Ha and Modi (1978) then provided one of the few analytical insights into three-dimensional heliocentric solar sail orbits. It was shown that the in-plane perturbations could be uncoupled from the orientation of the orbital plane via a choice of independent variable. Exact solutions were found for the threedimensional logarithmic spiral case with the inclination found to exhibit periodic behaviour. Series approximations were required to investigate the long term behaviour of a solar sail for non-logarithmic orbits. In terms of the orbital angles, the behaviour of the inclination and true latitude were investigated for near-circular initial orbits. The use of the independent variable here to uncouple the equations will form the basis of the approximate solution presented here.
This was followed by several papers over the following decades, most notably by Sauer (2000), Leipold (1999) and Vulpetti (1996Vulpetti ( , 1997Vulpetti ( , 1999, that used numerical methods to examine high performance solar sails and missions to the outer solar system. These studies in particular showed the enormous potential of solar sails. Missions to the outer solar system, including missions requiring significant inclination changes, could be accomplished in much shorter time scales than with conventional propulsion. Perturbation effects such as the Poynting-Robertson effect have also been studied in relation to these types of missions by Vázquez-Poritz 2011, 2013), showing decreased Heliocentric distance. Further examination of heliocentric orbits were almost entirely numerical.
Current research is heavily focussed on solar sails in the Earth-Sun or Earth-Moon circular restricted three-body problem. A large number of papers have been published on the use of a solar sail to generate and control artificial Lagrange points (McInnes et al. 1994;McInnes and Macdonald 2005;Macdonald et al. 2007;McInnes 2007, 2008a, b;Farres and Jorba 2008). While solar sails in previously reviewed studies have used large solar sails as a propulsion system for highly non-Keplerian orbits and deep space missions, current research uses much smaller solar sails for near-Earth missions. Now with larger solar sails set to be launched, focus is starting to shift back to the heliocentric solar sail orbits that dominated early research.
More recently, a geometric tool, the hodograph, has been used to analyse solar sail orbits. The hodograph is essentially a visualization of the tangential and radial velocity and it's mathematical treatment and usage in Keplerian orbits and orbital transfers were discussed by Szebehely (1964) and Battin (1999). It was first shown by Wokes et al. (2008) that the four-dimensional planar dynamics of a fixed Sun angle solar sail can be reduced down to a two-dimensional problem using the transformation to the hodograph. Trajectories could then be classified as either spiral, angular momentum reversal or hyperbolic based on the selection of initial conditions in this two-dimensional hodograph phase space. This was the first time the long term behaviour of two-dimensional heliocentric solar sail orbits had been described using an analytical model and we seek to build upon this work here by extending their approach to three-dimensional orbits. This paper will combine the hodograph results with the early work of Van Der Ha and Modi (1978) to present a full analytical description of three-dimensional heliocentric solar sail orbits. Although further work has been done on the angular momentum reversal trajectory Zeng et al. (2011), for this work the initial conditions near the logarithmic spiral orbit where chosen, as this is relevant to Earth escape orbits as their initial conditions are close to the logarithmic spiral. We extend the work of Van der Ha and Modi to include the longitude of ascending node and classify the possible long term behaviour of a fixed angle solar sail into distinct categories. Sections 2 and 3 will introduce the fixed angle solar sail and detail the hodographs use in trajectories. Section 4 derives the three-dimensional equations of motion showing how the hodograph solution couples with the orientation of the orbital plane. In Sect. 5 we demonstrate how the rotational symmetry of a heliocentric orbit allows us to simplify the equations of motion considerably by showing that any two orbits with identical hodograph parameters can be mapped to one another. We use this result in Sect. 6 to derive a complete analytical description of solar sail orbits near the logarithmic spiral initial conditions.

Force model for a fixed angle solar sail
In this section, we present a simple model for the forces on a fixed angle solar sail due to solar gravitational force and solar radiation pressure. Perturbation effects, such as the Poynting-Robertson effect as discussed in the introduction, will not be considered here. First, the definition of a fixed angle in three dimensions is provided. Letr be the radial unit vector connecting the Sun and solar sail andn the unit vector pointing normal to the sail. For a perfectly reflecting sail, a photon striking the sail will transfer double the momentum of a perfectly absorbing sail, with all the force directed alongn. The incident angle is the angle between the radial vector and sail normal vector and is denoted as α. Following McInnes (1999a), we also define a secondary angle, known as the clock angle δ which is defined as the angle between the angular momentum vector and the sail normal projected onto the plane formed byĥ andλ =r ×ĥ. Figure 1 shows these angles.
Following the derivation of McInnes (1999a), we have that the acceleration due to solar radiation pressure is given by where r is the radial distance, μ the standard gravitational parameter and β is a dimensionless parameter signifying the ratio of solar radiation force to gravitational force. We define β as  where L s is the solar luminosity at 1AU, equal to 3.839 × 10 26 W and c the speed of light. Definingλ =r ×ĥ, the sail normal vector can then be written aŝ n = cos αr + sin α sin δλ + sin α cos δĥ.
A fixed angle solar sail is then one for which the components of the sail normal vector are constant alongr,λ andĥ. Including solar gravity, the total acceleration can be written as where

The hodograph
The hodograph is a mathematical tool for looking at spacecraft velocity. In this section we provide a derivation along with an example. Following Wokes et al. (2008), the equations of motion for a fixed angle solar sail in a two-dimensional heliocentric orbit are given as where r is the radial distance, h = r 2ψ the angular momentum andṙ the radial velocity. In this two-dimensional case, ψ is just the true anomaly. This is however not true for the three-dimensional case as will be discussed later. We then use a change of variables in order to derive the hodograph reduction. Let where v and w are non-dimensional. Taking the derivative of v and w with respect to ψ (denoted as throughout this paper), we get We see then that the hodograph transformation reduces the system from four-dimensional down to two-dimensional as the derivatives of v and w only depend on v and w. The full planar dynamics can be recovered by including either the radial distance or magnitude of angular momentum, since a choice of one of these allows the other to be found through the equation for v. From Battin (1999), the eccentricity can also be written as We also normalize the variables such that μ = 1, the initial radial distance is 1 and one time unit is 2π. From Eq. (12), the differential equations for r and h along with their initial conditions are then This v − w system of equations exhibit two equilibrium points, (ṽ 1 ,w) and (ṽ 2 ,w) wherẽ These correspond to the logarithmic spiral trajectories when β = 0. When β = 0, the system reduces to Keplerian motion. In the hodograph phase space this is represented by a family of circles centred at v = 1, w = 0 with the radius being the eccentricity of the orbit. As such, a Keplerian orbit with zero eccentricity is simply the point (1, 0) in the phase space. This can be considered a realistic starting condition for a solar sail that has been released on an Earth escape orbit assuming the Earth is on a circular orbit around the Sun at 1AU. In order to develop analytical approximations for the full three-dimensional solar sail dynamics, we will first need to find approximations for v and w as their differential equations are not solvable analytically. Approximate solutions for the behaviour of v and w nearṽ 2 can be found using a linear approximation. Let v =ṽ 2 + x, w =w + y and a = k 2 2ṽ 2 where x and y are small deviations away from the equilibrium point. To first order, and neglecting any ak 2 terms, we find that where θ 0 is measured anti-clockwise from the horizontal nullcline and ρ 0 is the distance from the equilibrium pointṽ 2 , as shown in Fig. 2. Solving Eqs. (20) and (21) then yields the following solutions for v and w v =ṽ 2 + ρ 0 cos(ψ + θ 0 )e aψ (24) To complete the description of the in-plane trajectories, we will also derive the approximation for h as this is simpler than solving the equation for r . By substituting the above approximation for v into Eq. (17) and performing an expansion of 1/v in terms of ρ 0 , we find that From this we can recover r and so the complete dynamics of the two-dimensional solar sail orbit can be expressed in terms of ψ.
We show the accuracy of the v, w and h approximation in Fig. 3 by comparing the numerical solution of the equations of motion with our approximations for different values of β. We choose the initial values of v, w and h to be v 0 = 1, w 0 = 0 and h 0 = 1 to represent an Earth escape orbit. We let β range from 0.001 to 0.1 with α = arctan(1/ √ 2), δ = 1 and compare the solutions for an orbit lasting approximately 35 years. This value of alpha is chosen to maximize the tangential force (McInnes 1999b) and will be used throughout this paper. We use β = 0.1 as the maximum value here as this is twice what current technology will allow and so will exaggerate any errors. The mean error is calculated by summing the absolute value of the difference between the numerical solution and approximation at each integration step and dividing by the number of time steps, while the maximum error is simply the largest error recorded at any time step over a whole orbit. In doing this we can see that

Equations of motion
In this section we derive the full equations of motion for a solar sail in a three-dimensional heliocentric orbit. Define an inertial coordinate systemê x ,ê y ,ê z such that the spacecraft is initially situated on theê x axis and theê z axis is aligned with the initial angular momentum vector. The unit vectorsr,λ andĥ, as described in Sect. 2, form a non-uniformly rotating coordinate system with respect to this inertial system, we refer to this rotating system as the orbital frame. The orbital frame, and hence the orientation of the orbital plane, is described using the classical orbital elements; the longitude of ascending node, Ω, the inclination i and the true latitude λ = ω + ν where ω is the longitude of periapsis and ν the true anomaly. We refer to these three angles as the orbital angles. The orbital frame can be written in terms of the inertial coordinate system using three consecutive Euler rotations as follows ⎛ ⎝r λ where R is a matrix representing a 3-1-3 rotation using the angles Ω followed by i and finally λ. The angular velocity is then defined aṡ with the components given by The goal now is to find alternative expressions for these components and so compute the derivatives of the orbital angles. From the definition of angular momentum for a perturbed three-dimensional orbit as given by Brouwer and Clemence (1961), we have that and hence ω z = −h/r 2 . We will now derive the torque vector τ in two different ways. Doing this allows a comparison of components which will lead to alternative expressions for ω x and ω z . First, simply consider the definition of the torque vector which is τ = r × F with F the force on the solar sail. We find the force from the acceleration given by Eq. (5) and hence the torque is However, the torque can also be written as the derivative of the angular momentum and so we have τ =ḣ and hence τ = d dt hĥ =ḣĥ + hḣ.
Using the rotation matrix R from Eq. (27), we find the time derivative of the angular momentum unit vector in the rotating coordinate system to bė Substituting this into Eq. (32), we then compare the unit vectors components to find that ω y = 0 and Now, by combining Eqs. (28-30) with Eq. (34), ω z = −h/r 2 and ω y = 0, we can rearrange to solve for the derivatives of the orbital angles. By changing the independent variable to ψ, a factor h/r 2 can be removed from the derivatives and therefore we derive the following set of differential equations These six equations then fully describe the solar sail's three-dimensional orbit, with the inplane dynamics being described by the first three equations, and the orientation of the orbital plane and position of the solar sail in that plane described by the final three. For completeness, we use the definition of the hodograph variables from Battin (1999) to derive the argument of periapsis as The benefit of these equations is that they neatly split the kinematics and dynamics, with the orbital angles being coupled to the hodograph variables only through v. The problems with the equations are the singularities at i = 0, caused by the longitude of ascending node being undefined for a planar orbit. This causes problems when trying to find approximate solutions as λ and Ω change rapidly for small inclination values. In the next section we'll describe how we solve this problem by showing that the orbit generated by any initial set of orbital angles can be mapped to any other orbit with different orbital angles.

Rotational symmetry of the problem
There are six variables that are required to recover the full dynamics, giving us a six dimensional solution space. However, with both gravitational force and solar radiation pressure being dependent only on r and not the orientation in space, the problem admits rotational symmetry. We now use this to reduce the complexity of the solution space. Let F 1 denote an inertial frame. Consider an orbit with some arbitrary initial conditions, as measured in this frame, to be Let F 2 denote a second inertial frame which measures the initial conditions for the same orbit as with v 0 , w 0 and r 0 identical in both sets of initial conditions. Although both of these frames will measure the orbital angles differently as the orbit evolves, there is in fact only a single unique orbit here. Through rotational symmetry, it is possible to use some transformation to map the orbital angles measured by one inertial frame onto the orbital angles measured by another inertial frame. The question is then, If, at some point in time, the orbital angles, as measured in frame F 1 , are λ 1 , i 1 and Ω 1 , what are the orbital angles as measured in frame F 2 ?
Solving this will mean that we need only consider the dynamics of a single set of initial conditions for the orbital angles. Other orbits with arbitrary initial conditions can be obtained from this by applying rotations. There is a constant rotation matrix, R, defining the transformation from F 1 to F 2 which can be written as Taking the unit vectors of F 1 as coordinate axis of the identity matrix. The unit vectors of F 2 , denoted asx 2 ,ŷ 2 andẑ 2 are then the row vectors of R. The unit vectors of the orbital frame can be written as Fig. 4 Relationship between the two inertial frames and the orbital frame Figure 4 illustrates the relationship between F 1 , F 2 and the orbital frame.
We can now calculate the orbital angles, as measured in F 2 , to be the following cos λ 2 = n | n | ·r where n =ẑ 2 ×ĥ is the vector in the direction of the longitude of the ascending node. The orbital angles in F 2 can now be written in terms of only the orbital angles in F 1 . As a final reduction in complexity, we choose the initial radial distance to be unity. The three dimensional solar sail orbit is then defined only in terms of the initial hodograph variables, as any other orbit can be found through the rotation method described above. The initial condition solution space has therefore been reduced from six variables to two. These are supplemented by the sail parameters α, δ and β.

Analytic description of orbits near the log spiral
The rotational symmetry and angle mapping in Eq. (41) mean that if the equations of motion are solved for one set of initial orbital angles, they are solved for all sets. Consider an initial set of orbital angles With these initial conditions, the equations of motion given in Eq. (35) can be simplified. Assume that the variation in the inclination is small enough such that Using the approximation of the hodograph variables in Eq. (24) the differential equation for λ then becomes λ = 1 and hence, with λ(0) = 0, we find that Combining this with a first order expansion of 1/v from Eq. (24), we obtain the following differential equations for the inclination and longitude of ascending node: Solving this yields the following solutions where and To recover the full three-dimensional orbit, we can use either the radial distance or angular momentum derived earlier. The equations of motion are now solved for the given initial conditions meaning the analytic solution for any initial condition can be found through rotational symmetry. As an example to demonstrate the accuracy of the above solution, we use the same parameters as in Sect. 3 and compare the numerical solution and the approximate solution with β ranging from 0 to 0.1 (Fig. 5). We can see from Eqs. (48-51) that the long term behaviour of the orbital angles can be described simply in terms of v 0 and w 0 . The general behaviour for both angles is an oscillation around the initial condition withṽ 2 − v 0 andw − w 0 determining in which direction the drift away from this happens. For the inclination this neatly splits the hodograph phase space into two regions; v 0 <ṽ 2 and v 0 >ṽ 2 . Similarly for the longitude of ascending node the two regions are; w 0 <w and w 0 >w. The benefit of using this 90 • inclination frame is then apparent as the long term behaviour of a solar sail can be determined from one of four regions in the hodograph phase space which are found by combining the inclination and longitude of ascending node regions. As the problem admits rotational symmetry, there will also be four regions describing the behaviour of the angles in any other arbitrary frame. For the 90 • inclination frame, the behaviour is illustrated in Fig. 6. The initial conditions used to obtain these results were v 0 =ṽ 2 + 0.1 cos θ and w 0 =w + 0.1 sin θ for θ ∈ {45 • , 135 • , 225 • , 315 • } with the sail parameters being β = 0.2, α = arctan(1/ √ 2) and δ = 1.

Conclusions
This paper introduces an analytical approximation to three-dimensional solar sail orbits for trajectories starting close to the logarithmic spiral. It was shown that, by choosing a frame in which the initial inclination was 90 • , the equations of motion could be simplified considerably and accurate approximations to the orbital angles made. The rotational symmetry of the problem then allowed us to transform this approximation into any other inertial frame and hence solutions for any initial set of orbital angles could be found. The use of the 90 • inclination frame allowed the hodograph phase space to be divided into four distinct regions with each having a different long term behaviour for the inclination and longitude of ascending node. The general behaviour of these angles was shown to be an oscillation combined with a slowly increasing deviation away from the initial value, driven by the sail orientation and the distance from the logarithmic spiral initial condition.