Table of Contents Author Guidelines Submit a Manuscript
International Journal of Aerospace Engineering
Volume 2011 (2011), Article ID 476197, 13 pages
Research Article

Generating Solar Sail Trajectories in the Earth-Moon System Using Augmented Finite-Difference Methods

School of Aeronautics and Astronautics, Purdue University, 701 W. Stadium Avenue, West Lafayette, IN 47907-2045, USA

Received 8 November 2010; Accepted 1 February 2011

Academic Editor: Srinivas R. Vadali

Copyright © 2011 Geoffrey G. Wawrzyniak and Kathleen C. Howell. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Using a solar sail, a spacecraft orbit can be offset from a central body such that the orbital plane is displaced from the gravitational center. Such a trajectory might be desirable for a single-spacecraft relay to support communications with an outpost at the lunar south pole. Although trajectory design within the context of the Earth-Moon restricted problem is advantageous for this problem, it is difficult to envision the design space for offset orbits. Numerical techniques to solve boundary value problems can be employed to understand this challenging dynamical regime. Numerical finite-difference schemes are simple to understand and implement. Two augmented finite-difference methods (FDMs) are developed and compared to a Hermite-Simpson collocation scheme. With 101 evenly spaced nodes, solutions from the FDM are locally accurate to within 1740 km. Other methods, such as collocation, offer more accurate solutions, but these gains are mitigated when solutions resulting from simple models are migrated to higher-fidelity models. The primary purpose of using a simple, lower-fidelity, augmented finite-difference method is to quickly and easily generate accurate trajectories.

1. Introduction

When a permanent outpost on the Moon to support extended human expeditions is eventually established, the astronauts at the facility will require a continual communications link with the Earth. A leading candidate location for the lunar base is at the south pole, which is not always in view of either antennas on the Earth or space-based transceivers, such as the TDRSS satellites, in geosynchronous orbit. Therefore, at least two spacecraft in traditional polar orbits about the Moon are required for an Earth-Moon link [1].

In contrast to a constellation of multiple spacecraft, alternative communications strategies that rely on only one satellite do exist, using current or near-future technology. Advanced propulsion concepts, such as low-thrust ion engines, as well as solar sails, supply a force in addition to gravity and can actually offset an orbit from the Moon [24], that is, the orbit plane is displaced from the gravitational center and the spacecraft appears to hover above the surface. Such an orbit allows one spacecraft to be in view of a location on the lunar surface at all times. This dynamical system is complex due to the gravitational effects of the two primaries and, in the case of solar sails, the fact that the Sun, which influences the direction and magnitude of the resulting sail thrust vector, continually moves relative to a fixed Earth-Moon system. A previous approach, based on an understanding of the dynamical structure and using that knowledge to design an orbit, has been successful in developing some solar sail orbits in the vicinity of the lunar Lagrange points [59]. However, motion below one of the lunar poles requires an alternative strategy.

Trajectories can be represented as solutions to boundary value problems (BVPs) in terms of ordinary differential equations (ODEs). For a BVP, conditions are specified at the beginning and, also, at the end of the time domain. The specific values of the states at the extremes may or may not be fixed in a BVP with a periodicity constraint; periodicity only requires that the values at the extremes are equal. Differential-corrections schemes [10] and shooting methods (both adaptations of the initial value problem), collocation approaches, and finite-difference methods [1113] are common numerical processes to solve BVPs and are similar in that they all employ linear corrections to a nonlinear system (these three techniques are discussed in greater detail below). Variational methods are also used to solve BVPs when the problem possesses certain minimality properties [14]. Classical perturbation techniques have long supplied analytical approximations for trajectories [5, 1517]. Finally, recent efforts using evolutionary algorithms to solving optimization problems also provide an approach to solving boundary value problems [18, 19].

Shooting methods [10, 12, 2023] and collocation schemes [2, 2431] possess high degrees of accuracy and are frequently employed for trajectory generation and optimization. Shooting methods typically require some knowledge of the design space to initiate the procedure. Collocation schemes do not require a priori knowledge of the solution corresponding to the BVP, but the implementation is often neither intuitive nor straightforward. Furthermore, the partial derivatives of the constraints with respect to the states are nontrivial in collocation schemes and are often approximated numerically [2].

A more rudimentary method for solving BVPs is the finite-difference method (FDM), in which the derivatives that appear in the differential equations are replaced with their respective finite differences and evaluated at node points along the trajectory [12]. The solution process is iterative. The trajectory is discretized, and the equations that represent the relationships at the nodes are solved simultaneously. Using a central difference approximation, the smallest local errors associated with an FDM approach are proportional to Δ𝑡2. However, model accuracy and uncertainty may obviate any precision gains from a high-precision scheme; solutions from those methods may not persevere without additional control in a higher-fidelity model (e.g., a model that includes planetary ephemerides) or in actual flight. If a solution based on a high-precision approach (e.g., shooting or high-degree collocation) is migrated from a model with limited accuracy to one of higher accuracy (e.g., circular restricted three-body model versus ephemeris model or actual flight), it may be no more accurate in that high-fidelity model than a solution from a low-precision approach (e.g., finite-differencing). Moreover, the FDM is simple to formulate and offers an improved understanding of the design space. Partial derivatives are easily obtained via standard analytical techniques, without having to resort to numerical or automatic differentiation (each is often computationally expensive). Once the analyst is familiar with the feasible options, higher-fidelity methods for solving the BVP can be used to refine trajectories or combined with optimization techniques to produce optimal trajectories. Another option is to employ the FDM technique to generate trajectories that fit path constraints to quickly explore a broad design space. Using the MATLAB computing environment, Wawrzyniak and Howell [32] survey over 10 million combinations of initial guesses for the path and control profile, as well as a range of characteristic accelerations and path constraints appropriate to the lunar south pole coverage problem. The time to generate the millions of trajectories for the survey is approximately one week when the algorithms run on eight cores over five platforms.

In the following sections, an FDM is described and applied to solar sail orbits in the Earth-Moon circular restricted three-body (CR3B) system. Also included is a modified FDM for the same application, where the velocity is incorporated as part of the solution at each node along the discretized trajectory. Neither method is strictly a straightforward, “textbook” FDM [1113], since both are augmented with trajectory and control constraints. Additionally, in contrast to a simple two-point BVP, where the states are fixed at the extremes, the trajectory is required to be periodic; nowhere along the arc is the solution prespecified.

A description of the dynamical model in the Earth-Moon CR3B system is followed by an algorithm using the augmented FDMs for generating trajectories. The strategies are based on minimizing the difference between accelerations from the equations of motion (as well as the velocities as another alternative) and the corresponding values numerically derived from positions along a discretized trajectory. An error analysis is also included. A separate study uses these FDMs for surveying the solution space and assessing the required solar sail and spacecraft characteristics necessary for the lunar south pole (LSP) coverage problem [32].

2. System Dynamical Model

The problem that represents this application is defined within the context of the CR3B system, that is, the problem is formulated in a frame, 𝑅, that is rotating with respect to an inertial system, 𝐼. A CR3B model that incorporates the gravity contributions of two primary bodies is geometrically advantageous for understanding the problem. Consistent with McInnes [33], the nondimensional vector equation of motion for a spacecraft at a location 𝐫 relative to the barycenter (center of mass of the primaries) is𝑅𝐚+2𝐼𝝎𝑅×𝑅𝐯+𝑈(𝐫)=𝐚𝑠(𝑡),(1) where the first term is the acceleration relative to the rotating frame (more precisely expressed as 𝑅𝑑2𝐫/𝑑𝑡2, where the left superscript 𝑅 indicates a derivative in the rotating frame) and the second term is the corresponding Coriolis acceleration, which requires the velocity relative to the rotating frame, 𝑅𝐯 (more precisely 𝑅𝑑𝐫/𝑑𝑡). Vectors are denoted with boldface. Derivatives of the position vector, 𝑅𝐯 and 𝑅𝐚, are assumed to be relative to the rotating frame and, consequently, 𝑅 is dropped. The angular-velocity vector, 𝐼𝝎𝑅 (or 𝝎), relates the rate of change the rotating frame with respect to the inertial frame. The applied acceleration, from a solar sail in this case, is indicated on the right side by 𝐚𝑠(𝑡). The pseudogravity gradient, 𝑈(𝐫), combines centripetal and gravitational accelerations:𝑈(𝐫)=(𝝎×(𝝎×𝐫))+(1𝜇)𝑟31𝐫𝟏+𝜇𝑟32𝐫𝟐,(2) where 𝜇 represents the mass fraction of the smaller body, or 𝑚2/(𝑚1+𝑚2), and 𝑟1 and 𝑟2 are the distances from the larger and smaller bodies, respectively, that is, 𝑟1=(𝜇+𝑥)2+𝑦2+𝑧2,𝑟2=(𝜇+𝑥1)2+𝑦2+𝑧2.(3) Solar gravity is neglected in this model. At a distance of 1 AU, an appropriate distance to assume for a sailcraft in the Earth-Moon system, the applied acceleration from a solar sail is modeled as𝐚𝑠̂𝐮(𝑡)=𝛽𝓘(𝑡)2̂𝐮,(4) where ̂𝐮 is the sail-face normal, 𝓘(𝑡) is a unit vector in the sun-to-spacecraft direction, and 𝛽 is the sail characteristic acceleration in nondimensional units. These vectors appear in Figure 1. Observed from the rotating frame, 𝑅, the Sun moves in a clockwise direction about the fixed primaries. The sail mass, 𝑚3, is negligible compared to the masses of the Earth and Moon, which are 𝑚1 and 𝑚2, respectively. The term (̂𝓘(𝑡)𝐮) is also expressed as cos𝛼, where 𝛼 is the sail pitch angle, or the angle between the solar incidence direction and the sail normal.

Figure 1: Earth-Moon system model.

To generate the magnitude of the sail acceleration in dimensional units, 𝑎0, 𝛽 is multiplied by the system characteristic acceleration, 𝑎, which is the relationship between the dimensional and nondimensional acceleration in (1). In fact, 𝑎 is the ratio of the characteristic length, 𝐿 (384,400 km for the Earth-Moon distance), to the square of the characteristic time, 𝑡 (2𝜋𝑡=27.321 days), that is, 𝑎=𝐿𝑡2=2.7307mm/s2.(5) A recent sailcraft design for NASA's Space Technology competition (ST9) that was built by L'Garde possesses overall characteristic acceleration, 𝑎0, of 0.58mm/s2, while the characteristic acceleration of the sail and its support structure alone is closer to 1.70mm/s2 (0.212 to 0.623 in units of nondimensional acceleration, resp.) [34].

The sunlight direction is expressed relative to the rotating frame and is a function of time, that is,̂̂̂𝓘(𝑡)=cos(Ω𝑡)𝐱sin(Ω𝑡)𝐲+0𝐳,(6) where Ω is the ratio of the synodic rate of the Sun as it moves along its path to the system rate, approximately 0.9192. One physical constraint is imposed on the attitude of the spacecraft: the sail normal, ̂𝐮, which is coincident with the direction of the resultant force in an ideal model, is always directed away from the Sun. This is written mathematically aŝ𝓘(𝑡)𝐮cos𝛼max,(7) where 𝛼max is 90.

The sail modeled here is a perfectly reflecting, flat solar sail. Billowing is not incorporated in this force model; however, 𝛼max could be less than 90, as sail luffing (i.e., flapping) is assumed to occur at high pitch angles [35]. Higher fidelity models include optical models [33], parametric models that incorporate billowing in addition to optical effects [33, 36], and realistic models based on finite-element analysis that incorporates optical properties and manufacturing flaws [37]. Optical effects can represent a nonperfectly reflecting solar sail; some energy is absorbed, and some is reflected diffusely as well as specularly. An ideal sail reflects only specularly. In all of these models, the resulting acceleration from a solar sail is not perfectly parallel to the sail-face normal but, instead, is increasingly offset from the sail-face normal as the sail is pitched further from the sunlight direction [33]. Fully accounting for realistic solar sail properties attenuates the sail characteristic acceleration by nearly 25% and places an upper limit on the pitch angle between 50 and 60, depending on the properties of the sail [37]. Nevertheless, this analysis will employ an ideal sail to lend insight into the technology level that is required to solve the LSP coverage problem and into the use of augmented finite-difference methods for generating solar sail trajectories.

3. Augmented Finite-Difference Methods

The derivation for an augmented finite-difference method begins with generic, nonlinear, second-order, two-point BVP that can be represented as̈̇𝐫=𝐟(𝑡,𝐫,𝐫),(8) where 𝐫(𝑡1) and 𝐫(𝑡𝑛) are fixed at the extremes of the time span, [𝑡1,𝑡𝑛] in interval notation [1114]. The solution of (8) is a trajectory, 𝐫(𝑡). To solve the BVP, the classic FDM discretizes the domain into nodes, or epochs, 𝑡1,𝑡2,,𝑡𝑛 and replaces the derivatives in (8) by their respective finite differences. Central-difference approximations (CDAs) are commonly used because they are more accurate than forward or backward differences. The first and second time derivatives of 𝐫(𝑡) are, of course, the velocity and acceleration vectors and are approximated by CDAs as̃𝐯𝑖=𝐫𝑖+1𝐫𝑖1𝑡𝑖+1𝑡𝑖1,̃𝐚(9)𝑖𝐫=2𝑖+1𝐫𝑖𝑡𝑖𝑡𝑖1𝐫𝑖𝐫𝑖1𝑡𝑖+1𝑡𝑖𝑡𝑖𝑡𝑖1𝑡𝑖+1𝑡𝑖1𝑡𝑖+1𝑡𝑖,(10) where 𝐫(𝑡𝑖) is defined as 𝐫𝑖 and 𝑡𝑖(𝑡𝑖1,𝑡𝑖+1) and the symbol “” indicates a numerical approximation. Equation (10) arises from three central-differences: one for the velocity at 𝑡𝑖1/2, which is midway between 𝑡𝑖1 and 𝑡𝑖, another at 𝑡𝑖+1/2, midway between 𝑡𝑖 and 𝑡𝑖+1, and a third using the first two intermediate velocities resulting in the acceleration at 𝑡𝑖. The FDM does not require uniform node placement (mesh refinement and nonuniform node spacing may improve the accuracy of the FDM), but when the time steps between 𝑡𝑖1, 𝑡𝑖, and 𝑡𝑖+1 are fixed at a value of Δ𝑡, then the midpoint of [𝑡𝑖1,𝑡𝑖+1] is 𝑡𝑖 and (9) and (10) simplify tõ𝐯𝑖=𝐫𝑖+1𝐫𝑖1,̃𝐚2Δ𝑡(11)𝑖=𝐫𝑖+12𝐫𝑖+𝐫𝑖1Δ𝑡2.(12) Note that 𝐫𝑛1 is used for 𝐫0 when calculating ̃𝐯1 for a periodic solution. Similarly, 𝐫1 substitutes for 𝐫𝑛 in ̃𝐯𝑛1, so the problem does not include redundant constraints. The calculations for 𝐚1 and 𝐚𝑛1 are similar. The differences between the actual velocity and acceleration and their respective CDAs at 𝑡𝑖 are𝐯𝑖̃𝐯𝑖=Δ𝑡26𝐫𝜌𝑣,𝐚(13)𝑖̃𝐚𝑖=Δ𝑡2𝐫12iv𝜌𝑎,(14) where 𝜌𝑣,𝑎(𝑡𝑖1,𝑡𝑖+1) and 𝐫(𝑡) and 𝐫iv(𝑡) are continuous on [𝑡𝑖1,𝑡𝑖+1]. Derivations of these truncation errors are available in the literature [1113]. The approximations in (11) and (12) replace the first and second time derivatives in (8), and the result is an equation at each epoch of the form𝐫𝑖1+2𝐫𝑖+𝐫𝑖+1Δ𝑡2𝐫=𝐟𝑡,𝐫,𝑖+1𝐫𝑖1.2Δ𝑡(15) Equation (15) may still be nonlinear, but it can be linearized into a system of equations, one for each 𝑡𝑖, excluding the extremes, and solved iteratively. Due to the approximations in (11) and (12), the classic FDM results in a solution that will approximate each component of the true solution to the BVP by an error proportional to Δ𝑡2 at each node [11, 13].

For the current problem, the FDM is augmented to incorporate path and system constraints, as well as a control history. Two variations of this FDM are developed. In the first, the position is the only explicitly discretized trajectory state (FDM-R). The second FDM is formulated to directly solve for position and velocity along the discretized trajectory (FDM-RV).

The FDM-R formulation is based on approximating (1) via the CDAs for velocity and acceleration given by (11) and (12), respectively, that is,̃𝐚𝑖̃𝐯+2𝝎×𝑖𝐫+𝑈𝑖𝐚𝑠𝑡𝑖=0.(16) Note that the velocity term in (16) is replaced with its CDA, ̃𝐯𝑖. To solve the BVP, the trajectory is discretized at 𝑡1,𝑡2,𝑡𝑛. Position is expressed in terms of Cartesian coordinates at each node: 𝐫𝑖=𝑥𝑖𝑦𝑖𝑧𝑖T,(17) where “T” indicates a vector transpose. The control, 𝐮𝑖, that is, the sail pointing vector in this application, as well as a set of 𝑚 slack variables, 𝜂𝑖, are also included at each node. Note that the control, 𝐮𝑖, is not necessarily a unit vector until the FDM converges on a solution. The positions, control, and slack variables at node 𝑖 for the FDM-R are collected into the subvector, 𝐪𝑖=𝐫T𝑖𝐮T𝑖𝜼T𝑖T(6+𝑚),(18) where the length of 𝐪𝑖 is (6+𝑚). The subvectors corresponding to each node are collected into a column vector,𝐪𝐗=T1𝐪T2𝐪T𝑛1𝐪T𝑛T(6+𝑚)𝑛,(19) which represents the discretized trajectory as well as the associated control history and slack variables. The velocities associated with the trajectories generated by the FDM-R can be reconstructed with CDAs of 𝐫(𝑡) during postprocessing.

The FDM-RV process is similarly formulated, but there are two distinct sets of ODEs to solve, that is,̃𝐯𝑖𝐯𝑖̃𝐚=0,(20)𝑖+2𝝎×𝐯𝑖𝐫+𝑈𝑖𝐚𝑠𝑡𝑖=0.(21) Note the difference in the velocity terms in (16) and (21). With this FDM-RV alternative, the subvector is𝐪𝑖=𝐫T𝑖𝐯T𝑖𝐮T𝑖𝜼T𝑖T(9+𝑚),(22) since 𝐯𝑖 is now explicitly part of the solution. This formulation results in an 𝐗 vector with length (9+𝑚)𝑛.

Both the FDM-R and the FDM-RV algorithms constrain the difference between the evaluated EOMs at node 𝑖 and the associated numerically derived approximations to be zero. Additional constraints are required for periodicity, control-direction magnitude, and path characteristics. All of these constraints are dependent on 𝐗 and are collected in the column vector 𝐅(𝐗). An expansion of 𝐅(𝐗) about 𝐗𝑗 yields a linear approximation for 𝐅(𝐗𝑗+1), that is, 𝐅𝐗𝑗+1𝐗=𝐅𝑗𝐗+𝐷𝐅𝑗𝐗𝑗+1𝐗𝑗+𝒪Δ𝐗2,(23) where the 𝐷𝐅(𝐗𝑗) is the Jacobian, 𝜕𝐅(𝐗𝑗)/𝜕𝐗𝑗. Superscripts on vectors refer to iteration number. The initial guess is denoted 𝑗=0; a converged solution by 𝑗=𝑓. Partial derivatives used in the Jacobian are derived analytically, but are not specifically presented here. By assuming that 𝐗𝑗 is in the neighborhood of 𝐅(𝐗)=0 (i.e., 𝐗𝑗+1𝐗𝑗1), (23) is rearranged into a least-norm problem, that is,𝐗𝑗+1=𝐗𝑗𝐗𝐷𝐅𝑗T𝐗𝐷𝐅𝑗𝐗𝐷𝐅𝑗T1𝐅𝐗𝑗.(24) Equations (23) and (24) reflect the Newton-Raphson method of root finding for several variables. For efficient computation, 𝐷𝐅(𝐗) is preallocated and stored in memory as a sparse matrix [38]. Iterations based on (24) continue until some convergence tolerance is met, that is,𝐗𝑗+1𝐗𝑗𝐗𝑗tol.(25) When an initial guess is in the neighborhood of a solution, convergence is quadratic. Specifying tol1×107 is sufficient, since any further iteration approaches double precision. If the initial guess is not in the neighborhood of a solution, 𝐗𝑗+1 bears little resemblance to 𝐗𝑗; the step from 𝐗𝑗 to 𝐗𝑗+1 may even be chaotic. However, the new 𝐗𝑗+1 might be in the neighborhood of an alternate solution and subsequently lead to convergence. The converged discretized solution, 𝐗𝑓, that satisfies these constraints is the trajectory that solves the EOMs to within a theoretical error.

Algebraically, these FDMs are equivalent. In practice, the two formulations yield slightly different results and, sometimes, dramatically different results. This problem is sensitive to the initial guess; because the initial condition may or may not explicitly include velocity, the vector of initial guesses is not the same for the two FDMs. Additionally, the least-norm update in (24) minimizes the norm of the correction to the state vector; the state vectors for the two formulations are different, and the least-norm update for an initial guess that includes velocity states may move the solution to a different region when compared to an update that does not contain velocity states between iterations. The convergence criterion (25) also depends on the norm of the update to the state vector. Because the size of the state vectors differs between the methods, the vector norm of the update will also differ for converged solutions. The current formulations do not control specified states for convergence (e.g., position and not velocity), and other variations on the implementation of these methods will yield different results. Given that the objective is a simple and quick method to approximate the design space, neither strategy should be considered superior to the other; both achieve the stated goal and results from both are insightful.

At each epoch, the converged subvector 𝐪𝑖 differs from the true solution, 𝐪𝑖, due to the truncation error associated with the FDMs and the limited machine precision. The proceeding error analysis is described by Kincaid and Cheney [13, pages 589–592]. For a step size Δ𝑡 that is greater than some value, the error due to truncation in (9)–(14) will dominate the machine error. Assuming that there exists a true set of position states, 𝐫, that solve (1) and (13) as formulated for the FDM-R approach, (8) is written as𝐫𝑖1+2𝐫𝑖+𝐫𝑖+1Δ𝑡2Δ𝑡212𝐫iv𝜌𝑎=𝝎×𝐫𝑖+1𝐫𝑖1+Δ𝑡Δ𝑡26𝝎×𝐫𝜌𝑣𝑈𝐫𝑖+𝐚𝑠𝑡𝑖.(26) The pseudogravity gradient associated with the true solution, 𝐫𝑖, is expanded as a function of an approximate solution, 𝐫𝑖, that is,𝑈𝐫𝑖𝐫=𝑈𝑖+𝐌𝐫𝑖𝐞𝑖𝐞+𝒪T𝑖𝐞𝑖,(27) where 𝐞𝑖=𝐫𝑖𝐫𝑖 and 𝐌(𝐫𝑖) is the Hessian of 𝑈. Using the solution from the FDM-R algorithm, (15) is subtracted from (26), resulting inΔ𝑡2𝐞𝑖1+2𝐞𝑖+𝐞𝑖+1=Δ𝑡1𝐞𝝎×𝑖+1𝐞𝑖1𝐌𝑖𝐞𝑖+Δ𝑡2𝐡𝑖,(28) where 𝐡𝑖=(1/12)𝐫iv(𝜌𝑎)+(1/6)𝐫(𝜌𝑣). It is assumed that 𝐫(𝑡) is continuously differentiable to fourth order. Rearranging and multiplying by Δ𝑡2 yields(𝐈+Δ𝑡𝝎×)𝐞𝑖1+2𝐈Δ𝑡2𝐌𝑖𝐞𝑖+(𝐈Δ𝑡𝝎×)𝐞𝑖+1=Δ𝑡4𝐡𝑖,(29) where 𝐈 is the identity matrix and 𝝎× is the skew-symmetric gyroscopic matrix. Some terms in (29) can be simplified to a form that will be useful later, that is(𝐈+Δ𝑡𝝎×)+2𝐈Δ𝑡2𝐌𝑖+(𝐈Δ𝑡𝝎×)=Δ𝑡2𝐌𝑖.(30) Now, let 𝝀 correspond to the difference element 𝐞𝑖 that possesses the largest magnitude. Equation (29) reduces to an upper bound on 𝝀,(𝐈+Δ𝑡𝝎×)𝝀+2𝐈Δ𝑡2𝐌𝑖𝝀+(𝐈Δ𝑡𝝎×)𝝀Δ𝑡4𝐡𝑖.(31) Using (30), (31) simplifies toΔ𝑡2𝐌𝑖𝝀Δ𝑡4𝐡𝑖,(32)𝐌𝑖𝝀Δ𝑡2𝐡𝑖,(33)𝝀Δ𝑡2𝐌𝑖1𝐡𝑖.(34) The inequality in (34) represents three scalar equations. Thus, the position error corresponding to the FDM-R formulation is 𝒪(Δ𝑡2) as Δ𝑡0. To calculate the velocity from the FDM-R algorithm, the CDA from (11) is used, also with an error 𝒪(Δ𝑡2). The error analysis for the alternate FDM-RV algorithm incorporates the fact that the velocity, 𝐯=𝐟(𝑡,𝐫), possesses errors that can also be expressed by (13), and, therefore, the derivation continues with the same steps as those for the FDM-R formulation beginning with (26).

The augmented finite-difference methods sacrifice precision for simplicity. At each node, the error in the trajectory is proportional to Δ𝑡2, significantly less precise than common collocation methods [26, 39]. In this analysis, Δ𝑡=0.067, or 𝑛=101, which translates to an error of 0.452% of the Earth-Moon distance, approximately 1740 km, in each direction. To improve the precision of the FDM, the step size between nodes, Δ𝑡, can be decreased. The Jacobian 𝐷𝐅(𝐗) becomes quite large, imposing a significant computational cost. Furthermore, the smaller the step size, the larger the machine error [12]. Alternatively, Richardson extrapolation, or extrapolation to the limit, could be employed to further minimize the error of the solutions [12]. Incorporating a process to determine the specific step size that minimizes some combination of truncation and machine errors or incorporating an additional process, like Richardson extrapolation, adds complexity to a method that is formulated to be simple to understand and implement, as well as quick to yield results.

Prior to any analysis, the design space for this problem is not well known, so a trajectory that is precise to within Δ𝑡2 relative to an actual solution is meaningful. If the goal is a general understanding of the design space, accomplished in a relatively quick analysis, then these results can be very insightful.

4. Algebraic Constraint Vector: F(X)

At the core of the augmented finite-difference methods is the algebraic constraint vector, 𝐅(𝐗), which contains the algebraic approximations to the equations of motion and other constraints necessary for the desired periodic solar sail trajectory. In this formulation, each element in 𝐅(𝐗) should be zero to simultaneously satisfy the ODEs and path constraints. For a periodic solar sail trajectory, subject to path constraints 𝐠(𝐗),𝐅(𝐗)=Δ𝐚(𝐗)TΔ𝐯(𝐗)T𝐓(𝐗)T𝑦1𝐍(𝐗)T𝐠(𝐗)TT.(35) The presence of Δ𝐯(𝐗) depends on whether the formulation is the FDM-R or FDM-RV, and the length of the 𝐅(𝐗) vector is either (2+(4+𝑚)𝑛) or (2+(7+𝑚)𝑛), accordingly. The order of these elements within the vector is arbitrary and no difference in performance is apparent when 𝐅(𝐗) and 𝐗 are rearranged such that 𝐷𝐅(𝐗) is a sparse, banded-diagonal matrix. Using the configuration in (35), each element set appears as block diagonal in the corresponding Jacobian, 𝐷𝐅(𝐗). Each set is subsequently discussed, and the corresponding size of each set is indicated in a subscript outside of the braces.

The first element set in (35) is the difference in acceleration. The acceleration at a given epoch, 𝑖, as evaluated from the equations of motion, depends only on the position, velocity, and control, 𝐚𝑖=𝐟(𝐫𝑖,𝐯𝑖,𝐮𝑖). The numerically derived acceleration, ̃𝐚𝑖 from (10) or (12), depends on the state at multiple epochs. A valid trajectory possesses the same acceleration whether computed from the EOMs or from the CDAs, within the truncation error. Therefore, the difference between accelerations resulting from the EOMs and those from numerical calculation should nominally be zero for a converged solution; thus, the difference forms the first set of equality constraints in the composite constraint vector:𝐚Δ𝐚(𝐗)=1𝐫1,𝐯1,𝐮1̃𝐚1𝐚𝑛1||(𝐫,𝐯,𝐮)𝑛1̃𝐚𝑛13(𝑛1).(36) Likewise, the difference between the velocity from the EOMs and the numerically derived velocity is zero for a valid trajectory. Velocity from the EOMs at each epoch, 𝐯𝑖, is available from 𝐗 (see (19) and (22)). This difference forms a second set of equality constraints (only in the FDM-RV algorithm) in 𝐅(𝐗):𝐯Δ𝐯(𝐗)=1̃𝐯1𝐯𝑛1̃𝐯𝑛13(𝑛1),(37) where 𝐯𝑖 may be approximated as ̃𝐯𝑖. These differences in acceleration and velocity may also be denoted as defects and are located at the node corresponding to each 𝑡𝑖. In an initial guess for the trajectory state vector, the defects Δ𝐚𝑖 and Δ𝐯𝑖 at each node will most likely not be zero. The Newton-Raphson iteration process adjusts the path of the trajectory to resolve these differences.

The next constraints enforce periodicity and unit length of the control vector. To enforce periodicity, the goal is an originating and final state vector that are equal, which is represented as a set of constraints in the algebraic constraint vector, 𝐅(𝐗), that is,𝐓𝐫(𝐗)=𝑛𝐫1𝐯𝑛𝐯1𝐮𝑛𝐮1𝜼𝑛𝜼1(69+𝑚),(38) where 𝑚 is the number of path constraints at the first and last node, expressed as slack variables. The length of 𝐓(𝐗) depends on the application of an FDM-R or an FDM-RV approach. Although the only boundary condition is periodicity, for consistency the position at the boundary is constrained to be in the 𝑥𝑧-plane. Therefore, the 𝑦-coordinate at the first point along the trajectory is constrained, 𝑦1=0. In practice, without this constraint, 𝑦1 remained very close to zero, approximately 40 km from the 𝑥𝑧 plane once the solution is converged. With this constraint, 𝑦1 is reduced to 0.004 mm from the plane. Next, the magnitude of the control vector must remain of unit length in the iteration process𝐮𝐍(𝐗)=T1𝐮1𝐮1T𝑛1𝐮𝑛11(𝑛1).(39) Formulating the constraint in terms of 𝐮T𝑖𝐮𝑖1 as opposed to 𝐮T𝑖𝐮𝑖1 results in a simpler partial derivative while retaining the intended effect; this constraint is incorporated into the implementation of the FDM since the control history in this analysis is described in terms of a unit vector.

Path constraints are included as inequality constraints, which are converted to equality constraints by use of slack variables, a successful numerical adaptation from nonlinear programming [2, 40]. The slack variables are incorporated into 𝐗 and are associated with the other state elements at their particular epoch, 𝑖. In this analysis, the elevation angle constraint, 𝐸min, maintains the visibility of the spacecraft from an outpost located near the south pole of the Moon, and the spacecraft altitude constraint, 𝐴max, is imposed for radio power restrictions. Altitude is defined as the distance from the lunar south pole, that is,𝐴𝑖=𝑥𝑖1+𝜇2+𝑦2𝑖+𝑧𝑖+𝑅𝑚2.(40) The third path constraint requires that the sail-face normal, or control 𝐮𝑖, is always directed away from the Sun (the sunlight vector is 𝓘𝑖), or 𝛼max=90 in (7). Of the inequality constraints, only this attitude requirement is mandated. Together, these three inequality constraints are written as 𝐸min𝐸𝑖𝑧arcsin𝑖+𝑅𝑚𝐴𝑖𝐴,(41a)max𝐴𝑖,(41b)cos𝛼max𝓘𝑖𝐮𝑖.(41c)For the given problem and model, adding a path constraint to avoid the penumbra and umbra of the Earth or the Moon's shadows is unnecessary because of the elevation constraint. A shadow constraint could be added for another application or shadowing effects could be incorporated into the dynamical model directly [41]. Additional inequality path constraints could include limits on the body turn rates and accelerations governed by the attitude control system of the spacecraft. The associated slack variables are squared and added to the inequality constraints, resulting in the following equality constraints: 𝐠𝑖𝐫𝑖,𝜼𝑖=sin𝐸min+𝑧𝑖+𝑅𝑚𝐴𝑖+𝜂2𝐸,𝑖𝐴𝑖𝐴max+𝜂2𝐴,𝑖cos𝛼max𝓘𝑖𝐮𝑖+𝜂2,𝑖𝑚.(42) Combining the path constraints at all epochs results in 𝑚(𝑛1) total elements in 𝐠(𝐗). When using the path constraints in (42), 𝑚=3. The first two path constraints in (42) as well as the periodicity constraint in (38) are illustrated in Figure 2. The cos𝛼max constraint appears in Figure 1. The sailcraft in Figure 2 orbits below the Moon. Feasible solutions could exist that do not orbit below the Moon, but do orbit below either the 𝐿1 or the 𝐿2 point.

Figure 2: Path constraints for an orbit below the Moon (Moon image from

Only the first (𝑛1) epochs are required in formulating Δ𝐚, Δ𝐯, 𝐍, and 𝐠(𝐗) in the constraint vector 𝐅(𝐗). Due to periodicity, the first and 𝑛th epochs are identical. Activating the 𝑛th epoch yields a problem that is overconstrained and the Jacobian, 𝐷𝐅(𝐗), is not full rank.

As mentioned, the Jacobian is a sparse matrix. A diagram of the sparsity pattern based on partial derivatives of 𝐅(𝐗) with respect to 𝐗 for the FDM-RV formulation appears in Figure 3. Diagonals corresponding to the specific element sets in 𝐅(𝐗) are labeled in the figure. For a scenario employing 101 nodes, the size of the Jacobian is 1013 rows by 1212 columns where approximately 0.38% of the elements are nonzero. The sparsity pattern for the FDM-R formulation is similar to that of the FDM-RV, except that the diagonal corresponding to 𝐯(𝐗) does not exist in the FDM-R option. The size of the Jacobian for the FDM-R option associated with 101 nodes is 710 rows by 909 columns; approximately 0.62% of the entries are nonzero.

Figure 3: Sparsity pattern for the Jacobian 𝐷𝐅(𝐗) in the FDM-RV formulation.

The advantage of a finite-difference approach is its ease of implementation and the speed improvements enabled by its simplicity. Partial derivatives are easily accessible via analytical derivation, especially for an idealized force model such as (1). As complexity of the formulation increases, which is the case with collocation, for example, analysts rely on numerical or automatic differentiation to generate partial derivatives. MATLAB supplies the function NUMJAC.M, which numerically approximates a derivative using a forward difference approximation [42]. Third-party software for MATLAB exists for automatic (a.k.a. algorithmic) differentiation. TOMLAB Optimization's MAD (MATLAB Automatic Differentiation) suite employs automatic differentiation [43], and Shampine's PMAD (Poor Man's Automatic Differentiation) exploits complex step differentiation [44]. Both MAD and PMAD result in highly accurate derivatives, on the order of machine precision when compared to the partial derivatives determined analytically, while derivatives calculated by NUMJAC.M are accurate to 108. Because analytical derivatives are directly encoded into the algorithm, computation times are significantly faster when compared to the other options listed. Table 1 lists the computations times for generating 𝜕𝑈/𝜕𝐫, the most computationally intensive step in constructing an entry in the Jacobian, for the four differentiation strategies. Clearly, if analytical derivatives are easily available, they should be employed. Furthermore, if the subroutine to generate the derivatives analytically is programmed in a lower-level language (e.g., C or FORTRAN) and then complied as a MATLAB Executable (MEX), further time savings are realized (with compiled code the computation of 𝜕𝑈/𝜕𝐫 takes less than 0.001 seconds). The formulation of the FDMs for this analysis exploits the computational advantage of analytical derivatives and compiled MEX code.

Table 1: Computation times for various differentiation strategies.

In summary, a finite-difference method yields a solution for the path by replacing the path derivatives with their finite-difference approximations. The differences between the approximate derivatives and the derivatives determined by evaluating the equations of motion at a specific point along a path are minimized by iteratively solving a linearized system of equations. Inequality path constraints are added to this set of equality constraints by way of slack variables, and subsequently augmenting the constraint vector and linearized system of equations, resulting in a feasible path.

5. Results

The objective when employing a finite-difference method is to quickly generate a feasible trajectory where system constraints and sail characteristics are given. Because model simplification into the CR3B system ignores lunar librations, the inclined orbits of the primaries, additional perturbations to the orbit, and lunar surface features, a conservative minimum elevation angle of 𝐸min=15 is selected. A generous maximum-altitude constraint of one Earth-Moon distance (384,400 km) is also selected (in practice most results are an order of magnitude closer to the lunar south pole and this constraint is not active) to establish the constraints in (42).

As mentioned, a recent sail design possesses a characteristic acceleration of 0.58mm/s2 [34]. JAXA's IKAROS is the first mission to have flown using a sail as its only propulsive device, and NASA recently deployed the NanoSail-D2. Both of these sailcraft are relatively small (200 m2 and 10 m2, resp.) and are designed to demonstrate in-space deployment of the sail. Designs for larger sailcraft exist, and future solar sail spacecraft are likely to be hybridized with other propulsion devices [45]. However, it is instructive to demonstrate the FDMs with characteristic accelerations for sail spacecraft that may be available with near-future technological improvements. An assessment of the technology level (i.e., characteristic acceleration) is required for a sail-only LSP mission. The focus of this analysis is not to demonstrate that a sail-only LSP mission is feasible based on current technology. Instead, the results of this investigation demonstrate the use of finite-difference methods for the sample application of an LSP mission. Therefore, a high characteristic acceleration, that is, 1.70mm/s2, is used to generate example orbits.

To demonstrate the FDM-R and FDM-RV approaches, three sets of initial guesses for the path and control history are selected to initialize the process. The initial guesses corresponding to the slack variables are always established such that 𝐠0(𝐗)=𝟎 (the superscript “0” indicates an initial guess). For each orbit, the initial guess of the control history is one that maximizes the out-of-plane force contributed by the sail along the initial trajectory. Derived analytically by McInnes [33, pages 115–118, 223-224], the sail-pitch angle that maximizes the out-of-plane thrust below the Moon is 𝛼=35.26. The initial guess for the thrust vector is then̂𝐮0𝑖=cosΩ𝑡𝑖cos𝛼̂𝐱+sinΩ𝑡𝑖cos𝛼̂𝐲sin𝛼̂𝐳.(43) The initial guess for each trajectory is a circle offset below the Moon. Three sample initial guesses appear in Figure 4. The dark-blue circle with a radius of 59,000 km is offset below the Moon by 23,000 km. The radius of the light-blue circle is 14,000 km and the 𝑧-offset is 54,000 km. The dimensions of the red circle are 67,000 km by 75,000 km. As previously described, the system is formulated in the Earth-Moon circular restricted three-body system, where the reference frame is fixed to the Earth and the Moon. Therefore, the Sun is initially along the negative 𝑥-axis and moves in this reference frame clockwise around the Earth and the Moon, per (6). In this frame, the spacecraft is initially on the opposite side of the Moon from the Sun and moves clockwise around the Moon. The arrows on the orbits represent the initial guess for the direction of the sail-face normal, ̂𝐮, at the specified point along the orbit. Note that their direction is generally away from the Sun. The arrows are spaced in time by a little more than one day.

Figure 4: Three example initial guesses.

When the estimates from Figure 4 initialize the FDM-R algorithm, the three orbits in Figure 5 result. While the dark-blue orbit is not far from its initial approximation, the red orbit and the light-blue orbit are obviously dissimilar from their respective initial guesses. This indicates that a good initial guess for the finite-difference method is not essential for converging on a trajectory that solves the equations of motion. When the same initial guesses are used as inputs for the FDM-RV process, the trajectories in Figure 6 result. For both the FDM-R and FDM-RV algorithms, the dark-blue orbits converged in less than 10 iterations. The other two orbits required between 15 and 20 iterations for both methods. When individual iterations are examined for these two orbits, it is observed that the general shape of the orbit emerges quickly from the initial guess and the large number of iterations results from small attitude and path changes at the extremities in the ±𝑦 directions to accommodate the elevation constraint. All solutions converge quadratically over the final four or five iterations.

Figure 5: Three converged orbits from the FDM-R process.
Figure 6: Three converged orbits from the FDM-RV process.

For all trials of the FDM-R and FDM-RV formulations (and in the collocation formulation described in the next section) it is observed that the corrections to the components of the initial guess of the control history are not smooth at the first and 𝑛th point; this phenomenon is apparent in the control history from the seventh-degree Gauss-Lobatto collocation formulation presented in Figure 7 of Ozimek et al. [2]. All corrections to the discretized path and attitude profiles are done simultaneously, and the first point and 𝑛th points are the only two points in this simulation that are constrained to be equal to each other to ensure periodicity. Because the remaining points in the attitude profile do appear to be continuous and smooth, a postfit modification is applied such that the first and 𝑛th points are shifted by interpolating with neighboring points in the periodic profile. This postfit modification is only applied when the Newton-Raphson solver is not within the radius of convergence (i.e., when the convergence criterion in (25) exceeds 0.1). Control profiles from subsequent iterations appear smooth thereafter.

Figure 7: Three converged orbits from the Hermite-Simpson collocation process.

Although the dark-blue and red orbits appear to be similar in both the FDM-R and FDM-RV results, they differ by up to 8000 km and 1800 km, respectively, at certain epochs along the respective paths. The difference in the light-blue orbits in Figures 5 and 6 strongly suggests that the different methods, FDM-R and FDM-RV, can produce different solutions to the equations of motion using the same initial guess. While the theoretical errors in the solutions from the respective algorithms are equivalent, the differences in the formulation of the FDM (either -R or -RV) drive the evolution of the trajectory from iteration to iteration. It is not appropriate to suggest that the solution from one method is more “correct” than the other; both methods theoretically solve the equations of motion to within their errors bounds, and both methods produce orbits that are not obvious when the design space is not well understood.

6. Comparison to a Collocation Method

As stated, prior to any analysis, the design space for this problem is unknown. Entry into the design space via analytic techniques is difficult because of coupling of the sail attitude and thrust vector. However, based on the FDM algorithms, solutions precise to 𝒪(Δ𝑡2), or 1740 km for this step size, in each direction emerge. In work by Ozimek et al. [2], a Hermite-Simpson collocation method for the LSP coverage problem that has a theoretical local error on the order of 𝒪(Δ𝑡5), or 0.5 km for the same step size as that from the FDM algorithms. The Hermite-Simpson collocation method is adapted here as a basis to evaluate the accuracy of the FDM algorithms.

Using the same initial guesses as those represented in Figure 4 to seed the collocation method, the three orbits in Figure 7 result. The three orbits from the Hermite-Simpson approach are clearly similar to those produced by the FDM-RV process and only differ by, at most, 3000 km for the red orbit, 450 km for the light-blue orbit, and 70 km for the dark-blue orbit. The dark-blue and red paths from the Hermite-Simpson process are also similar to the dark-blue and red paths from the FDM-R approach (differing by, at most, 8000 km and 3500 km, resp.), but the light-blue paths are significantly different between the two strategies. Since the theoretical local error in the Hermite-Simpson method is on the order of 𝒪(Δ𝑡5), or 0.5 km for this step size, these differences between the similar orbits are not surprising [30]. Given that the initial guesses for the three sample orbits do not resemble their associated converged trajectories, it is remarkable that the converged trajectories from all three numerical processes are similar.

There is one noteworthy difference between the Hermite-Simpson collocation method employed by Ozimek et al. and the adaptation used for this analysis. Ozimek et al. employ a small step along the imaginary axis to calculate the partial derivatives related to the defects. For the current analysis, a simple central difference approximation (CDA) represents the same partial derivatives. Both are numerical schemes for computing difficult partial derivatives and subject to roundoff errors, but the CDA is also subject to truncation error, while the imaginary-step method is not. Thus, the error for the CDA employed in the Hermite-Simpson collocation method is minimized at a step size greater than zero (on the order of 108), but the error for the imaginary-step method approaches the roundoff error as the step size approaches zero (on the order of machine precision) [46]. This difference in precision should not have large ramifications for the process to converge, but must be noted. Furthermore, because the partial derivatives are numerically formulated for the Hermite-Simpson collocation scheme in this analysis, the three orbits in Figure 7 took approximately 2 seconds to generate. The orbits in Figures 5 and 6 took approximately 1 second to generate.

The results from any procedure can be used to seed another, usually more precise, process. The solutions produced by the FDM-R and FDM-RV algorithms are used to initialize the Hermite-Simpson scheme. All orbits reconverge quadratically in three or four steps when the Hermite-Simpson strategy is incorporated. The differences in position and velocity, as well as the direction of the sail-face normal, along the 29.5 day trajectories appear in Figure 8. Each color in the figure corresponds to the respective “initial” orbits appearing in Figures 5 and 6. The differences are well within the theoretical errors for the finite-difference methods, and, as such, the augmented finite-difference methods produce accurate and appropriate reference orbits for preliminary mission analysis.

Figure 8: Differences between the solutions from the Hermite-Simpson method and their initial guesses supplied by the FDM-R (a) and FDM-RV (b) solutions.

7. Conclusions

Augmented finite-difference methods are useful tools for generating trajectories in poorly understood dynamical regimes, such as the mechanics of flying a solar sail in a multibody system. These schemes are simple to understand and implement, notably in the presence of path constraints. While the theoretical errors of finite-difference methods may be larger than other similar methods that examine the trajectory as a whole, such as collocation, finite-difference methods can provide a reasonable entrance to the design space of a complicated nonlinear problem. They also quickly return results in a regime where little intuition concerning the trajectory and sail-angle history exist. Because of its speed and simplicity, this method may also serve as the basis for generating a large survey of orbit options. Once a viable trajectory is uncovered by the finite-difference method, other techniques may be employed to refine it further. The FDM approach developed for this analysis need not be applied to just solar sail problems, but should be applicable to a large variety of mission design problems.


The authors would like to thank Purdue Professor William Crossley for his explanation of optimization techniques in the literature used to solve similar problems and Professor Ahmed Sameh, also from Purdue University, for his insights on numerical methods used to solve boundary-value problems. Former Purdue graduate students Daniel Grebow and Zubin Olikara were invaluable to the initial development and understanding of the numerical methods examined in this paper. The authors also are grateful for the comments from the Editor and reviewers.


  1. A. Ely and E. Lieb, “Constellations of elliptical inclined lunar orbits providing polar and global coverage,” in Advances in the Astronautical Sciences, vol. 123, pp. 1447–1462, Lake Tahoe, Calif, USA, 2005, paper no. AAS 05-343. View at Google Scholar
  2. M. T. Ozimek, D. J. Grebow, and K. C. Howell, “Design of solar sail trajectories with applications to lunar south pole coverage,” Journal of Guidance, Control, and Dynamics, vol. 32, no. 6, pp. 1884–1897, 2009. View at Publisher · View at Google Scholar · View at Scopus
  3. G. G. Wawrzyniak and K. C. Howell, “The solar sail lunar relay system: an application of solar sails in the Earth–Moon system,” in Proceedings of the 59th International Astronautical Congress, Glasgow, Scotland, September 2008, paper no. IAC-08.C1.3.14.
  4. D. J. Grebow, M. T. Ozimek, and K. C. Howell, “Design of optimal low-thrust lunar pole-sitter missions,” in Advances in the Astronautical Sciences, vol. 134, pp. 741–760, Savannah, Ga, USA, 2009, paper no. AAS 09-148. View at Google Scholar
  5. R. W. Farquhar, “The control and use of libration point satellites,” Tech. Rep. TR-R-346, National Aeronautics and Space Administration, February 1970. View at Google Scholar
  6. C. R. McInnes, “Solar sail trajectories at the lunar L2 Lagrange point,” Journal of Spacecraft and Rockets, vol. 30, no. 6, pp. 782–784, 1993. View at Google Scholar · View at Scopus
  7. C. R. McInnes, A. J. McDonald, J. F. L. Simmons, and E. W. MacDonald, “Solar sail parking in restricted three-body systems,” Journal of Guidance, Control, and Dynamics, vol. 17, no. 2, pp. 399–406, 1994. View at Google Scholar · View at Scopus
  8. J. Simo and C. R. McInnes, “Solar sail trajectories at the Earth–Moon Lagrange points,” in Proceedings of the 59th International Astronautical Congress, Glasgow, Scotland, 2008, paper no. IAC-08.C1.3.13.
  9. J. Simo and C. R. McInnes, “Solar sail orbits at the Earth-Moon libration points,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 12, pp. 4191–4196, 2009. View at Publisher · View at Google Scholar · View at Scopus
  10. K. C. Howell and H. J. Pernicka, “Numerical determination of Lissajous trajectories in the restricted three-body problem,” Celestial Mechanics, vol. 41, no. 1–4, pp. 107–124, 1987. View at Publisher · View at Google Scholar · View at Scopus
  11. H. B. Keller, Numerical Methods for Two-Point Boundary-Value Problems, Blaisdell, Waltham, Mass, USA, 1968.
  12. S. D. Conte and C. de Boor, Elementary Numerical Analysis, McGraw-Hill, New York, NY, USA, 1980.
  13. D. R. Kincaid and E. W. Cheney, Numerical Analysis: Mathematics of Scientific Computing, American Mathematical Society, Providence, RI, USA, 2002.
  14. J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, Springer, New York, NY, USA, 3rd edition, 2002.
  15. E. W. Brown, An Introductory Treatise on Lunar Theory, Cambridge University Press, London, UK, 1896.
  16. D. L. Richardson, “Analytic construction of periodic orbits about the collinear points,” Celestial Mechanics, vol. 22, no. 3, pp. 241–253, 1980. View at Publisher · View at Google Scholar
  17. S. H. Strogatz, Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering, Westview Press, Cambridge, Mass, USA, 2001.
  18. D. Izzo, “Global optimization and space pruning for spacecraft trajectory design,” in Spacecraft Trajectory Optimization, B. A. Conway, Ed., Cambridge University Press, New York, NY, USA, 2010. View at Google Scholar
  19. M. Pontani and B. A. Conway, “Swarming theory applied to space trajectory optimization,” in Spacecraft Trajectory Optimization, B. A. Conway, Ed., Cambridge University Press, New York, NY, USA, 2010. View at Google Scholar
  20. R. S. Wilson and K. C. Howell, “Trajectory design in the Sun-Earth-Moon system using lunar gravity assists,” Journal of Spacecraft and Rockets, vol. 35, no. 2, pp. 191–198, 1998. View at Google Scholar · View at Scopus
  21. B. G. Marchand, K. C. Howell, and R. S. Wilson, “Improved corrections process for constrained trajectory design in the n-body problem,” Journal of Spacecraft and Rockets, vol. 44, no. 4, pp. 884–897, 2007. View at Publisher · View at Google Scholar · View at Scopus
  22. O. von Stryk and R. Bulirsch, “Direct and indirect methods for trajectory optimization,” Annals of Operations Research, vol. 37, no. 1, pp. 357–373, 1992. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  23. J. T. Betts, “Survey of numerical methods for trajectory optimization,” Journal of Guidance, Control, and Dynamics, vol. 22, no. 2, pp. 193–207, 1998. View at Google Scholar · View at Scopus
  24. R. D. Russell and L. F. Shampine, “A collocation method for boundary value problems,” Numerische Mathematik, vol. 19, no. 1, pp. 1–28, 1972. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  25. E. J. Doedel, “Finite difference collocation methods for nonlinear two point boundary value problems,” SIAM Journal on Numerical Analysis, vol. 16, pp. 173–185, 1979. View at Google Scholar
  26. C. R. Hargraves and S. W. Paris, “Direct trajectory optimization using nonlinear programming and collocation,” Journal of Guidance, Control, and Dynamics, vol. 10, no. 4, pp. 338–342, 1987. View at Google Scholar · View at Scopus
  27. P. J. Enright and B. A. Conway, “Optimal finite-thrust spacecraft trajectories using collocation and nonlinear programming,” Journal of Guidance, Control, and Dynamics, vol. 14, no. 5, pp. 981–985, 1991. View at Google Scholar · View at Scopus
  28. P. J. Enright and B. A. Conway, “Discrete approximations to optimal trajectories using direct transcription and nonlinear programming,” Journal of Guidance, Control, and Dynamics, vol. 15, no. 4, pp. 994–1002, 1992. View at Google Scholar · View at Scopus
  29. S. Tang and B. A. Conway, “Optimization of low-thrust interplanetary trajectories using collocation and nonlinear programming,” Journal of Guidance, Control, and Dynamics, vol. 18, pp. 599–604, 1995. View at Google Scholar
  30. A. L. Herman and B. A. Conway, “Direct optimization using collocation based on high-order Gauss-Lobatto quadrature rules,” Journal of Guidance, Control, and Dynamics, vol. 19, no. 3, pp. 592–599, 1996. View at Google Scholar · View at Scopus
  31. R. G. Melton, “Comparison of direct optimization methods applied to solar sail problems,” in Proceedings of the AIAA/AAS Astrodynamics Specialists Conference and Exhibit, Monterey, Calif, USA, August 2002, paper no. AIAA 2002-4728.
  32. G. G. Wawrzyniak and K. C. Howell, “Investigating the design space for solar sail trajectories in the Earth—Moon system using augmented finite-difference methods,” Purdue University Document, School of Aeronautics and Astronautics, February 2011. View at Google Scholar
  33. C. R. McInnes, Solar Sailing: Technology, Dynamics and Mission Applications. Space Science and Technology, Springer, New York, NY, USA, 1999.
  34. D. Lichodziejewski, B. Derbes, D. Sleight, and T. Mann, “Vacuum deployment and testing of a 20m solar sail system,” in Proceedings of the 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Newport, RI, USA, May 2006, paper no. AIAA-2006-1705.
  35. B. Derbes, D. Lichodziejewski, and G. Veal, “A “yank and yaw” control system for solar sails,” in Advances in the Astronautical Sciences, vol. 19, pp. 2893–2907, 2004. View at Google Scholar
  36. L. Rios-Reyes and D. J. Scheeres, “Solar-sail navigation: estimation of force, moments, and optical parameters,” Journal of Guidance, Control, and Dynamics, vol. 30, no. 3, pp. 660–668, 2007. View at Publisher · View at Google Scholar · View at Scopus
  37. B. Campbell, An analysis of thrust of a realistic solar sail with focus on a flight validation missionin a geocentric orbit, Ph.D. Thesis, The George Washington University, Washington, DC, USA, 2010.
  38. The MathWorksTM, Inc, “MATLAB R, Version 7.11.0 (R2010b),” August 2010,
  39. A. L. Herman, Improved collocation methods with application to direct trajectory optimization, Ph.D. Thesis, University of Illinois at Urbana-Champaign, 1995.
  40. J. T. Betts, Practical Methods for Optimal Control Using Nonlinear Programming, SIAM, Philadelphia, Pa, USA, 2001.
  41. C. R. Ortiz Longo and S. L. Rickman, “Method for the calculation of spacecraft umbra and penumbrashadow terminator points,” Tech. Rep. 3547, Lyndon B. Johnson Space Center, Houston,Tex, USA, 1995. View at Google Scholar
  42. D. E. Salane, “Adaptive routines for forming Jacobians numerically,” Tech. Rep. SAND86-1319, Sandia National Laboratories, Albuquerque, NM, USA, 1986. View at Google Scholar
  43. S. A. Forth and M. Edvall, User Guide for MAD—A Matlab Automatic Differentiation Package, TOMLAB/MAD, Version 1.4 The Forward Mode, Engineering Systems Department, Defence College of Management & Technology, Cranfield University, Shrivenham, UK, 2007.
  44. L. F. Shampine, “Accurate numerical derivatives in MATLAB,” ACM Transactions on Mathematical Software, vol. 33, no. 4, Article ID 1268781, 2007. View at Publisher · View at Google Scholar · View at Scopus
  45. M. Macdonald and C. McInnes, “Solar sail mission applications and future advancement,” in Proceedings of the 2nd International Symposium on Solar Sailing, pp. 1–26, New York City College of Technology, City University of New York, New York, NY, USA, July 2010.
  46. J. Martins, I. Kroo, and J. Alonso, “An automated method for sensitivity analysis using complex variables,” in Proceedings of the 38th Aerospace Sciences Meeting and Exhibit, Reno, Nev, USA, January 2000, paper no. AIAA 2000-0689.