# Some New Locally Optimal Control Laws for Sailcraft Dynamics in Heliocentric Orbits

Authors: F. A. Abd El-Salam
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The concept of solar sailing and its developing spacecraft is presented. The gravitational and solar radiation forces are considered. The effect of source of radiation pressure and the force due to coronal mass ejections and solar wind on the sailcraft configurations is modeled. Some analytical control laws with some mentioned input constraints for optimizing sailcraft dynamics in heliocentric orbit using lagrange’s planetary equations are obtained. Optimum force vector in a required direction is maximized by deriving optimal sail cone angle. Ignoring the absorbed and diffusely reflected parts of the radiation, some special cases are obtained. New control laws that maximize thrust to obtain certain required maximization in some particular orbital element are obtained.


---

## Body

## 1. Introduction

The theoretical concept of solar sails is old-standing and dates back to Johannes Kepler when he observed that comet tails point away from the Sun and he suggested that the Sun caused this propulsion effect. After that it transferred into practice by the Soviet pioneers of astronautics.  Solar sails are a form of spacecraft propulsion using a combination of light (radiation pressure) and high speed ejected gasses from a star (e.g., solar wind and coronal mass ejection) to push large ultra-thin mirrors to high speeds.  This concept of using solar radiation pressure as a means of propulsion for space vehicles was first introduced in the 1920s by the father of Russian astronautics, Tsiolkovsky [1] and Tsander [2]. It is also used in the spacecraft attitude dynamics. Rizvi [3] he developed a control method for the solar sail normal vector to trace a desired circular coning trajectory at orbit rate. He finally concludes that the control torques can be applied to the sailcraft to enable orbit rate cone tracing of the sail normal and yield the desired orbital effects.In his book, Wright [4] presents a detailed analysis on some possible solar sail applications. During his time at the Jet Propulsion Laboratory (JPL), Wright was actively involved in the planning of a rendezvous mission to comet Halley using solar sail technology. In 1977, a solar electric propulsion concept was selected instead, primarily because of technology maturity. Not long thereafter, the Halley rendezvous mission was dropped by NASA. All these concepts were first experienced with the Mariner 10 mission to Mercury and Venus. Mariner 10 was also the first spacecraft to use a gravity assist trajectory, accelerating as it entered the gravitational influence of Venus, then using the planet's gravity field to move onto a slightly different course to reach Mercury. Since then, there have been several attempts to realize a solar sail mission.The sail concept is performed by gaining momentum from an ambient source, solar electromagnetic radiation.  Using momentum gained only by reflecting ambient sunlight, the sail slowly but continuously accelerated to accomplish a wide-range of potential missions. Light sails could also be driven by energy beams to extend their range of operations, which is strictly beam sailing rather than solar sailing. Solar sailcrafts offer the possibility of low-cost operations combined with long operating lifetimes. Since they have few moving parts and use no propellant, they can potentially be used numerous times for delivery of payloads.The mechanism of this concept refering to the momentum carried by individual photons is extremely small. Thus, to provide a suitably large momentum transfer, we require the sail to have a large surface, while maintaining mass as low a as possible.  At best, a solar sail will experience only 9 N of force per square km of sail located in Earth orbit.  Adding the impulse due to incident and reflected photons, it is found that the thrust vector is directed normally to the surface of the sail; hence, by controlling the orientation of the sail relative to the Sun, we can gain or lose orbital angular momentum.A true solar sail has yet to fly; however, significant steps have been taken since the beginning of the 1990s. In February 1993, under the guidance of Vladimir Syromiatnikov, the Russian Space Regatta Consortium deployed a 20 m spinning reflector, Znamya, from a progress supply vehicle.  Observed from the Mir station, this deployment showed that spin deployment could be controlled by passive means.  In May 1996, a large deployable reflector was demonstrated during the shuttle mission STS-77.  The 14 m inflatable antenna experiment was primarily designed as a radio frequency reflector, but the promise of inflatable technology towards solar sailing was clearly demonstrated despite mission anomalies.  In Köln, in December 1999, the German space agency, DLR, in association with ESA deployed a square, 20 m solar sail.  This deployment now forms the basis for a future DLR/ESA in-orbit deployment demonstration perhaps in early 2006. Furthermore, in August 2004, the Institute of Space and Astronautical Science in Japan, ISAS, deployed two solar sails in space from a S-310 sounding rocket. NASA and several private enterprises all seeking to advance solar sailing are also conducting further work. It is thus clear that the technology is currently undergoing a revival in interest and a renewed drive towards flight status.The literature on sailing dynamics, development, and attitude control is wealthy and it is of great interest to sketch some important ideas on these problems. The first solar sail trajectories were calculated by Tsu [5] and London [6]. Tsu investigated various means of propulsion and showed that in many cases solar sails show superior performance when compared to chemical and ion propulsion systems. He used approximated heliocentric equations of motion to obtain spiraling trajectories. London presented similar spiral solutions for Earth-Mars transfers with constant sail orientation using the exact equations of motion. Optimal solar sail trajectories were first computed by Zhukov and Lebedev [7] for interplanetary missions between coplanar circular orbits. In 1980, Jayaraman [8] published similar minimum-time trajectories for transfers between the Earth and Mars. Two years later, Wood et al. [9] presented an analytical proof to show that the orbital transfer times obtained by Jayaraman [8] were incorrect due to the incorrect application of a transversality condition of variational calculus and an erroneous control law. Powers et al. [10] and Powers and Coverstone [11] obtained results similar to those reported in Wood's paper and obtained solutions for transfers to synchronous orbits. The more general time-optimal control problem of three-dimensional, inclined and elliptic departure and rendezvous planet orbits was discussed by Sauer [12]. Hughes and McInnes [13] used genetic algorithms and sequential quadratic programming to obtain interplanetary trajectories via a direct method. Dachwald [14] presented a novel approach based on evolutionary neurocontrollers (ENC) to calculate optimal solar sail trajectories for interplanetary missions.Wie [15] considered the orbital dynamics of a solar sail in the Earth-Sun circular restricted 3-body problem. The equations of motion of the sail are given by a set of nonlinear autonomous ordinary differential equations, which are nonconservative due to the noncentral nature of the force on the sail. They found there are equilibria admitting homoclinic paths, where the stable and unstable invariant manifolds are identical. As well as these periodic orbits about these equilibria also admit homoclinic paths.Waters and McInnes [16] developed an attitude control systems for solar sail spacecraft are presented.  He analysed a sailcraft in an Earth centered elliptic orbit, with particular emphasis on the significant effect of a solar-pressure disturbance torque (caused by an uncertain center-of-mass and center-of-pressure offset).Gong et al. [17] investigated the time-optimal interplanetary transfer trajectories to a circular orbit of given inclination and radius. They derived optimal control law from the principle of maximization. An indirect method is used to solve the optimal control problem by selecting values for the initial adjoint variables. The conditions for the existence of the time-optimal transfer are dependent on the lightness number of the sail and the inclination and radius of the target orbit. A numerical method is used to obtain the boundary values for the time-optimal transfer trajectories. For the cases where no time-optimal transfer trajectories exist, first-order necessary conditions of the optimal control are proposed to obtain feasible solutions. The results show that the transfer time decreases as the minimum distance from the Sun decreases during the transfer duration. For a solar sail with a small lightness number, the transfer time may be evaluated analytically for a three-phase transfer trajectory. The analytical results are compared with previous results and the associated numerical results. The transfer time of the numerical result is smaller than the transfer time from previous results and is larger than the analytical result.In contrast to ordinary ballistic navigation, the mission analysis for solar sail trajectories is not a simple task. The difference is that ballistic navigation is determined by a finite sequence of manoeuvres, with each of them being fully defined by Keplerian theory; 5 orbital parameters (the time they are executed, the intensity of the pulse, and the 3 angles determining its direction). But solar sailing, on the contrary, is a continuous process and the mathematics for reaching final conditions is far more complex, even more than in the case of electric propulsion, since the thrust intensity is related to the sail orientation: the orientation of the force vector applying to a perfectly reflective solar sail is normal to the sail, in the antisolar direction. Its intensity is proportional to the square cosine of the angle between the normal to the sail and the sunline.


## 2. Modeling the Source of Radiation Pressure

The absorbed energy heats the sail, which reradiates that energy from the front and rear surfaces. A solar sail is an oriented surface such that the acceleration experienced by the surface is a function of the surfaces attitude. The variation of the solar radiation pressure with solar distance can be approximated by an inverse square variation relationship; this assumption breaks down at low solar radius, when the finite angular size of the solar disc must be considered.  The modeling of the source of radiation pressure is distinct and independent from the modeling of solar radiation pressure force, which is dependent on the optical properties of the sail surface.Initially, we assume the solar disc has uniform brightness, such that an element of the solar disc will appear equally bright when viewed from any aspect angle. Thus, the specific intensityIr is time independent and isotropic across the solar disc. Therefore, the solar radiation pressure exerted on a radially oriented, perfectly reflecting sail at a heliocentric distance r, can be written as
(1)P(r)=2c∫0∞∫02π∫0θ0Ircos2θsinθdrdθdφ,
where is the velocity of light and θ0 is the angular radius of the solar disc that is given by θ0=sin-1(RSun/r) as shown in Figure 1. The geometric symmetry about the azimuth angle φ leads to P independent of r. Thus, (1) reduces to
(2)P(r)=4π3cI0[1-[1-(RSunr)2]3/2]=2πcI0[RSunr]2+𝒪(RSunr)4,I0=LS4π2RS2,
where I0 is the frequency integrated-specific intensity.Figure 1
Solar radiation pressure due to a finite solar disc.


## 3. The Force Model on the Sailcraft Configurations

### 3.1. The Force due to Solar Radiation Pressure

Solar radiation exerts a pressure on theith surface of area Ai of the solar sail due to photons impinging on a surface in space. If a fraction ρa of the impinging photons is absorbed, a fraction ρs is specularly reflected, and a fraction ρd is diffusely reflected by a surface, then we have  ρa+ρs+ρd=1. The force due to solar radiation pressure (in brief SRP) acting on such a flat Lambertian surface located one astronomical unit (1 AU) from the Sun is modeled as
(3)FSRP=PAi{ρa(e^s·n^)e^s+2ρs(e^s·n^)2n^(e^s+23n^)+ρd(e^s·n^)(e^s+23n^)}=PAi(e^s·n^)×{(ρa+ρd)e^s+(2ρs(e^s·n^)+23ρd)n^},
where P=4.563×10-6 N/m2 is the nominal solar-radiation pressure constant 1 AU from the Sun, A is the surface area, n^ is a unit vector normal to the sail film, and e^s is a unit vector pointing from the Sun to the sail, as shown in Figure 2. The solar radiation pressure varies inversely with the square of the distance from the Sun. The unit vector e^s, as is clear from the geometry depicted in Figure 1, can be written as a linear combination of n^ and a unit vector transverse to it, t^; that is, e^s=cosαn^+cosαt^. Therefore, (3) can be written as
(4)FSRP=PAi{[(1+ρs)cos2α+23ρdcos2α]n^+[23(1-ρs)cosαsinα]t^}.
Now, we can rewrite (4) in the coordinate system e^se^te^w, where e^t is a unit vector perpendicular to e^s, and e^w completes the right-handed coordinate system (STW). Thus, the thrust unit vector (sail normal) in the three-dimensional space is defined by
(5)n^=cosαe^s+sinαcosδe^t+sinαsinδe^w.
Now, the force due solar radiation pressure becomes
(6)FSRP=PAicosα{(ρa+ρs)e^s+(2ρscosα+23ρd)·[cosαe^s+sinαcosδe^t+sinαsinδe^w]23}.Figure 2
Definition of sail control angles.Carrying out the dot products, we get(7)FSRP=PAicosα{(1-ρs+232ρscos2α)e^s+(232ρscosαsinαcosδ+23ρdsinαcosδ)e^t+(232ρscosαsinαsinδ+23ρdsinαsinδ)e^w}.


### 3.2. The Forces due to Coronal Mass Ejections and Solar Wind

Mass is carried away from the Sun into interplanetary space, primarily in the form of protons, through several mechanisms: solar wind, solar flares, and coronal mass ejections. The majority of the solar mass carried away is in the solar wind, a persistent stream of plasma particles generated in the photosphere of the Sun that flows out in all directions. At the Earth, the speed of the solar wind is typically~400 km/s with a particle density of about 10 protons/cm3, Breen et al. [18]. When Mariner 10 flew by Mercury in 1974 it measured a solar wind speed of 423 km/s and a quiescent solar wind particle density of about 60 protons/cm3.Assuming the following simplifying assumptions: (a) spherical envelope for the angular span of ejections, and (b) all particle collisions with the sailcraft are inelastic. The perturbing accelerations on the sailcraft, due to the coronal mass ejections, and due to solar wind respectively are given by Abd El-Salam [19]:
(8)FSW+CME=[(Aim)SWnSWMPvSW2+(Aim)CMEMCMEτCMEvCME2]×(e^s·n^)es=[(Aim)SWnSWMPvSW2+(Aim)CMEMCMEτCMEvCME2]×cosαes,
where (Ai/m)SW, nSW and vSW are the ith area-mass ratio exposed to the solar wind stream, particle density, and the velocity of the solar wind and MP is the mass of the proton, while (Ai/m)CME, τCME, MCME and vCME are the area-mass ratio exposed to the coronal mass ejections, the overall volume, the mass, and the velocity of the coronal mass ejections, respectively. Adding (8) and (7), we obtain the force model due to the solar radiation pressure, solar wind, and the coronal mass ejection:
(9)F=Pcosα{[Ai(1-ρs+2ρscos2α)(Aim)SW+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]es+A(232ρscosαsinαcosδ+23ρdsinαcosδ)et+A(232ρscosαsinαsinδ+23ρdsinαsinδ)ew}.


## 3.1. The Force due to Solar Radiation Pressure

Solar radiation exerts a pressure on theith surface of area Ai of the solar sail due to photons impinging on a surface in space. If a fraction ρa of the impinging photons is absorbed, a fraction ρs is specularly reflected, and a fraction ρd is diffusely reflected by a surface, then we have  ρa+ρs+ρd=1. The force due to solar radiation pressure (in brief SRP) acting on such a flat Lambertian surface located one astronomical unit (1 AU) from the Sun is modeled as
(3)FSRP=PAi{ρa(e^s·n^)e^s+2ρs(e^s·n^)2n^(e^s+23n^)+ρd(e^s·n^)(e^s+23n^)}=PAi(e^s·n^)×{(ρa+ρd)e^s+(2ρs(e^s·n^)+23ρd)n^},
where P=4.563×10-6 N/m2 is the nominal solar-radiation pressure constant 1 AU from the Sun, A is the surface area, n^ is a unit vector normal to the sail film, and e^s is a unit vector pointing from the Sun to the sail, as shown in Figure 2. The solar radiation pressure varies inversely with the square of the distance from the Sun. The unit vector e^s, as is clear from the geometry depicted in Figure 1, can be written as a linear combination of n^ and a unit vector transverse to it, t^; that is, e^s=cosαn^+cosαt^. Therefore, (3) can be written as
(4)FSRP=PAi{[(1+ρs)cos2α+23ρdcos2α]n^+[23(1-ρs)cosαsinα]t^}.
Now, we can rewrite (4) in the coordinate system e^se^te^w, where e^t is a unit vector perpendicular to e^s, and e^w completes the right-handed coordinate system (STW). Thus, the thrust unit vector (sail normal) in the three-dimensional space is defined by
(5)n^=cosαe^s+sinαcosδe^t+sinαsinδe^w.
Now, the force due solar radiation pressure becomes
(6)FSRP=PAicosα{(ρa+ρs)e^s+(2ρscosα+23ρd)·[cosαe^s+sinαcosδe^t+sinαsinδe^w]23}.Figure 2
Definition of sail control angles.Carrying out the dot products, we get(7)FSRP=PAicosα{(1-ρs+232ρscos2α)e^s+(232ρscosαsinαcosδ+23ρdsinαcosδ)e^t+(232ρscosαsinαsinδ+23ρdsinαsinδ)e^w}.


## 3.2. The Forces due to Coronal Mass Ejections and Solar Wind

Mass is carried away from the Sun into interplanetary space, primarily in the form of protons, through several mechanisms: solar wind, solar flares, and coronal mass ejections. The majority of the solar mass carried away is in the solar wind, a persistent stream of plasma particles generated in the photosphere of the Sun that flows out in all directions. At the Earth, the speed of the solar wind is typically~400 km/s with a particle density of about 10 protons/cm3, Breen et al. [18]. When Mariner 10 flew by Mercury in 1974 it measured a solar wind speed of 423 km/s and a quiescent solar wind particle density of about 60 protons/cm3.Assuming the following simplifying assumptions: (a) spherical envelope for the angular span of ejections, and (b) all particle collisions with the sailcraft are inelastic. The perturbing accelerations on the sailcraft, due to the coronal mass ejections, and due to solar wind respectively are given by Abd El-Salam [19]:
(8)FSW+CME=[(Aim)SWnSWMPvSW2+(Aim)CMEMCMEτCMEvCME2]×(e^s·n^)es=[(Aim)SWnSWMPvSW2+(Aim)CMEMCMEτCMEvCME2]×cosαes,
where (Ai/m)SW, nSW and vSW are the ith area-mass ratio exposed to the solar wind stream, particle density, and the velocity of the solar wind and MP is the mass of the proton, while (Ai/m)CME, τCME, MCME and vCME are the area-mass ratio exposed to the coronal mass ejections, the overall volume, the mass, and the velocity of the coronal mass ejections, respectively. Adding (8) and (7), we obtain the force model due to the solar radiation pressure, solar wind, and the coronal mass ejection:
(9)F=Pcosα{[Ai(1-ρs+2ρscos2α)(Aim)SW+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]es+A(232ρscosαsinαcosδ+23ρdsinαcosδ)et+A(232ρscosαsinαsinδ+23ρdsinαsinδ)ew}.


## 4. Optimum Force Vector

Referring to Figure2, we see the orientation of the solar sail, and so the thrust force vector is described relative to the sun-line by the sail pitch angle, α, and clock angle, δ. We note that the sail control angles αandδ can be defined as either
(10)0∘≤α≤90∘,0∘≤δ≤360∘or-90∘≤α≤0∘,0∘≤δ≤360∘.
Now, the sail thrust vector is defined by the cone and clock angles in the radial-transverse-normal (STW) frame. In order to optimize the sail control angles, we define a required direction, ξ^, along which the component of the sail thrust is to be maximized:
(11)ξ^=cosα~e^s+sinα~cosδ~e^t+sinα~sinδ~e^w.
The force in this required direction, namely, ξ^-direction is
(12)Fξ=F·ξ^=Pcosα{[(Aim)CME(MCMEτCME)vCME2)Ai(1-ρs+2ρscos2α)(Aim)SW+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)](es·ξ^)+Ai(2ρscosαsinαcosδ+23ρdsinαcosδ)×(et·ξ^)+Ai(2ρscosαsinαsinδ+23ρdsinαsinδ)×(ew·ξ^)[Ai(1-ρs+2ρscos2α)(Aim)SW}.Using (11), we can obtain
(13)Fξ=Pcosα{[Ai(1-ρs+2ρscos2α)(Aim)SW+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]cosα~+Ai(232ρscosαsinαcosδ+23ρdsinαcosδ)sinα~cosδ~+Ai(232ρscosαsinαsinδ+23ρdsinαsinδ)sinα~sinδ~},Fξ=P{23Ai(1-ρs)cosαcosα~+2ρscos3αcosα~+[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]×cosα~cosα+Ai(232ρscos2αsinαcosδsinα~cosδ~+23ρdcosαsinαcosδsinα~cosδ~)+Ai(232ρscos2αsinαsinδsinα~sinδ~+23ρdcosαsinαsinδsinα~sinδ~)}.
Arranging the terms and using the trigonometric identities, we get
(14)Fξ=PAi(1-ρs)cosαcosα~+2PρsAi×(cos3αcosα~+cos2αsinαsinα~cos(δ-δ~))+P[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]×cosα~cosα+23PAiρdcosαsinαsinα~cos(δ-δ~).
Consider the following definitions:
(15)𝒜1=AiP(1-ρs)+[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2],𝒜2=2ρsAiP,𝒜3=23ρdAiP.
Equation (14) can be simplified to
(16)Fξ=𝒜1cosαcosα~+𝒜2(cosαcosα~+sinαsinα~cos(δ-δ~))cos2α+𝒜3cosαsinαsinα~cos(δ-δ~).
Differentiating (16) with respect to the cone angle and finding the turning points, (∂Fξ/∂α), then we can obtain the optimal sail cone angle, which maximizes the force in the required direction as
(17)∂Fξ∂α=0=-𝒜1sinαcosα~+𝒜2cosα~×(-3cos2αsinα+[cosα-3sin2αcosα]tanα~cos(δ-δ~))+𝒜3[cos2α-sin2α]sinα~cos(δ-δ~)=0.
Dividing on 3cos2αsinα yields
(18)-𝒜13(1+tan2α)cosα~+𝒜2×cosα~(-1+tanα~cos(δ-δ~)3cosαsinα-tanαtanα~cos(δ-δ~)tanα~cos(δ-δ~)3cosαsinα)+𝒜331sinα[(1-tan2α)]sinα~cos(δ-δ~)=0,-𝒜13cosα~-𝒜13tan2αcosα~-𝒜2cosα~+tanα~cos(δ-δ~)3cosαsinα𝒜2cosα~-tanαtanα~cos(δ-δ~)𝒜2cosα~+𝒜331sinαsinα~cos(δ-δ~)-𝒜331sinαtan2αsinα~cos(δ-δ~)=0,
which can be simplified to
(19)(13cosαsinα)×{tanα~cos(δ-δ~)𝒜2cosα~-𝒜1cosα~cosαsinα-3𝒜2cosα~cosαsinα-𝒜1tan2αcosα~cosαsinα+𝒜3sinα~cos(δ-δ~)cosα-𝒜3tan2αsinαsinα~cos(δ-δ~)-3𝒜2tanαtanα~cos(δ-δ~)cosα~cosαsinα}=0.
Equating the numerator to zero and rearranging the terms yield
(20)tanα~cos(δ-δ~)𝒜2cosα~-(𝒜1+3𝒜2)cosα~cosαsinα-𝒜1tan2αcosα~cosαsinα+𝒜3sinα~cos(δ-δ~)cosα-𝒜3tan2αsinαsinα~cos(δ-δ~)-3𝒜2tanαtanα~cos(δ-δ~)cosα~cosαsinα=0.
Multiplying both sides by  1/cos2α=1+tan2α  yields
(21)𝒜2tanα~cos(δ-δ~)cosα~(1+tan2α)-(𝒜1+3𝒜2)cosα~tanα-𝒜1cosα~tan3α+𝒜3sinα~cos(δ-δ~)(1+tan2α)cosα-𝒜3sinα~cos(δ-δ~)(1+tan2α)tanαsinα-3𝒜2tanα~cos(δ-δ~)cosα~tan2α=0,𝒜2tanα~cos(δ-δ~)cosα~+𝒜2tanα~cos(δ-δ~)cosα~tan2α-(𝒜1+3𝒜2)cosα~tanα-𝒜1cosα~tan3α+𝒜3sinα~cos(δ-δ~)cosα+𝒜3sinα~cos(δ-δ~)tan2αcosα-𝒜3sinα~cos(δ-δ~)tanαsinα-𝒜3sinα~cos(δ-δ~)tan3αsinα-3𝒜2tanα~cos(δ-δ~)cosα~tan2α=0.
Rearranging the terms yields
(22)-𝒜3sinα~cos(δ-δ~)(tanαsinα-cosα)×[tan2α+1]-𝒜1cosα~tan3α-2𝒜2tanα~cos(δ-δ~)cosα~tan2α-(𝒜1+3𝒜2)cosα~tanα+𝒜2tanα~cos(δ-δ~)cosα~=0.
Dividing on cosα~ yields
(23)-𝒜3tanα~cos(δ-δ~)(tanαsinα-cosα)×[tan2α+1]-𝒜1tan3α-2𝒜2tanα~cos(δ-δ~)tan2α-(𝒜1+3𝒜2)tanα+𝒜2tanα~cos(δ-δ~)=0
This is a transcendental equation which is so difficult to be solved analytically, but in the following, we will consider some special cases.


### 4.1. Special Case I

Ignoring the absorbed and diffusely reflected parts of the radiation, of the coefficients (15) become
(24)𝒜1=[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]𝒜2=2ρsA,𝒜3=0-𝒜1tan3α-2𝒜2tanα~cos(δ-δ~)tan2α-(𝒜1+3𝒜2)tanα+𝒜2tanα~cos(δ-δ~)=0.
Dividing on -𝒜1  yields
(25)tan3α+2𝒜2𝒜1tanα~cos(δ-δ~)tan2α+(𝒜1+3𝒜2)𝒜1tanα-𝒜2𝒜1tanα~cos(δ-δ~)=0.
Let us now define these new coefficients:
(26)ℬ1=2𝒜2𝒜1tanα~cos(δ-δ~),ℬ2=(𝒜1+3𝒜2)𝒜1,ℬ3=-𝒜2𝒜1tanα~cos(δ-δ~).
Now, the cubic equation (25) becomes
(27)tan3α+ℬ1tan2α+ℬ2tanα+ℬ3=0,
which has the solution in the maximizing direction ξ(28)(tanα1*)ξ=-13ℬ1+𝔹+ℬ12-3ℬ29𝔹,(29)(tanα2,3*)ξ=-13ℬ1-12𝔹-ℬ12-3ℬ218𝔹±i32ℬ1[𝔹-ℬ12-3ℬ29𝔹],
where
(30)𝔹=(-16ℬ1ℬ2ℬ3+127ℬ23+14ℬ22)1/216ℬ1ℬ2-12ℬ3-127ℬ13+(127ℬ13ℬ3-1108ℬ12ℬ22-16ℬ1ℬ2ℬ3+127ℬ23+14ℬ22)1/2)1/3.
Since tanα2,3* given by (29) are imaginary solutions, therefore, these two roots are rejected.   And thus the only considered root is the first root  tanα1*.


### 4.2. Special Case II

In addition to the simplification addressed in the special Case I, If we consider also the solar wind and the coronal mass ejection does not depend on the cone angleα. Then, (27) can be reduced to the following quadratic equation:
(31)-2𝒜2tanα~cos(δ-δ~)tan2α-3𝒜2tanα+𝒜2tanα~cos(δ-δ~)=0.
The solution of this equation yields the required cone angle α=α* to maximize the force vector in the required direction:
(32)(tanα*)ξ=-3±9+8tan2α~cos2(δ-δ~)4tanα~cos(δ-δ~).


## 4.1. Special Case I

Ignoring the absorbed and diffusely reflected parts of the radiation, of the coefficients (15) become
(24)𝒜1=[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]𝒜2=2ρsA,𝒜3=0-𝒜1tan3α-2𝒜2tanα~cos(δ-δ~)tan2α-(𝒜1+3𝒜2)tanα+𝒜2tanα~cos(δ-δ~)=0.
Dividing on -𝒜1  yields
(25)tan3α+2𝒜2𝒜1tanα~cos(δ-δ~)tan2α+(𝒜1+3𝒜2)𝒜1tanα-𝒜2𝒜1tanα~cos(δ-δ~)=0.
Let us now define these new coefficients:
(26)ℬ1=2𝒜2𝒜1tanα~cos(δ-δ~),ℬ2=(𝒜1+3𝒜2)𝒜1,ℬ3=-𝒜2𝒜1tanα~cos(δ-δ~).
Now, the cubic equation (25) becomes
(27)tan3α+ℬ1tan2α+ℬ2tanα+ℬ3=0,
which has the solution in the maximizing direction ξ(28)(tanα1*)ξ=-13ℬ1+𝔹+ℬ12-3ℬ29𝔹,(29)(tanα2,3*)ξ=-13ℬ1-12𝔹-ℬ12-3ℬ218𝔹±i32ℬ1[𝔹-ℬ12-3ℬ29𝔹],
where
(30)𝔹=(-16ℬ1ℬ2ℬ3+127ℬ23+14ℬ22)1/216ℬ1ℬ2-12ℬ3-127ℬ13+(127ℬ13ℬ3-1108ℬ12ℬ22-16ℬ1ℬ2ℬ3+127ℬ23+14ℬ22)1/2)1/3.
Since tanα2,3* given by (29) are imaginary solutions, therefore, these two roots are rejected.   And thus the only considered root is the first root  tanα1*.


## 4.2. Special Case II

In addition to the simplification addressed in the special Case I, If we consider also the solar wind and the coronal mass ejection does not depend on the cone angleα. Then, (27) can be reduced to the following quadratic equation:
(31)-2𝒜2tanα~cos(δ-δ~)tan2α-3𝒜2tanα+𝒜2tanα~cos(δ-δ~)=0.
The solution of this equation yields the required cone angle α=α* to maximize the force vector in the required direction:
(32)(tanα*)ξ=-3±9+8tan2α~cos2(δ-δ~)4tanα~cos(δ-δ~).


## 5. Heliocentric Orbital Dynamics

In the first instance, the Sun-centered orbital dynamics of solar sail spacecraft will be described, since this is traditionally where solar sailing is applied and where the bulk of the mission concepts exists, due to the largeΔv requirements.  The vector equation of motion of a solar sail spacecraft moving in a heliocentric orbit is defined by (33), where r is the position vector of the spacecraft with respect to the Sun at time t.  The gravitational parameter of the Sun is defined by μ.  For an ideal sail, the thrust vector is aligned along the sail normal direction, n:
(33)d2rdt2+μr2r^=P{23[(Aim)CMEAicosα(1-ρs+2ρscos2α)+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]e^s+Aicosα(232ρscosαsinαcosδ+23ρdsinαcosδ)e^t+Aicosα(232ρscosαsinαsinδ+23ρdsinαsinδ)e^w}.
For the simple analysis to follow, it is appropriate to show how this equation is represented in two-dimensional polar coordinates, where θ is the azimuth angle from the x-axis, r is the distance of the spacecraft from the Sun.  The pitch or cone angle, α, is defined as the angle between the sail normal and the radial vector.Consider(34)r¨-rθ˙2-rφ˙2sin2θ=-μr2+βμr2PAi×cosα[(Aim)CME(1-ρs+2ρscos2α)+(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2],1r(r2θ˙)·-rφ˙2sinθcosθ=βμr2PAcosα(232ρscosαsinαcosδ+23ρdsinαcosδ),1rsinθ(r2φ˙sin2θ)·=βμr2Aicosα(2ρscosαsinαsinδ+23ρdsinαsinδ).
The parameter, β, is known as the sail lightness number and is defined as the ratio of the local solar radiation pressure acceleration produced by the sail to the local solar gravitational acceleration.  This number is a useful parameter since it is independent of the solar distance and defines the performance of the sail.  Another useful performance parameter is the sail characteristic acceleration, ac, which is the solar radiation pressure-induced acceleration of an ideal sail, pitched face on to the Sun (α=0) at 1 AU from the Sun.  It is easy to calculate that a characteristic acceleration of 1.0 mm s-2 corresponds to a dimensionless sail lightness number of 0.1686.  As has been discussed, the solar sail characteristic acceleration is dependent on the surface area and reflectivity of the sail assembly and the mass of the spacecraft.


## 6. Lagrange’s Planetary Equations

The state of a spacecraft can be described by a vector of 6 orbital elements, namely, semimajor axis,a, eccentricity, e, inclination, i, argument of perihelion, ω, right ascension of the ascending node, Ω, true anomaly, f, or any other time elements. Theses 6 elements are equivalent to 6 Cartesian position and velocity components. To measure the rate of change in these elements, we use the very famous system of 6 first-order differential equations known as Lagrange planetary equations; see (35). They are used to solve the equations motion of the sailcraft. These equations are particularly useful when we want to maximize the rate of change of a particular orbital element.  This is useful when we want to modify one orbital element, while leaving the other time-averaged elements unchanged.  One form of this system is given by
(35)dadf=2pr2μ(1-e2)2[Sesinf+Tpr],dedf=r2μ[Ssinf+T(1+rp)cosf+Trpe],didf=r3μpcos(f+ω)W,dΩdf=r3μpsinisin(f+ω)W,dωdf=-dΩdfcosi+r2μe[-Scosf+T(1+rp)sinf],dtdf=r2μp[1-r2μe[Scosf-T(1+rp)sinf]],
where a,e,i,Ω,ω,andf are the usual Keplerian orbital elements, semimajor axis, eccentricity, inclination, longitude of the ascending node, and true anomaly, respectively, p is the semilatus rectum, μ the gravitational parameter, and n is the orbital mean motion.  The components of the solar sail thrust are denoted by S,T,W radial, transverse, and normal, respectively.All Lagrange planetary equations can be written in compact form as only one equation as follows:(36)d𝒳df=F·Ξ𝒳,
where  𝒳(≡a,e,i,Ω,ω,f) denotes any element of the Keplerian orbital elements.F,Ξ𝒳  are the force vector and primer vector of optimization, respectively. These vectors can be conveniently written as
(37)F=Se^s+Te^t+We^w,Ξ𝒳=Ξs𝒳e^s+Ξt𝒳e^t+Ξw𝒳e^w,
where
(38)S=βμr2PAi×cosα[(1-ρs+2ρscos2α)+(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2],T=βμr2PAi×cosα(2ρscosαsinαcosδ+23ρdsinαcosδ),W=βμr2Ai×cosα(2ρscosαsinαsinδ+23ρdsinαsinδ).The primer vector in the direction of maximizing the semimajor axis is(39)Ξa=Ξsae^s+Ξtae^t=2pr2μ(1-e2)2[esinfe^s+pre^t].
The primer vector in the direction of maximizing the eccentricity is
(40)Ξe=Ξsee^s+Ξtee^t=r2μ[sinfe^s+(1+rp(e+cosf))e^t].
The primer vector in the direction of maximizing the inclination is
(41)Ξi=Ξwie^w=r3μpcos(f+ω)e^w.
The primer vector in the direction of maximizing the longitude of the ascending node is
(42)ΞΩ=ΞwΩe^w=r3μpsinisin(f+ω)e^w.
The primer vector in the direction of maximizing the argument of periapsis is
(43)Ξω=Ξsωe^s+Ξtωe^t+Ξwωe^w=r2μe[-erpsinisin(f+ω)cosie^s-cosfe^t+(1+rp)sinfe^w].
Using this system, we can maximize directly any orbital elements or any other dynamical orbital parameters such as, the radius of periapsis and apoapsis.  For example, if the radius of apoapsis, rA=a(1+e) is differentiated, we obtain
(44)drAdf=dadf(1+e)+adedf.


## 7. Some New Locally Optimal Control Laws

To maximize the thrust to obtain certain required maximization in some particular orbital element, we set(∂/∂α)(d𝒳/df)=0 and thus obtain a new set of locally optimal control laws. These control laws cannot guarantee global optimality, and they are often termed closed-loop methods.  Global optimality requires the use of numerical methods, and even then, the true optimum solution is hard to attain.


### 7.1. Optimal Control Law for the Rate of Change of Semimajor Axis

The primer vector components maximize the first Lagrange planetary equation for the rate of change of semimajor axis.In what follows, we will find an optimum cone angle that maximizes the rate of change of semimajor axis:(45)∂∂α(dadf)=2pr2μ(1-e2)2[esinf(∂S∂α)+pr(∂T∂α)]=0⟹[esinf(∂S∂α)+pr(∂T∂α)]=0.
The force in this required direction is given by
(46)dadt=F·Ξa=P2pr2μ(1-e2)2×cosα{[(Aim)CME(MCMEτCME)vCME2)]Ai(1-ρs+2ρscos2α)+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×(es·Ξa)+Ai(2ρscosαsinαcosδ+23ρdsinαcosδ)×(et·Ξa)+Ai(2ρscosαsinαsinδ+23ρdsinαsinδ)×(ew·Ξa)[(Aim)CME(MCMEτCME)vCME2)]Ai(1-ρs+2ρscos2α)},(47)dadt=P2pr2μ(1-e2)2×cosα{(et·pre^t)[(Aim)CMEAi(1-ρs+2ρscos2α)+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×(es·esinfe^s)+Ai(2ρscosαsinαcosδ+23ρdsinαcosδ)×(et·pre^t)},dadt=P2pr2μ(1-e2)2×cosα{[(Aim)CMEAi(1-ρs+2ρscos2α)+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]esinf+Aipr(232ρscosαsinαcosδ+23ρdsinαcosδ)}.
Consider the following definitions:
(48)𝒜1=AiP(1-ρs)+[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]𝒜2=2ρsAiP,𝒜3=23ρdAiP.
Equation (46) can be simplified to
(49)dadt=2pr2μ(1-e2)2×[{𝒜3prcosαsinαcosδ}𝒜1esinfcosα+𝒜2(esinfcosα+prsinαcosδ)cos2α+𝒜3prcosαsinαcosδ].
If we differentiate (49) with respect to the cone angle and find the turning points, (∂Fξ/∂α), then we can obtain the optimal sail cone angle, which maximizes the force in the required direction as follows:
(50)∂∂α(dadt)=-𝒜1esinfsinα+𝒜2e×sinf(-3cos2αsinα+pr×[cosα-3sin2αcosα]cosδpr)+𝒜3pr[cos2α-sin2α]cosδ=0.
Dividing on 3cos2αsinα yields
(51)∂∂α(dadt)=-𝒜13(1+tan2α)esinf+𝒜2e×sinf(-1+pcosδ3rcosαsinα-ptanαcosδr)+𝒜3p3rsinα[1-tan2α]cosδ=0=-𝒜13esinf-𝒜13esinftan2α-𝒜2esinf+𝒜2esinfpcosδ3rcosαsinα-𝒜2pesinftanαcosδr+𝒜3p3rsinαcosδ-𝒜3p3rsinαtan2αcosδ=0,
which can be written as
(52)∂∂α(dadt)=(13cosαsinα)×{𝒜2eprsinfcosδ-𝒜1esinfcosαsinα-𝒜1esinftan2αcosαsinα-3𝒜2esinfcosαsinα-3𝒜2presinfcosδtanαcosαsinα+𝒜3prcosδcosα-𝒜3prcosδtan2αcosα}=0.
Equating the numerator to zero and rearranging the terms yield
(53)∂∂α(dadt)=𝒜2eprsinfcosδ-(𝒜1+3𝒜2)esinfcosαsinα-𝒜1esinftan2αcosαsinα-3𝒜2presinfcosδtanαcosαsinα+𝒜3prcosδcosα-𝒜3prcosδtan2αcosα=0.
Multiplying both sides by  1/cos2α=1+tan2α  yields
(54)∂∂α(dadt)=𝒜2eprsinfcosδ(1+tan2α)-(𝒜1+3𝒜2)esinftanα-𝒜1esinftan3α-3𝒜2presinfcosδtan2α+𝒜3prcosδ(1+tan2α)cosα-𝒜3prcosδtanα(1+tan2α)sinα=0.
Finally, we can obtain
(55)=-𝒜1tan3-2𝒜2prcosδtan2α+𝒜3presinfcosδ(1+tan2α)×[cosα-tanαsinα]-(𝒜1+3𝒜2)tanα+𝒜2prcosδ=0.
This is a transcendental equation which is so difficult to be solved analytically, but in the following, we will consider some special cases.


### 7.2. Special Case I

Ignoring the absorbed and diffusely reflected parts of the radiation, of the coefficients (47) become
(56)𝒜1=[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]𝒜2=2ρsAi,𝒜3=0.
Thus, (55) becomes
(57)∂∂α(dadt)=tan3α+2𝒜2𝒜1prcosδtan2α+(1+3𝒜2𝒜1)tanα-𝒜2𝒜1prcosδ=0.
Let us now define the new coefficients:
(58)𝒦1=2𝒜2𝒜1prcosδ,𝒦2=(1+3𝒜2𝒜1),𝒦3=-𝒜2𝒜1prcosδ.
Equation (57) becomes clear as a cubic equation:
(59)tan3α+𝒦1tan2α+𝒦2tanα+𝒦3=0,
which has the same solution as given by (27), but upon replacing the coefficients  ℬ's with the newly defined coefficients 𝒦's we obtain
(60)(tanα1*)a=-13𝒦1+𝕂+𝒦12-3𝒦29𝕂,(tanα2,3*)a=-13𝒦1-12𝕂-𝒦12-3𝒦218𝕂±i32𝒦1[𝕂-𝒦12-3𝒦29𝕂],
where
(61)𝕂=(16𝒦1𝒦2-12𝒦3-127𝒦13+(127𝒦13𝒦3-1108𝒦12𝒦22-16𝒦1𝒦2𝒦3+127𝒦23+14𝒦22)1/2)1/3.
Since tanα2,3* are imaginary solution, therefore, these two roots are rejected.   And thus the only considered root is the first root  tanα1*.


### 7.3. Special Case II

Ignoring the absorbed and diffusely reflected parts of the radiation and assuming that the solar wind and the coronal mass ejection do not depend on the cone angle, then (57) can be reduced to the following quadratic equation:
(62)dadt=F·Ξa=P2pr2μ(1-e2)2×{[(Aim)CME2ρsAicosαcos2α+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)](es·Ξa)+Aicosα(2ρscosαsinαcosδ)(et·Ξa)+Aicosα(2ρscosαsinαsinδ)×(ew·Ξa)(Aim)CME},dadt=P2pr2μ(1-e2)2×{[(Aim)CME2ρsAicosαcos2α+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×(es·esinfe^s)+Aicosα(2ρscosαsinαcosδ)×(et·pre^t)(Aim)CME},
which can be simplified to
(63)dadt=P2pr2μ(1-e2)2×{[(Aim)CME2ρsAicosαcos2α+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×esinf+Aiprcosα(2ρscosαsinαcosδ)}.
Consider the following definitions:
(64)𝒜0=(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2,𝒜2=2ρsAiP.
Equation (63) can be simplified to
(65)dadt=2pr2μ(1-e2)2×[pr𝒜0esinf+𝒜2(esinfcosα+prsinαcosδ)cos2α].
Differentiating  (65) with respect to α  and equating the result to zero, we obtain
(66)∂∂α(dadt)=-3esinfcos2αsinα+prcosαcosδ-3prsin2αcosαcosδ=0.
Dividing on 3cos2αsinα yields
(67)∂∂α(dadt)=-esinf+13prsinαcosαcosδ-3prtanαcosδ=0,
which can be written as
(68)(13cosαsinα){-3esinfcosαsinα+prcosδ-3prtanαcosαsinαcosδ}=0.
Equating the numerator to zero and rearranging the terms yield
(69)2(1+ecosfesinf)cosδtan2α+3tanα-1+ecosfesinfcosδ=0,
which has the solution of quadratic equations
(70)(tanα*)a=-3±9+8((1+ecosf)/esinf)2cos2δ4((1+ecosf)/esinf)cosδ.
This equation represents the control law that maximizes the rate of change of semimajor axis. It computes the actual sail pitch/cone angle profile necessary to increase the semimajor axis at a maximum rate (Figure 3).(a) The optimum cone angle that maximizes the rate of change of semimajor axis when taking the positive sign in the numerator. (b) The optimum cone angle that maximizes the rate of change of semimajor axis when taking the negative sign in the numerator.
(a)
(b)


### 7.4. Optimal Control Law for the Rate of Change of Eccentricity

When comparing (70) and (69), we observe that the required direction to maximize certain variable is given by tanα~=(1+ecosf)/esinf=(p/r)/esinf, which represents  the transverse component divided by the radial component of the primer vector. Thus, in case of eccentricity, set
(71)tanα~=ΞteΞse=cosf+(r/p)(e+cosf)sinf=(2+ecosf)cotf(1+ecosf)+ecscf(1+ecosf).
The maximum rate of change of the orbital eccentricity follows directly, taking into account the simplification addressed in special Case II (Figure 4), as
(72)(tanα*)e=(-3±9+8((2+ecosf)cotf+ecscf(1+ecosf))2cos2δ)×(4((2+ecosf)cotf+ecscf(1+ecosf))cosδ)-1.(a) The optimum cone angle that maximizes the rate of change of eccentricity when taking the positive sign in the numerator. (b) The optimum cone angle that maximizes the rate of change of eccentricity when taking the negative sign in the numerator.
(a)
(b)


### 7.5. Optimal Control Law for the Rate of Change of Inclination

Changes in the out-of-plane orbital elements, such as inclination and right ascension of the ascending node, can be effected by the use of simple switching functions.  For maximum rate of change of inclination, the solar sail thrust can be directed alternately above and below the orbit plane every half orbit by thesign function which has +1 or −1.In what follows, we will find an optimum cone angle that maximizes the rate of change of inclination:(73)∂∂α(didf)=r3μpsign[cos(f+ω)]∂∂αW=0,
from which one can obtain
(74)3ρsρd[1-2tan2α]=1cosα[tan2α-1].
Setting x=tan2α,  the above equation can be written as
(75)x3+(-1-ρd3ρs)x2+(2ρd-3ρs3ρs)x+(1-ρd3ρs)=0,
which has the solution
(76)(tan2α1*)i=19ρs(1+ρd)+(ρ~0+ρ~11/2)1/3+19ρs2[(1+ρd)2-9ρs(2ρd-3ρs)]+181ρs[16(-1-ρd)(2ρd-3ρs)-12(1-ρd)-127(-1-ρd)3](ρ~0+ρ~11/2)-1/3,(tan2α2,3*)i=19ρs(1+ρd)-12(ρ~0+ρ~11/2)1/3-19×18ρs2[(1+ρd)2-9ρs(2ρd-3ρs)]×(ρ~0+ρ~11/2)-1/3±i36ρs(1+ρd)×[(ρ~0+ρ~11/2)1/3-181ρs2×[(1+ρd)2-9ρs(2ρd-3ρs)]×(ρ~0+ρ~11/2)-1/3181ρs2],
where
(77)ρ~0=145ρs2(1+ρd)(2ρd-3ρs)-16ρs(1-ρd)+127×27ρs3(1+ρd)3,ρ~1=1729ρs4(1+ρd)3(9-9ρd)-1162ρs(1-ρd2)(2ρd-3ρs)+1729ρs3(2ρd-3ρs)3+136ρs2(2ρd-3ρs)2-18748(1+ρd)2(1-ρd)2.
As a special case when ignoring the diffusive part of radiation, we obtain
(78)3ρs[1-2tan2α]=0,⟹(tanα*)i=12⟹αi*=35.26°.
This result is in agreement with the previous results.  If we substitute the control law into the inclination equation in Lagrange planetary equations and integrate over one orbit, we can obtain the change in inclination per orbit.  We find that the change in inclination per orbit is independent of orbit radius and only depends on the sail lightness number  (straight number relation) as follows:
(79)Δi=88.2β(degreesperorbit).
However, closer orbits to the Sun have shorter orbit periods, and so the time to achieve an overall inclination change is shorter. This optimal angle is 35.26°, which enables us to maximize angular momentum.


### 7.6. Optimal Control Law for the Rate of Change of Ascending Node

The procedure is the same as in optimal control law for inclination(80)∂∂α(dΩdf)=r3μpsinisign[sin(f+ω)]∂∂αW=0⟹3ρsρd[1-2tan2α]=1cosα[tan2α-1],
which has the same solutions as those maximizing the inclination and is in good agreement with the previous results. If we substitute this control law into the ascending node in Lagrange planetary equations and integrate over one orbit, we can obtain the change in ascending node per orbit. We find that the change in ascending node per orbit is independent of orbit radius and only depends on the sail lightness number as follows:
(81)ΔΩ=(88.2βsini)(degreesperorbit).
However, closer orbits to the Sun have shorter orbit periods, and so the time to achieve an overall inclination change is shorter (Figure 5).Figure 5
The optimum change in the longitude in the ascending nodeΔΩ versus the inclination sin i and sail lightness number β.


### 7.7. Optimal Control Law for the Rate of Change of Argument of Perihelion

When comparing (70) and (69), we observe that the required direction to maximize a certain variable is given by tanα~=(1+ecosf)/esinf=(p/r)/esinf, which represents  the transverse component divided by the radial component of the primer vector. Thus, in case of argument of perihelion (Figure 6), we obtain
(82)tanα~=ΞtωΞsω=(1+(r/p))sinfcosf=(1+(rp))tanf=(2+ecosf1+ecosf)tanf,(tanα*)ω=-3±9+8(((2+ecosf)/(1+ecosf))tanf)2cos2δ4(((2+ecosf)/(1+ecosf))tanf)cosδ.(a) The optimum cone angle that maximizes the rate of change of argument of perihelion when taking the positive sign in the numerator. (b) The optimum cone angle that maximizes the rate of change of argument of perihelion when taking the negative sign in the numerator.
(a)
(b)


## 7.1. Optimal Control Law for the Rate of Change of Semimajor Axis

The primer vector components maximize the first Lagrange planetary equation for the rate of change of semimajor axis.In what follows, we will find an optimum cone angle that maximizes the rate of change of semimajor axis:(45)∂∂α(dadf)=2pr2μ(1-e2)2[esinf(∂S∂α)+pr(∂T∂α)]=0⟹[esinf(∂S∂α)+pr(∂T∂α)]=0.
The force in this required direction is given by
(46)dadt=F·Ξa=P2pr2μ(1-e2)2×cosα{[(Aim)CME(MCMEτCME)vCME2)]Ai(1-ρs+2ρscos2α)+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×(es·Ξa)+Ai(2ρscosαsinαcosδ+23ρdsinαcosδ)×(et·Ξa)+Ai(2ρscosαsinαsinδ+23ρdsinαsinδ)×(ew·Ξa)[(Aim)CME(MCMEτCME)vCME2)]Ai(1-ρs+2ρscos2α)},(47)dadt=P2pr2μ(1-e2)2×cosα{(et·pre^t)[(Aim)CMEAi(1-ρs+2ρscos2α)+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×(es·esinfe^s)+Ai(2ρscosαsinαcosδ+23ρdsinαcosδ)×(et·pre^t)},dadt=P2pr2μ(1-e2)2×cosα{[(Aim)CMEAi(1-ρs+2ρscos2α)+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]esinf+Aipr(232ρscosαsinαcosδ+23ρdsinαcosδ)}.
Consider the following definitions:
(48)𝒜1=AiP(1-ρs)+[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]𝒜2=2ρsAiP,𝒜3=23ρdAiP.
Equation (46) can be simplified to
(49)dadt=2pr2μ(1-e2)2×[{𝒜3prcosαsinαcosδ}𝒜1esinfcosα+𝒜2(esinfcosα+prsinαcosδ)cos2α+𝒜3prcosαsinαcosδ].
If we differentiate (49) with respect to the cone angle and find the turning points, (∂Fξ/∂α), then we can obtain the optimal sail cone angle, which maximizes the force in the required direction as follows:
(50)∂∂α(dadt)=-𝒜1esinfsinα+𝒜2e×sinf(-3cos2αsinα+pr×[cosα-3sin2αcosα]cosδpr)+𝒜3pr[cos2α-sin2α]cosδ=0.
Dividing on 3cos2αsinα yields
(51)∂∂α(dadt)=-𝒜13(1+tan2α)esinf+𝒜2e×sinf(-1+pcosδ3rcosαsinα-ptanαcosδr)+𝒜3p3rsinα[1-tan2α]cosδ=0=-𝒜13esinf-𝒜13esinftan2α-𝒜2esinf+𝒜2esinfpcosδ3rcosαsinα-𝒜2pesinftanαcosδr+𝒜3p3rsinαcosδ-𝒜3p3rsinαtan2αcosδ=0,
which can be written as
(52)∂∂α(dadt)=(13cosαsinα)×{𝒜2eprsinfcosδ-𝒜1esinfcosαsinα-𝒜1esinftan2αcosαsinα-3𝒜2esinfcosαsinα-3𝒜2presinfcosδtanαcosαsinα+𝒜3prcosδcosα-𝒜3prcosδtan2αcosα}=0.
Equating the numerator to zero and rearranging the terms yield
(53)∂∂α(dadt)=𝒜2eprsinfcosδ-(𝒜1+3𝒜2)esinfcosαsinα-𝒜1esinftan2αcosαsinα-3𝒜2presinfcosδtanαcosαsinα+𝒜3prcosδcosα-𝒜3prcosδtan2αcosα=0.
Multiplying both sides by  1/cos2α=1+tan2α  yields
(54)∂∂α(dadt)=𝒜2eprsinfcosδ(1+tan2α)-(𝒜1+3𝒜2)esinftanα-𝒜1esinftan3α-3𝒜2presinfcosδtan2α+𝒜3prcosδ(1+tan2α)cosα-𝒜3prcosδtanα(1+tan2α)sinα=0.
Finally, we can obtain
(55)=-𝒜1tan3-2𝒜2prcosδtan2α+𝒜3presinfcosδ(1+tan2α)×[cosα-tanαsinα]-(𝒜1+3𝒜2)tanα+𝒜2prcosδ=0.
This is a transcendental equation which is so difficult to be solved analytically, but in the following, we will consider some special cases.


## 7.2. Special Case I

Ignoring the absorbed and diffusely reflected parts of the radiation, of the coefficients (47) become
(56)𝒜1=[(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2]𝒜2=2ρsAi,𝒜3=0.
Thus, (55) becomes
(57)∂∂α(dadt)=tan3α+2𝒜2𝒜1prcosδtan2α+(1+3𝒜2𝒜1)tanα-𝒜2𝒜1prcosδ=0.
Let us now define the new coefficients:
(58)𝒦1=2𝒜2𝒜1prcosδ,𝒦2=(1+3𝒜2𝒜1),𝒦3=-𝒜2𝒜1prcosδ.
Equation (57) becomes clear as a cubic equation:
(59)tan3α+𝒦1tan2α+𝒦2tanα+𝒦3=0,
which has the same solution as given by (27), but upon replacing the coefficients  ℬ's with the newly defined coefficients 𝒦's we obtain
(60)(tanα1*)a=-13𝒦1+𝕂+𝒦12-3𝒦29𝕂,(tanα2,3*)a=-13𝒦1-12𝕂-𝒦12-3𝒦218𝕂±i32𝒦1[𝕂-𝒦12-3𝒦29𝕂],
where
(61)𝕂=(16𝒦1𝒦2-12𝒦3-127𝒦13+(127𝒦13𝒦3-1108𝒦12𝒦22-16𝒦1𝒦2𝒦3+127𝒦23+14𝒦22)1/2)1/3.
Since tanα2,3* are imaginary solution, therefore, these two roots are rejected.   And thus the only considered root is the first root  tanα1*.


## 7.3. Special Case II

Ignoring the absorbed and diffusely reflected parts of the radiation and assuming that the solar wind and the coronal mass ejection do not depend on the cone angle, then (57) can be reduced to the following quadratic equation:
(62)dadt=F·Ξa=P2pr2μ(1-e2)2×{[(Aim)CME2ρsAicosαcos2α+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)](es·Ξa)+Aicosα(2ρscosαsinαcosδ)(et·Ξa)+Aicosα(2ρscosαsinαsinδ)×(ew·Ξa)(Aim)CME},dadt=P2pr2μ(1-e2)2×{[(Aim)CME2ρsAicosαcos2α+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×(es·esinfe^s)+Aicosα(2ρscosαsinαcosδ)×(et·pre^t)(Aim)CME},
which can be simplified to
(63)dadt=P2pr2μ(1-e2)2×{[(Aim)CME2ρsAicosαcos2α+((Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2)]×esinf+Aiprcosα(2ρscosαsinαcosδ)}.
Consider the following definitions:
(64)𝒜0=(Aim)SWnSWMPvSW2+(Aim)CME(MCMEτCME)vCME2,𝒜2=2ρsAiP.
Equation (63) can be simplified to
(65)dadt=2pr2μ(1-e2)2×[pr𝒜0esinf+𝒜2(esinfcosα+prsinαcosδ)cos2α].
Differentiating  (65) with respect to α  and equating the result to zero, we obtain
(66)∂∂α(dadt)=-3esinfcos2αsinα+prcosαcosδ-3prsin2αcosαcosδ=0.
Dividing on 3cos2αsinα yields
(67)∂∂α(dadt)=-esinf+13prsinαcosαcosδ-3prtanαcosδ=0,
which can be written as
(68)(13cosαsinα){-3esinfcosαsinα+prcosδ-3prtanαcosαsinαcosδ}=0.
Equating the numerator to zero and rearranging the terms yield
(69)2(1+ecosfesinf)cosδtan2α+3tanα-1+ecosfesinfcosδ=0,
which has the solution of quadratic equations
(70)(tanα*)a=-3±9+8((1+ecosf)/esinf)2cos2δ4((1+ecosf)/esinf)cosδ.
This equation represents the control law that maximizes the rate of change of semimajor axis. It computes the actual sail pitch/cone angle profile necessary to increase the semimajor axis at a maximum rate (Figure 3).(a) The optimum cone angle that maximizes the rate of change of semimajor axis when taking the positive sign in the numerator. (b) The optimum cone angle that maximizes the rate of change of semimajor axis when taking the negative sign in the numerator.
(a)
(b)


## 7.4. Optimal Control Law for the Rate of Change of Eccentricity

When comparing (70) and (69), we observe that the required direction to maximize certain variable is given by tanα~=(1+ecosf)/esinf=(p/r)/esinf, which represents  the transverse component divided by the radial component of the primer vector. Thus, in case of eccentricity, set
(71)tanα~=ΞteΞse=cosf+(r/p)(e+cosf)sinf=(2+ecosf)cotf(1+ecosf)+ecscf(1+ecosf).
The maximum rate of change of the orbital eccentricity follows directly, taking into account the simplification addressed in special Case II (Figure 4), as
(72)(tanα*)e=(-3±9+8((2+ecosf)cotf+ecscf(1+ecosf))2cos2δ)×(4((2+ecosf)cotf+ecscf(1+ecosf))cosδ)-1.(a) The optimum cone angle that maximizes the rate of change of eccentricity when taking the positive sign in the numerator. (b) The optimum cone angle that maximizes the rate of change of eccentricity when taking the negative sign in the numerator.
(a)
(b)


## 7.5. Optimal Control Law for the Rate of Change of Inclination

Changes in the out-of-plane orbital elements, such as inclination and right ascension of the ascending node, can be effected by the use of simple switching functions.  For maximum rate of change of inclination, the solar sail thrust can be directed alternately above and below the orbit plane every half orbit by thesign function which has +1 or −1.In what follows, we will find an optimum cone angle that maximizes the rate of change of inclination:(73)∂∂α(didf)=r3μpsign[cos(f+ω)]∂∂αW=0,
from which one can obtain
(74)3ρsρd[1-2tan2α]=1cosα[tan2α-1].
Setting x=tan2α,  the above equation can be written as
(75)x3+(-1-ρd3ρs)x2+(2ρd-3ρs3ρs)x+(1-ρd3ρs)=0,
which has the solution
(76)(tan2α1*)i=19ρs(1+ρd)+(ρ~0+ρ~11/2)1/3+19ρs2[(1+ρd)2-9ρs(2ρd-3ρs)]+181ρs[16(-1-ρd)(2ρd-3ρs)-12(1-ρd)-127(-1-ρd)3](ρ~0+ρ~11/2)-1/3,(tan2α2,3*)i=19ρs(1+ρd)-12(ρ~0+ρ~11/2)1/3-19×18ρs2[(1+ρd)2-9ρs(2ρd-3ρs)]×(ρ~0+ρ~11/2)-1/3±i36ρs(1+ρd)×[(ρ~0+ρ~11/2)1/3-181ρs2×[(1+ρd)2-9ρs(2ρd-3ρs)]×(ρ~0+ρ~11/2)-1/3181ρs2],
where
(77)ρ~0=145ρs2(1+ρd)(2ρd-3ρs)-16ρs(1-ρd)+127×27ρs3(1+ρd)3,ρ~1=1729ρs4(1+ρd)3(9-9ρd)-1162ρs(1-ρd2)(2ρd-3ρs)+1729ρs3(2ρd-3ρs)3+136ρs2(2ρd-3ρs)2-18748(1+ρd)2(1-ρd)2.
As a special case when ignoring the diffusive part of radiation, we obtain
(78)3ρs[1-2tan2α]=0,⟹(tanα*)i=12⟹αi*=35.26°.
This result is in agreement with the previous results.  If we substitute the control law into the inclination equation in Lagrange planetary equations and integrate over one orbit, we can obtain the change in inclination per orbit.  We find that the change in inclination per orbit is independent of orbit radius and only depends on the sail lightness number  (straight number relation) as follows:
(79)Δi=88.2β(degreesperorbit).
However, closer orbits to the Sun have shorter orbit periods, and so the time to achieve an overall inclination change is shorter. This optimal angle is 35.26°, which enables us to maximize angular momentum.


## 7.6. Optimal Control Law for the Rate of Change of Ascending Node

The procedure is the same as in optimal control law for inclination(80)∂∂α(dΩdf)=r3μpsinisign[sin(f+ω)]∂∂αW=0⟹3ρsρd[1-2tan2α]=1cosα[tan2α-1],
which has the same solutions as those maximizing the inclination and is in good agreement with the previous results. If we substitute this control law into the ascending node in Lagrange planetary equations and integrate over one orbit, we can obtain the change in ascending node per orbit. We find that the change in ascending node per orbit is independent of orbit radius and only depends on the sail lightness number as follows:
(81)ΔΩ=(88.2βsini)(degreesperorbit).
However, closer orbits to the Sun have shorter orbit periods, and so the time to achieve an overall inclination change is shorter (Figure 5).Figure 5
The optimum change in the longitude in the ascending nodeΔΩ versus the inclination sin i and sail lightness number β.


## 7.7. Optimal Control Law for the Rate of Change of Argument of Perihelion

When comparing (70) and (69), we observe that the required direction to maximize a certain variable is given by tanα~=(1+ecosf)/esinf=(p/r)/esinf, which represents  the transverse component divided by the radial component of the primer vector. Thus, in case of argument of perihelion (Figure 6), we obtain
(82)tanα~=ΞtωΞsω=(1+(r/p))sinfcosf=(1+(rp))tanf=(2+ecosf1+ecosf)tanf,(tanα*)ω=-3±9+8(((2+ecosf)/(1+ecosf))tanf)2cos2δ4(((2+ecosf)/(1+ecosf))tanf)cosδ.(a) The optimum cone angle that maximizes the rate of change of argument of perihelion when taking the positive sign in the numerator. (b) The optimum cone angle that maximizes the rate of change of argument of perihelion when taking the negative sign in the numerator.
(a)
(b)


## 8. Conclusion

We can conclude the new contributions of the present work as follows. The force model on the sailcraft configurations is amended with the solar wind and coronal mass ejection forces. The solar radiation pressure model is modified so as to include the absorbed, the specularly reflected, and the diffusely reflected into the model. Some analytical control laws with some mentioned input constraints for optimizing sailcraft dynamics in Heliocentric orbit using Lagrange’s planetary equations are mathematically explored. Optimum force vector in a required direction is maximized by deriving optimal sail cone angle. Ignoring the absorbed and diffusely reflected parts of the radiation, some special cases are obtained. New control laws that maximize thrust to obtain certain required maximization in some particular orbital element are obtained.



---

*Source: 353056-2013-11-05.xml*