Analytical model for starshade formation flying with applications to exoplanet direct imaging observation scheduling

Gabriel J. Soto; Dmitry Savransky; Rhonda Morgan

doi:10.1117/1.JATIS.7.2.021209

11 January 2021 Analytical model for starshade formation flying with applications to exoplanet direct imaging observation scheduling

Gabriel J. Soto, Dmitry Savransky, Rhonda Morgan

Author Affiliations +

Journal of Astronomical Telescopes, Instruments, and Systems, Vol. 7, Issue 2, 021209 (January 2021). https://doi.org/10.1117/1.JATIS.7.2.021209

Abstract

We present an analytical model for the desired kinematics of the starshade-telescope relative motion during exoplanet direct imaging observations. We combine this model with an existing deadbanding strategy published by the NASA JPL S5 Team to define a dynamics framework for deadbanding simulations. Global results of these simulations show that the fuel usage and the number of observation interruptions vary as a function of the target star ecliptic coordinates and time, meaning there exist optimal times to observe particular targets. We combine these results with the telescope pointing constraints due to the relative position of the Sun and other bright solar system objects. We show that optimally scheduling an observation could result in up to 30 more min of integration time and 26 fewer interruptions per observation, improvements of almost 300% in some cases. We also show how phasing the start time of the telescope on its halo orbit is paramount for ensuring optimal observations, providing up to 68 additional min and 31 fewer interruptions per observation. Choosing an optimal halo phasing can also increase, for some near-ecliptic target stars, the fraction of a year that the target is observable from a few percent to more than 30%.

1. Introduction

Creating optimal observation schedules for exoplanet direct imaging with a starshade requires careful consideration of mission constraints. Of the two primary starshade flight modes—formation flying and slewing—the former requires higher precision control. To image exo-Earths at contrast levels on the order of $10^{- 10}$ ,¹^,² the starshade must fly in very tight formation with the telescope-star line of sight (LOS) at a distance of about a dozen Earth radii. The centroid of the starshade cannot stray more than 1 m from the LOS to the star. This constraint prevents diffraction effects from degrading the contrast levels and ruining the sensing of the Arago spot.³^–⁵ This level of precision control, fortunately, has been demonstrated both analytically and in experimental simulations.⁶^–⁸ The starshade, in many designs, carries out this precision control using a bipropellant engine: it is allowed to drift from its desired positioning and fires its thruster only when it reaches the lateral tolerance of 1 m from the LOS.⁴^,⁹ However, every burst of propellant used to maneuver the starshade back toward its desired trajectory creates a transient plume that reflects light and prevents a successful observation. The bi-prop engine creates discrete plumes that disperse quickly rather than the continuous cloud formed by solar electric propulsion. Many starshade designs instead use a separate SEP engine to conduct slew maneuvers.⁴^,¹⁰^,¹¹ Propellant cloud interruptions from discrete, bi-prop maneuvers add overhead to the observation time. Optimal trajectories around the desired starshade position should maximize the drift time between pulses. It is also important to know these overheads for different target stars as the telescope and starshade orbit the Sun in complicated patterns. This will help us to predict optimal times where overheads are low for a particular target star of interest. We, therefore, categorize the aims of this study into three major areas:

1. present a thorough, analytical model for the telescope-star LOS and the nominal starshade kinematics;
2. develop a dynamics model based on deadbanding strategies outlined by the JPL S5 Team in Ref. 9 to simulate observations; and
3. perform observation simulations with a large target list to determine when and where optimal observation conditions exist.

First, we develop a model to describe the complex geometry of the starshade relative to the telescope. Previous works, including Sirbu et al.,¹² have developed equations for the relative position based on vector components. We have instead elected to use an analytical model using Euler angles to define the LOS from the telescope to a target star solely as a function of time and the ecliptic coordinates of the star. Doing so yields concise formulas for the required pointing of the telescope as well as the necessary starshade kinematics for perfect formation flying—the condition where the starshade sits perfectly on the LOS at its constant separation distance. We also use this model to define the injection velocities required to start an observation. These velocities serve as boundary conditions for solving slew trajectories.¹³ These analytical provide a clear look into the variables that affect the kinematics of the starshade.

Although these formulas define perfect formation flying for the starshade, the starshade is not required to maintain these conditions exactly. Gravitational and other disturbance forces will pull it away from the desired track toward some offset trajectory. Rather than continuously correct these offsets—which would require more fuel usage and ruin observations with continuous plumes—we can instead allow the starshade to drift from the defined desired track onto an offset trajectory and only perform corrective maneuvers when it reaches the deadbanding limits set by telescope contrast requirements. We study the dynamics of these offset trajectories using the framework developed by Flinois et al.⁹ from the JPL S5 Team¹⁴ which assumes that, for short trajectories, lateral differential forces are constant. Flinois et al.⁹ demonstrated Monte Carlo simulations of the offset dynamics using high-fidelity guidance, navigation, and control algorithms and integrators. We conduct similar simulations of the offset dynamics, first by defining the necessary frames using our analytical starshade model and then integrating our equations of relative motion within a Python environment using standard scipy and numpy packages.¹⁵ Our open-source software provides accessibility and the ability to simulate deadbanding of the starshade during observations with a wide set of parameters, including any target star coordinate at any point throughout the mission time and any telescope orbit.

Our model and simulations provide results that are particularly useful as heuristics for exoplanet direct imaging simulations. Full end-to-end direct imaging missions for starshades can be simulated using EXOSIMS.¹⁶^,¹⁷ At every decision step, a scheduling algorithm selects the next best star to observe. Many metrics are considered in this scheduling step: how likely a star is to have an orbiting exoplanet that is observable to your instrument,¹⁸^–²⁰ whether each star is observable due to the relative location of the Sun or other bright solar system objects,¹³^,²⁰ and the required fuel needed to slew to that next star,¹³^,²¹ among others. The station-keeping fuel usage has not been used as one of these metrics. The station-keeping fuel used during each observation was previously calculated with a simplistic model within EXOSIMS. Increasing the precision of the station-keeping fuel determination is important because long integration times are required to photometrically detect and spectrally characterize the atmospheres of smaller, dimmer exoplanets. The station-keeping fuel use is dominated by the lateral forces pushing the starshade off the desired LOS trajectory. The lateral component of the disturbance acceleration varies with the LOS configuration relative to the Sun, Earth, and Moon. This causes variation in the fuel use and required number of thruster firings as a function of time (or its position in orbit) and the location of the target star on the sky. We determine optimal times for observing particular stars. These are combined with keepout constraints to ensure that whenever stars are observable, the starshade can optimally observe the target. The phasing of the telescope on its Sun–Earth L2 halo orbit—where the telescope is positioned at the mission start time—is also key to ensuring optimally timed exoplanet observations.

In the work presented, we first define the reference frames and general dynamics used to analyze motion in Sec. 2. We then give an overview of the starshade formation flying dynamics, defining all forces acting on the starshade in Sec. 3. We then develop expressions for the desired starshade formation flying kinematics and parallax correction acceleration in Sec. 4. Keepout angles are defined for the telescope in Sec. 5. In Sec. 6, we integrate the relative dynamics of the starshade and connect them to the model given by Flinois et al.⁹ We list the algorithms used to simulate deadbanding within the lateral tolerances. Finally, we present global trends of key metrics in Sec. 7, such as daily fuel usage and thruster firings per observation, as a function of the target star chosen and the mission elapsed time. We also show the relationship between the keepout zones and halo orbit phasing, presenting an optimal phasing for given target stars that will always lead to optimal observations.

2. Dynamical Background

The formation flying dynamics of the starshade are derived in Sec. 3 using the Sun–Earth–Moon model of the circular restricted three body problem (CR3BP).²² In this model, the Earth and Moon are treated as a single entity located at their center of mass. The Earth–Moon barycenter and the Sun are then assumed to orbit their common Sun–Earth–Moon barycenter (SEMB) on circular orbits. This model captures most of the spacecraft dynamics and enables computational efficiency when conducting simulations to produce formation-keeping costs.

2.1.

Definition of Reference Frames

We first define an inertial frame $I$ with basis vectors ${\hat{e}}_{1}$ , ${\hat{e}}_{2}$ , ${\hat{e}}_{3}$ and origin $O$ coinciding with the SEMB point. Position components in the $I$ frame are represented along each axis as $(X, Y, Z)$ using the convention in Ref. 22. Scalar time derivatives of components will be denoted with single dots $(\dot{X}, \dot{Y}, \dot{Z})$ and double dots $(\ddot{X}, \ddot{Y}, \ddot{Z})$ for velocities and accelerations, respectively. In this frame, the Sun and the combined Earth–Moon system both move in circular orbits about $O$ . We represent the Earth–Moon system barycenter as the point $O^{'}$ shown in Fig. 1. We then define the rotating frame $R$ , which is central to the CR3BP model, relative to our inertial frame. The $R$ rotates with the Sun and $O^{'}$ at a constant rate and shares the same origin $O$ with basis vectors $\hat{x}$ , $\hat{y}$ , $\hat{z}$ . Position components in the $R$ frame are represented along each axis as $(x, y, z)$ . Scalar time derivatives of these components will be denoted the same way as for the inertial components. The $I$ and $R$ frames are initially aligned at some reference epoch which we assume to be March 20, 2034. On this date, the vernal equinox coincides with the alignment of the lunar orbital nodes with the ecliptic and is a useful reference for ephemerides approximations.

Fig. 1

Diagram of telescope-starshade-target star vector configuration in the inertial frame relative to origin $O$ . Orbital path of the reference halo orbit is shown both in the inertial frame and as a closed orbit relative to the Sun–Earth rotating frame (which is rotated from the inertial frame by an angle $t$ ). Star coordinates are defined by ecliptic longitude ( $λ$ ) and latitude ( $β$ ).

We use the CR3BP canonical units for all values unless otherwise specified. The distance unit is equal to the Sun- $O^{'}$ separation, or 1 AU, and the time unit is scaled such that the mean motion of the Sun and $O^{'}$ orbits is 1. As a consequence, the period of these orbits is equal to $2 π$ and the gravitational parameter $G$ is 1. The mass unit is the sum of the primary masses and we define the dimensionless parameter $μ^{*}$ , given as

Eq. (1)

μ^{*} \equiv \frac{m_{\oplus} + m_{☾}}{m_{⊙} + m_{\oplus} + m_{☾}} .

The individual masses for the primary bodies are then defined as $μ_{⊙} \equiv 1 - μ^{*}$ and $μ_{\oplus + ☾} \equiv μ^{*}$ . The individual Earth and Moon mass parameters $μ_{\oplus}$ and $μ_{☾}$ are their masses scaled by the sum of the primary masses. With this parameter, we define the positions of the Sun and $O^{'}$ , respectively, as $r_{⊙ / O} = - μ^{*} \hat{x}$ and $r_{O^{'} / O} = (1 - μ^{*}) \hat{x}$ in components of the rotating frame. In inertial frame components, the positions vary as a function of time as

Eq. (2)

{[r_{⊙ / O}]}_{I} = - μ^{*} [\begin{matrix} \cos t \\ \sin t \\ 0 \end{matrix}],

and

Eq. (3)

{[r_{O^{'} / O}]}_{I} = (1 - μ^{*}) [\begin{matrix} \cos t \\ \sin t \\ 0 \end{matrix}],

where

t

is the time elapsed from the reference epoch.

We also assume that any star in our catalog or target list will be stationary relative to this inertial frame $I$ for the duration of the mission. The position of a given star $i$ is defined relative to the origin $O$ as

Eq. (4)

{[r_{i / O}]}_{I} = \frac{1}{ϖ} [\begin{matrix} \cos β \cos λ \\ \cos β \sin λ \\ \sin β \end{matrix}]

in inertial frame components. The ecliptic latitutde

β

is measured from the ecliptic plane and the ecliptic longitude

λ

is measured from the vernal equinox at our reference epoch. The parallax angle

ϖ

is defined as

\tan ϖ = \frac{1}{S_{i}}

where

S_{i}

is the stellar distance from

O

in canonical units. Parallax angles and other stellar parameters are catalogued and calculated using the astropy Python package.²³ We can approximate

\tan ϖ \approx ϖ

. Here we take

O

as an approximation of the solar system barycenter. In reality, the solar system barycenter and the defined position of

O

has varied on average by about 0.005 AU over the last 20 years according to ephemerides from the JPL Horizons system.²⁴ However, the worst possible viewing scenario—on a star at 1 pc located at an ecliptic pole—would produce an

\sim 3.7 mas

directional error. A star at 30 pc would create an error of 0.17 mas. When compared to a telescope field of view of

12 \times 12 arc \sec^{2}

for HabEx at visible wavelengths, we can assume that treating stellar spherical coordinates relative to

O

has a negligible effect for our simulation purposes.¹¹ Proper motion can be applied to the model for high-fidelity star tracking, but the effects will be negligible for the mission times considered. Other effects like stellar aberration should be considered in the future for a higher fidelity model.

2.2.

General Dynamic Model

The motion of any spacecraft in either the inertial or rotating frame will be driven by their respective equations of motion. In the inertial frame, the inertial acceleration of some particle $P$ is governed, through Newton’s Second Law of Motion by

Eq. (5)

{{}^{I}a}_{P / O} = \sum f_{P} = g_{P / ⊙} + g_{P / (\oplus + ☾)},

where the superscript

I

represents a vector derivative with respect to the inertial frame,

\sum f_{P}

is the summation of specific forces on the particle

P

,

g_{P / ⊙}

is the gravitational specific force on

P

from the Sun, and

g_{P / (\oplus + ☾)}

is the combined gravitational specific force on

P

from the Earth and Moon. These forces are

Eq. (6)

g_{P / ⊙} = - μ_{⊙} \frac{r_{P / ⊙}}{{‖ r_{P / ⊙} ‖}^{3}},

and

Eq. (7)

g_{P / (\oplus + ☾)} = - μ_{\oplus + ☾} \frac{r_{P / O^{'}}}{{‖ r_{P / O^{'}} ‖}^{3}},

noting that

r_{P / ⊙} = r_{P / O} - r_{⊙ / O}

and

r_{P / O^{'}} = r_{P / O} - r_{O^{'} / O}

. We define the inertial frame components of

P

as

{[r_{P / O}]}_{I} \equiv {[X_{P}, Y_{P}, Z_{P}]}^{T}

to generate the relative positions to the Sun and

O^{'}

using Eqs. (2) and (3). The resultant acceleration of

P

from Eq. (5) is, therefore, dependent only on the position of

P

and time. We refer to these equations of motion, namely Eq. (5), as the inertial CR3BP (or ICR3BP) since they share the same nomenclature and canonical units.

The equations of motion in the rotating frame, the CR3BP for the same particle $P$ , are given similarly as

Eq. (8)

{}^{R}{a_{P / O}} = ν_{P / O} + \sum f_{P},

where the superscript

R

represents a vector derivative with respect to the rotating frame

R

and

ν_{P / O}

encompasses the noninertial accelerations, which in rotating frame components is

Eq. (9)

{[ν_{P / O}]}_{R} = [\begin{matrix} x_{P} + 2 {\dot{y}}_{P} \\ y_{P} - 2 {\dot{x}}_{P} \\ 0 \end{matrix}] .

The specific gravity forces are the same as in Eq. (5) but are written in rotating frame components. The $R$ frame components of both sides of Eq. (8) are

Eq. (10)

[\begin{matrix} {\ddot{x}}_{P} \\ {\ddot{y}}_{P} \\ {\ddot{z}}_{P} \end{matrix}] = [\begin{matrix} x_{P} + 2 {\dot{y}}_{P} + \frac{(- μ - x_{P}) (1 - μ)}{{‖ r_{P / ⊙} ‖}^{3}} + \frac{(1 - μ - x_{P}) μ}{{‖ r_{P / O^{'}} ‖}^{3}} \\ y_{P} - 2 {\dot{x}}_{P} - \frac{(1 - μ) y_{P}}{{‖ r_{P / ⊙} ‖}^{3}} - \frac{μ y_{P}}{{‖ r_{P / O^{'}} ‖}^{3}} \\ - \frac{(1 - μ) z_{P}}{{‖ r_{P / ⊙} ‖}^{3}} - \frac{μ z_{P}}{{‖ r_{P / O^{'}} ‖}^{3}} \end{matrix}] .

Since the Sun and $O^{'}$ are stationary in the $R$ frame, Eq. (8) is time-independent and therefore an autonomous system.

2.3.

Halo Orbit of the Telescope

The space telescope is assumed to be on a halo orbit about the Sun–(Earth + Moon) $L_{2}$ .²⁵^,²⁶ The distance of $L_{2}$ from the Earth–Moon barycenter is $‖ r_{L_{2} / O^{'}} ‖ = 0.010075$ AU. Halo orbits are generated using the rotating frame equations of motion in Eq. (10). Initial conditions that lead to a periodic halo orbit are found using an iterative differential correction process combined with a single shooting algorithm.²¹^,²⁷ The reference halo generated for our simulations is shown in Fig. 2. It has a period of about 179 days, with amplitudes of 800,000 km in the $y$ direction, 400,000 km in the $x$ direction and also above the ecliptic plane. We calculate sufficient points throughout its period to create an interpolant which allows us to effectively know the full telescope trajectory—position, velocity, and acceleration—as a function of time. In rotating frame components and derivatives, these are

Eq. (11)

{[\begin{matrix} r_{T / O} (t) \\ {{}^{R}v}_{T / O} (t) \\ {{}^{R}a}_{T / O} (t) \end{matrix}]}_{R} = {[x_{T} (t), y_{T} (t), z_{T} (t), {\dot{x}}_{T} (t), {\dot{y}}_{T} (t), {\dot{z}}_{T} (t), {\ddot{x}}_{T} (t), {\ddot{y}}_{T} (t), {\ddot{z}}_{T} (t)]}^{T},

where the subscript

T

denotes the telescope. The first six terms, the telescope position and velocity components, are computed by integrating Eq. (10) in first-order form; the last three acceleration components are the resultant accelerations from directly evaluating Eq. (10) with the previous six terms. We do not consider station-keeping of the telescope on this halo orbit, assuming that the telescope is responsible for maintaining that reference halo trajectory.

Fig. 2

Reference halo orbit for the space telescope relative to $L_{2}$ . Telescope positions at different halo orbit phasings are shown in red.

Although generating the reference halo orbit in Fig. 2, we assumed the starting point at our reference epoch to be at the southern-most point of the halo. Since the CR3BP equations are autonomous, we can choose a different injection point on the halo orbit for the telescope, independent of time. We implement this variation of injection points as a halo orbit phase time $Δ t_{P}$ . Different halo orbit phasings are shown in Fig. 2 over half the halo period. Although this may affect the initial transfer fuel costs of the telescope and starshade from Earth, there exist many trajectory options with marginal $Δ v$ changes.²⁸ The effects of varying this phase for different missions are discussed in Sec. 7.4.

Though it would be convenient to frame our model in the CR3BP rotating frame and take advantage of its time independence during numerical integration, the starshade station-keeping problem is inherently time dependent—we must either keep track of the moving Sun, Earth, and Moon in the inertial frame or the moving stars in the rotating frame. We choose to develop our station-keeping models in the inertial frame to ease the additions of perturbation forces and because the line of sight definitions are more intuitive in inertial frame components. The telescope’s trajectory on the halo orbit, in the inertial frame as a function of the computed rotational frame components, is

Eq. (12)

{[\begin{matrix} r_{T / O} \\ {}^{I}{v_{T / O}} \\ {{}^{I}a}_{T / O} \end{matrix}]}_{I} = [\begin{matrix} X_{T} (t) \\ Y_{T} (t) \\ Z_{T} (t) \\ {\dot{X}}_{T} (t) \\ {\dot{Y}}_{T} (t) \\ {\dot{Z}}_{T} (t) \\ {\ddot{X}}_{T} (t) \\ {\ddot{Y}}_{T} (t) \\ {\ddot{Z}}_{T} (t) \end{matrix}] = [\begin{matrix} x_{T} (t) \cos t - y_{T} (t) \sin t \\ x_{T} (t) \sin t + y_{T} (t) \cos t \\ z_{T} (t) \\ {\dot{x}}_{T} (t) \cos t - {\dot{y}}_{T} (t) \sin t - Y_{T} (t) \\ {\dot{x}}_{T} (t) \sin t + {\dot{y}}_{T} (t) \cos t + X_{T} (t) \\ {\dot{z}}_{T} (t) \\ {\ddot{x}}_{T} (t) \cos t - {\ddot{y}}_{T} (t) \sin t - 2 {\dot{Y}}_{T} (t) + X_{T} (t) \\ {\ddot{x}}_{T} (t) \sin t + {\ddot{y}}_{T} (t) \cos t + 2 {\dot{X}}_{T} (t) + Y_{T} (t) \\ {\ddot{z}}_{T} (t) \end{matrix}] .

These time histories for the inertial halo trajectories are continuous and can also be interpolated over a finite mission time.

3. Starshade Formation Flying Dynamics

Science requirements for a starshade observation dictate that the starshade be positioned at a constant separation distance from the telescope during an observation. This ensures the best possible contrast at the desired observing wavelength. We, therefore, define a desired position $D$ for the starshade and state that $r_{D / T}$ be constant along the LOS to a star. In reality, the distance from the LOS can vary from this desired location. The actual position of the starshade $S$ is $r_{S / O}$ . Both vectors are shown in Fig. 1. The purpose of formation-keeping maneuvers is then to keep the difference between the desired and actual positions, namely

Eq. (13)

r_{S / D} = r_{S / O} - r_{D / O},

within the lateral limits that still guarantee high-contrast imaging. These are 1 m laterally from the LOS and 250 km axially along the LOS. Differentiating the positional difference vector in Eq. (13) twice with respect to time results in

Eq. (14)

{{}^{I}a}_{S / D} = {{}^{I}a}_{S / O} - {{}^{I}a}_{D / O} .

Expanding the acceleration of the desired starshade position and grouping terms results in

Eq. (15)

{{}^{I}a}_{S / D} = ({{}^{I}a}_{S / O} - {{}^{I}a}_{T / O}) - {{}^{I}a}_{D / T} .

The dynamics of the formation flying error—the deviation of the actual starshade from its desired trajectory due to the disturbance acceleration ${{}^{I}a}_{S / D}$ —consists of the relative dynamics $({{}^{I}a}_{S / O} - {{}^{I}a}_{T / O})$ between the starshade and telescope modified by the acceleration ${{}^{I}a}_{D / T}$ of the desired relative position relative to the telescope. We use Eq. (15) to numerically integrate the deviation of the starshade trajectory and implement formation keeping costs in Sec. 6.2. We derive and discuss the three individual terms of Eq. (15) before discussing simulations. The inertial acceleration of the telescope ${{}^{I}a}_{T / O}$ is known as a function of time from Eq. (12). The acceleration of the actual starshade ${{}^{I}a}_{S / O}$ is discussed in more detail in Sec. 3.1. The desired starshade acceleration ${{}^{I}a}_{D / T}$ , often called the parallax correction, is derived from geometric constraints on the starshade kinematics in Sec. 4 and is discussed specifically in Sec. 4.5.

3.1.

Summed Forces on the Starshade

We develop equations for the starshade directly in the inertial frame using an enhanced version of the ICR3BP in Eq. (5) to find ${{}^{I}a}_{S / O}$ . To add higher fidelity to the ICR3BP, we split Eq. (7) into separate Earth and lunar gravity specific forces. A solar radiation pressure (SRP) specific force is also added as a perturbing force. The full equations of motion of the starshade is therefore

Eq. (16)

{{}^{I}a}_{S / O} = \sum f_{S} = g_{S / ⊙} + g_{S / \oplus} + g_{S / ☾} + p_{S / ⊙},

where

p_{S / ⊙}

is the specific force on

S

due to SRP. The gravitational specific force from the Sun on the starshade

S

g_{S / ⊙}

is the same as in Eq. (2). The Earth and lunar gravity specific forces are derived from a more detailed model in Sec. 3.2. The new set of equations in Eq. (16) will be referred to as the enhanced ICR3BP.

3.2.

Earth and Lunar Gravity

We add more detail to the Earth–Moon model by including their individual orbital positions about $O^{'}$ . We assume a circular orbit for the Moon about $O^{'}$ with radius $r$ and inclination $i$ . The precession of the lunar nodes relative to the ecliptic plane is represented by a constant lunar nodal rate ${\dot{Ω}}$ of 1 revolution per 18.59 years, converted to canonical units. We assume that the lunar nodes cross the ecliptic at our reference epoch. The lunar ephemerides are then

Eq. (17)

{[r_{☾ / O^{'}}]}_{I} = - r_{☾} [\begin{matrix} \sin (n_{☾} t) \sin ({\dot{Ω}}_{☾} t) \cos i_{☾} + \cos (n_{☾} t) \cos ({\dot{Ω}}_{☾} t) \\ \sin (n_{☾} t) \cos ({\dot{Ω}}_{☾} t) \cos i_{☾} - \cos (n_{☾} t) \sin ({\dot{Ω}}_{☾} t) \\ \sin (n_{☾} t) \sin i_{☾} \end{matrix}] .

where

r_{☾} = 384,748 km

,

n_{☾} = \frac{2 π}{T_{☾}}

,

i = 5.15 \deg

, and

T_{☾}

is the lunar orbital period of 29.53 days converted to canonical units. The full lunar orbit is calculated as

r_{☾ / O} = r_{☾ / O^{'}} + r_{O^{'} / O}

. The lunar gravity specific force on the starshade is

Eq. (18)

g_{S / ☾} = - μ_{☾} \frac{r_{S / ☾}}{{‖ r_{S / ☾} ‖}^{3}}

noting that the gravitational parameter

G

is equal to 1 because of canonical units. We can similarly define an orbit for the Earth, though more simply, as

Eq. (19)

{[r_{\oplus / O^{'}}]}_{I} = - r_{\oplus} [\begin{matrix} \cos (n_{\oplus} t) \\ \sin (n_{\oplus} t) \\ 0 \end{matrix}],

where

n_{\oplus} = n_{☾}

and

r_{\oplus} = \frac{μ_{☾}}{μ_{\oplus}} r_{☾} = 4730 km

. The full lunar orbit is calculated as

r_{\oplus / O} = r_{\oplus / O^{'}} + r_{O^{'} / O}

. The gravitational specific force due to the Earth on the starshade is then

Eq. (20)

g_{S / \oplus} = - μ_{\oplus} \frac{r_{S / \oplus}}{{‖ r_{S / \oplus} ‖}^{3}} .

Both forces are time dependent and have varying magnitude ranges as shown in Table 1.

Table 1

Maximum specific force magnitudes on the desired starshade position. The maximum lateral and axial components are separately calculated, directions are relative to the LOS.

	Total acceleration	Lateral component	Axial component
$‖ g_{D / ⊙} ‖$ ( $mm / s^{2}$ )	5.85	5.84	5.85
$‖ g_{D / \oplus} ‖$ ( $μ m / s^{2}$ )	318.49	282.43	318.48
$‖ g_{D / ☾} ‖$ ( $μ m / s^{2}$ )	8.15	6.91	8.15
$‖ p_{D / ⊙} ‖$ ( $μ m / s^{2}$ )	3.31	0.02	3.31
$‖ {{}^{I}a}_{T / O} ‖$ ( $mm / s^{2}$ )	6.10	6.10	6.10

3.3.

Solar Radiation Pressure Force

Due to the size of the starshade, which can have a radius as large as 72 m for the HabEx design, we implement an SRP force on the starshade. The model of the SRP force is given by Refs. 29 and 30. The SRP specific force is

Eq. (21)

p_{S / ⊙} = 2 \frac{P A}{m_{S}} \cos α [b_{1} {\hat{r}}_{S / ⊙} + (b_{2} \cos α + b_{3}) \hat{n}],

where

m_{S}

is the mass of the starshade,

α

is the pitch angle of the starshade pointing relative to the Sun-starshade line,

P

is the solar radiation pressure at

L_{2}

,

A

is the starshade area based on the starshade radius

R_{S}

, and

b_{1}

,

b_{2}

, and

b_{3}

are the optical parameters whose starshade-relevant values are taken from Ref. 31. The solar radiation pressure

P

is an inverse square law

Eq. (22)

P = P_{0} {(\frac{r_{0}}{r_{S / ⊙}})}^{2},

where

P_{0} = 4.563 μ N / m^{2}

and

r_{0} = 1

AU. The pitch angle can be computed from the relation

Eq. (23)

\cos α = {\hat{r}}_{S / ⊙} \cdot \hat{n} .

Assuming no tilting of the starshade during an observation, the normal vector

n

is assumed to point along the LOS to a star as is explained in Sec. 4.

4. Desired Kinematics of the Starshade

The last term to consider from Eq. (15) is the parallax correction term. These are geometrical constraints placed on the starshade desired motion to ensure its trajectory keeps up with the changing LOS from the telescope to a star. We first define the geometry of the LOS in general terms and use that definition to derive the desired starshade kinematics in the inertial frame.

We now have analytic formulas for the LOS from the telescope to any target star assuming a known orbit for the telescope in an inertial frame. Knowing the look vector ${\hat{b}}_{3}$ , we can now define the full kinematics—position, velocity, and acceleration—of the starshade on its nominal or fully aligned trajectory.

4.1.

Tracking the Line of Sight Using Euler Angles

The LOS from the telescope $T$ to a target star $i$ is defined using the vector triad shown in Fig. 1, also written as

Eq. (24)

r_{i / O} = r_{i / T} + r_{T / O} .

The inertial position of the star $r_{i / O}$ is given in Eq. (4) and the telescope position $r_{T / O}$ is known from Eq. (12). We now solve for the relative position of the star to the telescope $r_{i / T}$ and do so component-wise as has been done previously in Ref. 32. We have elected, however, to further introduce Euler angles to define the target star LOS from the telescope. We develop analytical formulas for these LOS tracking angles and their derivatives. This helps us to define tracking rates for the telescope to aid mission planning and design. These Euler angles also help us rigorously to define the lateral and axial directions in which starshade perturbations will occur and provide formulation that leads to greater insight of the starshade formation flying problem.

We use a 3-2-3 Euler angle sequence to define the LOS from the telescope to the target star as illustrated in Fig. 3. We first rotate the inertial frame $I$ about its third axis by an azimuth angle $θ$ , effectively a rotation on the ecliptic plane to create an auxiliary frame $A$ . We then rotate this frame about its second axis by a polar angle $ϕ$ , creating a new frame $B$ . We complete the Euler angle set in Sec. 6.1 once we define the forces on the starshade. These first two Euler angles are sufficient to define the direction from the telescope to the target star relative to the inertial frame. The relative position of the target star to the telescope in inertial frame components is fully defined as

Eq. (25)

{[r_{i / T}]}_{I} = S_{T} [\begin{matrix} \sin ϕ \cos θ \\ \sin ϕ \sin θ \\ \cos ϕ \end{matrix}],

where

S_{T}

is the distance from the telescope to the target star. Note that all variables, components, and time derivatives of components are given in canonical units unless otherwise specified. We combine Eqs. (24) and (25) to solve for the Euler angles and the telescope-star distance:

Eq. (26)

S_{T} = \frac{\sin β - ϖ Z_{T}}{ϖ \cos ϕ},

Eq. (27)

\tan θ = \frac{ϖ Y_{T} - \sin λ \cos β}{ϖ X_{T} - \cos β \cos λ},

Eq. (28)

\tan ϕ = \frac{\sqrt{{(ϖ X_{T} - \cos β \cos λ)}^{2} + {(ϖ Y_{T} - \sin λ \cos β)}^{2}}}{\sin β - ϖ Z_{T}},

and all are given as functions of known or measurable values: time, the ecliptic coordinates and parallax distance of the target star, and the position of the telescope in inertial frame components. Note that we have dropped the

(t)

notation from the telescope coordinates only for clarity. An unfortunate consequence is the existence of a singularity in the angular kinematics at

ϕ = 0

and

ϕ = π

. However, current mission concepts only consider stars up to

\sim 30 pc

away as candidates for exoplanet detections, which means we do not encounter targets located exactly at the ecliptic poles. Moreover, our selected Sun keepout angles, discussed in Sec. 5, exclude ecliptic coordinates near the poles. If, in the future, target stars near the ecliptic are selected, one could redefine the Euler angle set to place the singularities elsewhere (perhaps changing the definition of

ϕ

to be a latitude rather than colatitude).

Fig. 3

Auxiliary frames used to define the target star LOS from the telescope. Two Euler angles are used: an azimuth angle ( $θ$ ) and a polar angle ( $ϕ$ ).

4.2.

Line of Sight Tracking Rates

The tracking rates for the LOS to a given target star are a function of star coordinates and time. Differentiating the Euler angles in time yields the angular rate of change of the azimuthal angle

Eq. (29)

\dot{θ} = \cos^{2} θ \frac{d}{d t} (\frac{ϖ Y_{T} - \sin λ \cos β}{ϖ X_{T} - \cos β \cos λ}),

Eq. (30)

\dot{θ} = ϖ (- {\dot{X}}_{T} \sin θ + {\dot{Y}}_{T} \cos θ),

and the angular rate of change of the polar angle

Eq. (31)

\dot{ϕ} = \cos^{2} ϕ \frac{d}{d t} [\frac{\sqrt{{(ϖ X_{T} - \cos β \cos λ)}^{2} + {(ϖ Y_{T} - \sin λ \cos β)}^{2}}}{\sin β - ϖ Z_{T}}],

Eq. (32)

\dot{ϕ} = ϖ [\frac{\cos ϕ ({\dot{X}}_{T} \cos θ + {\dot{Y}}_{T} \sin θ) + {\dot{Z}}_{T}}{\sin ϕ}],

both as functions of time dependent variables—the trajectory of the telescope and Euler angles—as well as the time-independent parallax angle of the target star. The tracking rates scale with parallax, meaning they are smaller for target stars that are further away. We differentiate again to find the angular accelerations of the tracking angles:

Eq. (33)

\ddot{θ} = ϖ [({\ddot{Y}}_{T} - {\dot{X}}_{T} \dot{θ}) \cos θ - ({\ddot{X}}_{T} + {\dot{Y}}_{T} \dot{θ}) \sin θ],

Eq. (34)

\ddot{θ} = \frac{ϖ}{2} [ϖ {\dot{X}}_{T}^{2} \sin (2 θ) - 2 ϖ {\dot{X}}_{T} {\dot{Y}}_{T} \cos (2 θ) - ϖ {\dot{Y}}_{T}^{2} \sin (2 θ) - 2 {\ddot{X}}_{T} \sin θ + 2 {\ddot{Y}}_{T} \cos θ],

and

Eq. (35)

\ddot{ϕ} = \frac{ϖ}{\sin ϕ} {[({\ddot{X}}_{T} \cos θ + {\ddot{Y}}_{T} \sin θ - {\dot{X}}_{T} \dot{θ} \sin θ + {\dot{Y}}_{T} \dot{θ} \cos θ) \cos ϕ - ({\dot{X}}_{T} \cos θ + {\dot{Y}}_{T} \sin θ) \dot{ϕ} \sin ϕ + {\ddot{Z}}_{T}] - {\dot{ϕ}}^{2} \cos ϕ} .

We have also given approximate values for the second derivatives after eliminating all $ϖ^{2}$ terms since these terms are very small and approach machine precision. The angular accelerations are then

Eq. (36)

\ddot{θ} \approx ϖ ({\ddot{Y}}_{T} \cos θ - {\ddot{X}}_{T} \sin θ),

Eq. (37)

\ddot{ϕ} \approx ϖ [\frac{\cos ϕ ({\ddot{X}}_{T} \cos θ + {\ddot{Y}}_{T} \sin θ) + {\ddot{Z}}_{T}}{\sin ϕ}] .

An important note is that, while we are assuming a halo orbit for the telescope, these equations work for any telescope orbit so long as its time history is known.

4.3.

Positioning the Starshade

The desired trajectory of the starshade $D$ relative to the telescope $T$ is derived using our geometric LOS definitions as a function of the known telescope inertial trajectory, star location, and time. The desired placement of the starshade is perfectly along the LOS. Adherence to the axial distance constraints along the LOS enables observations in specific science bands based on the telescope and starshade size.⁷ In reality, the actual starshade position $S$ can vary by up to 250 km axially but only 1 m laterally from the LOS.⁴ We assume, for the desired trajectory $D$ , a constant separation distance $s$ and measure deviations from the desired trajectory using Eq. (15). The starshade, therefore, can be located relative to the telescope using the $B$ -frame which, by definition, points toward the target star as follows:

Eq. (38)

r_{D / T} = s {\hat{b}}_{3} .

The position vector of the starshade from the origin of the inertial frame, in inertial frame components, is therefore

Eq. (39)

{[r_{D / O}]}_{I} = {[r_{D / T} + r_{T / O}]}_{I} = [\begin{matrix} s \sin ϕ \cos θ + X_{T} \\ s \sin ϕ \sin θ + Y_{T} \\ s \cos ϕ + Z_{T} \end{matrix}] .

Knowing the position vector of the starshade as a function of time and star coordinates allows us to define initial conditions when integrating the equations of motion in Eq. (15).

4.4.

Starshade Velocities in Inertial and Rotating Frames

Before fully deriving the desired starshade velocity, we must first define the angular velocities of each auxiliary frame in Fig. 3. The angular velocity of the $B$ frame in the inertial frame is given by

Eq. (40)

{{}^{I}ω}^{B} = {{}^{I}ω}^{A} + {{}^{A}ω}^{B} = \dot{θ} {\hat{a}}_{3} + \dot{ϕ} {\hat{a}}_{2},

or, in component form:

Eq. (41)

{[{{}^{I}ω}^{B}]}_{I} = [\begin{matrix} - \dot{ϕ} \sin θ \\ \dot{ϕ} \cos θ \\ \dot{θ} \end{matrix}] .

Due to our assumption of constant separation distance, the starshade velocity relative to the telescope in the $B$ -frame is zero, or ${{}^{B}v}_{D / T} = 0$ . We use this to derive the inertial velocity relative to the telescope as

Eq. (42)

{{}^{I}v}_{D / T} = {{}^{B}v}_{D / T} + {{}^{I}ω}^{B} \times r,

or

Eq. (43)

{[{{}^{I}v}_{D / T}]}_{I} = s [\begin{matrix} \dot{ϕ} \cos ϕ \cos θ - \dot{θ} \sin ϕ \sin θ \\ \dot{ϕ} \cos ϕ \sin θ + \dot{θ} \sin ϕ \cos θ \\ - \dot{ϕ} \sin ϕ \end{matrix}] .

The magnitude of this relative inertial velocity is shown in Fig. 4 over half the reference halo orbit period in SI units. Finally, we complete the vector triad with the telescope’s inertial velocity to get

Eq. (44)

{{}^{I}v}_{D / O} = {{}^{I}v}_{D / T} + {{}^{I}v}_{T / O},

where the inertial velocity of the telescope

{{}^{I}v}_{T / O}

is given in Eq. (12). This final inertial velocity, given in inertial components, is the desired starshade velocity for formation flying during an observation.

Fig. 4

Magnitude in SI units of the (a)–(c) inertial velocity and (d)–(f) acceleration of the starshade relative to the telescope on a reference halo orbit as a function of time after mission start. Target stars are placed at distance of 1 pc from origin $O$ . Both colorscales are logarithmic. Relative velocities are small and the relative accelerations can be considered negligible for stars further than 1 pc.

Throughout this study, we will primarily use kinematic derivatives in the inertial frame; however, kinematics relative to the rotating frame $R$ of the Sun and Earth are still of interest for solving other starshade dynamics problems, primarily starshade slews between observations. Note that because of our defined canonical units the angular velocity of the rotating frame relative to the inertial frame is ${{}^{I}ω}^{R} = - {{}^{R}ω}^{I} = {[0,0, 1]}_{I}^{T}$ . We can derive the rotating frame velocity necessary to begin an observation on any star (i.e., the injection velocity at the boundaries of the slew trajectory) by first computing the velocity of the nominal starshade relative to the telescope

Eq. (45)

{{}^{R}v}_{D / T} = {{}^{I}v}_{D / T} + {{}^{R}ω}^{I} \times r_{D / T},

or

Eq. (46)

{[{{}^{R}v}_{D / T}]}_{I} = s [\begin{matrix} \dot{ϕ} \cos ϕ \cos θ - \dot{θ} \sin ϕ \sin θ + \sin ϕ \sin θ \\ \dot{ϕ} \cos ϕ \sin θ + \dot{θ} \sin ϕ \cos θ - \sin ϕ \cos θ \\ - \dot{ϕ} \sin ϕ \end{matrix}],

given in inertial frame components. We can convert to

R

-frame components by rotating through an angle

t

about the

\hat{z}

in Fig. 1. The full rotating frame velocity of the starshade relative to the origin is

Eq. (47)

{{}^{R}v}_{D / O} = {{}^{R}v}_{D / T} + {{}^{R}v}_{T / O},

where

{{}^{R}v}_{T / O}

is defined in Eq. (11). These calculations ensure that the proper injection velocity is used in the slew boundary value problem rather than assuming the starshade matches inertial velocities with the telescope as in Refs. 13 and 21. In Sec. 6.2, we will define the optimal injection velocity assuming we allow offset motion from the desired trajectory.

4.5.

Parallax Correction Acceleration for Starshade

Finally, we can derive the full, desired starshade acceleration in the inertial frame as well as the parallax correction term for Eq. (15). First, we define the angular acceleration vector

Eq. (48)

{{}^{I}ω}^{B} \equiv^{I} \frac{d}{d t} ({{}^{I}ω}^{B}) =^{I} \frac{d}{d t} ({{}^{I}ω}^{A} + {{}^{A}ω}^{B}) = \ddot{θ} {\hat{a}}_{3} + (\ddot{ϕ} {\hat{a}}_{2} - \dot{ϕ} \dot{θ} {\hat{a}}_{1}),

or

Eq. (49)

{[{}^{I}{\dot{ω}}^{B}]}_{I} = [\begin{matrix} - \ddot{ϕ} \sin θ - \dot{ϕ} \dot{θ} \cos θ \\ \ddot{ϕ} \cos θ - \dot{ϕ} \dot{θ} \sin θ \\ \ddot{θ} \end{matrix}],

which is equivalent to the time derivative of Eq. (40). Noting, again, that the desired starshade position relative to the telescope is constant in the

B

-frame, we find that

Eq. (50)

{{}^{I}a}_{D / T} =^{I} \frac{d}{d t} ({{}^{I}v}_{D / T}) =^{B} \frac{d}{d t} ({{}^{I}v}_{D / T}) + {{}^{I}ω}^{B} \times ({{}^{I}v}_{D / T}) .

We make the following substitutions:

Eq. (51)

{{}^{I}a}_{D / T} =^{B} \frac{d}{d t} ({{}^{B}v}_{D / T} + {{}^{I}ω}^{B} \times r_{D / T}) + {{}^{I}ω}^{B} \times ({{}^{B}v}_{D / T} + {{}^{I}ω}^{B} \times r_{D / T}),

Eq. (52)

{{}^{I}a}_{D / T} = {}^{I}{\dot{ω}}^{B} \times r_{D / T} + {{}^{I}ω}^{B} \times {{}^{I}ω}^{B} \times r_{D / T} .

We then calculate the desired starshade acceleration in the inertial frame relative to the telescope, or the parallax correction acceleration:

Eq. (53)

{[{{}^{I}a}_{D / T}]}_{I} \approx s [\begin{matrix} \ddot{ϕ} \cos ϕ \cos θ - \ddot{θ} \sin ϕ \sin θ \\ \ddot{ϕ} \cos ϕ \sin θ + \ddot{θ} \sin ϕ \cos θ \\ - \ddot{ϕ} \sin ϕ \end{matrix}]

after eliminating terms on the order of

{\dot{θ}}^{2}

and

{\dot{ϕ}}^{2}

, since they inherently contain

ϖ^{2}

terms and would, therefore, be negligible for stars at distances

> 1 pc

. To get the desired starshade acceleration relative to the origin, we need to add the inertial telescope acceleration:

Eq. (54)

{{}^{I}a}_{D / O} = {{}^{I}a}_{D / T} + {{}^{I}a}_{T / O},

which is given in Eq. (12). We have found, however, that for the specific range of star parallaxes in our target list (namely,

> 1 pc

) the relative acceleration of the starshade to the telescope in Fig. 4 is orders of magnitude smaller than the telescope’s acceleration, which is on the order of

6 mm / s^{2}

. These formation flying accelerations would be more notable for closer objects, such as bodies in our solar system.

5. Keepout Angles

With our new Euler angle definitions, we define analytic formulas for the telescope keepout angles. Keepout angles limit the telescope pointing at a particular time in its mission: any star that is within a given angular separation from the Sun or another bright object is not observable, as light from the foreground object would contaminate the observation, and requires the telescope to wait until it is within viewing range. Upper bounds on these pointing angles are also imposed due to incident light requirements on the solar panels as well as the danger of contamination by reflected light from the starshade during observations.⁵ Keeping track of the dynamic keepout regions is vital for observation scheduling and is especially important when considering overheads and shared telescope time with other science instruments.

To define the keepout regions for some bright body in the solar system $B$ , we assume full knowledge of its ephemerides in the inertial frame as functions of time: ${[r_{B / O}]}_{I} \equiv {[x_{B}, y_{B}, z_{B}]}_{I}^{T}$ . We define the angle between the target star and bright body as

Eq. (55)

\cos κ_{B} = \frac{r_{B / T} \cdot {\hat{b}}_{3}}{‖ r_{B / T} ‖},

where

r_{B / T} = r_{B / O} - r_{T / O}

, so that

Eq. (56)

\cos κ_{B} = \frac{(x_{B} - x_{T}) \sin ϕ \cos θ + (y_{B} - y_{T}) \sin ϕ \sin θ + (z_{B} - z_{T}) \cos ϕ}{\sqrt{{(x_{B} - x_{T})}^{2} + {(y_{B} - y_{T})}^{2} + {(z_{B} - z_{T})}^{2}}} .

This expression is a function only of the ecliptic coordinates of the target star and time (if the orbit of both the telescope and bright body are known). We reiterate the fact that we can substitute any orbit for the telescope in this formulation. Further simplifications are made in accordance with the enhanced ICR3BP framework by assuming simple orbits for the Sun, Earth, and Moon. The ephemerides of each are given in Eqs. (2), (19), and (17) and make Eq. (56) more explicitly a function of time. Similar simplifications can be made for the other planets of the solar system, though they are not considered here.

Our implementation of keepout constraints is demonstrated in Fig. 5. We analytically compute the angles from the telescope LOS to the Sun, Earth, and Moon as a function of time and star coordinates. We then check to see if each of these angles fit within the upper and lower bounds $κ_{\min}$ and $κ_{\max}$ found in Table 2. Stars that meet the criteria

Eq. (57)

κ_{⊙, \min} < κ_{⊙} (t) < κ_{⊙, \max},

Eq. (58)

κ_{\oplus, \min} < κ_{\oplus} (t) < κ_{\oplus, \max},

Eq. (59)

κ_{☾, \min} < κ_{☾} (t) < κ_{☾, \max},

at a specific mission time

t

are shown in the unshaded regions of Fig. 5. Stars that do not meet the criteria are within the yellow, blue, and gray regions of the Sun, Earth, and Moon keepouts, respectively. We implement

κ_{\max}

for the Sun only to account for reflected light from the starshade. This model can be updated to include a

κ_{\max}

term for Earthshine. As a conservative case, the minimum keepout constraints for all bright bodies is set by the telescope scarf angle, assumed to be 45 deg. However, it is uncertain based on the current literature how the phase function of the Earth and Moon would diminish light contamination on the telescope. For these reasons, we include two separate classifications of minimum keepout angle

κ_{\min}

for the Earth and Moon: case 1 (C1) with an optimistic minimum keepout angle of 5 deg and case 2 (C2) with more conservative minimum keepout constraints set by the telescope scarf angle. These time-dependent keepout constraints will affect the timings of observations and are used in Sec. 7 to filter results.

Fig. 5

Keepout regions in ecliptic coordinates 30 days after mission start. Only targets in the unshaded regions are observable at the time shown. Yellow, blue, and gray regions depict keepout zones due to the Sun, Earth, and Moon, respectively. Paths throughout the sky relative to the telescope are shown for the three bodies.

Table 2

Solar system keepout constraints and parameters in canonical units.

Body	κmin,C1 (deg)	κmin,C2 (deg)	κmax (deg)	Mass	Orbital radius (rB)	ω
Sun	45	45	83	0.999997	$3.040433 \times 10^{- 6}$	1
Earth	5	45	180	$3.040433 \times 10^{- 6}$	0.999997	1
Moon	5	45	180	$3.6923866 \times 10^{- 8}$	$2.571881 \times 10^{- 3}$	12.37

6. Simulating Deadbanding Maneuvers

We now have a full model for the starshade dynamics relative to the desired trajectory. Deviations from the desired trajectory are numerically integrated to study when and how often correction maneuvers must be used to maintain lateral and axial tolerances during an observation. From the right-hand side of Eq. (15), we define a lateral disturbance acceleration:

Eq. (60)

Δ a_{S} \equiv (\sum f_{S} - {{}^{I}a}_{T / O}) - {{}^{I}a}_{D / T},

where the parallax correction term can be neglected. The final equations of motion are then

Eq. (61)

{{}^{I}a}_{S / D} = Δ a_{S}

showing that the disturbance acceleration on the starshade

S

drives its deviation from its desired trajectory

{{}^{I}a}_{S / D}

.

6.1.

Parabolic Trajectory Approximations

We can make a further simplification of the relative dynamics in Eq. (61) by assuming that the specific forces at the actual starshade location $S$ are similar to those at the desired location $D$ or

Eq. (62)

\sum f_{S} \approx \sum f_{D} .

Because these specific forces evolve over long timescales, noting that the period of the halo orbit is 6 months, we can assume that the forces are constant over the relatively shorter observation times and rewrite Eq. (61) as

Eq. (63)

{{}^{I}a}_{S / D} \approx Δ a_{D},

where

Eq. (64)

Δ a_{D} \equiv \sum f_{D} - {{}^{I}a}_{T / O}

and is constant throughout the simulated observation time. Of primary interest are the lateral disturbances relative to the LOS which drive motion in the

{\hat{b}}_{1} - {\hat{b}}_{2}

plane. We define a new frame

C

to better study the motion lateral to the LOS as shown in Fig. 6. The frame is centered at the desired position of the starshade

D

and rotated by a final roll angle

ψ

to complete the Euler angle set relative to the inertial frame

I

. We define the axial and lateral components of the lateral disturbance acceleration as

Eq. (65)

a_{A} \equiv Δ a_{D} \cdot {\hat{b}}_{3},

Eq. (66)

a_{L} \equiv Δ a_{D} - a_{A} {\hat{b}}_{3},

where the axial direction points along the telescope-target LOS and the lateral direction lies on a plane perpendicular to the LOS. The roll angle is defined as

Eq. (67)

\tan ψ = \frac{a_{L} \cdot {\hat{b}}_{1}}{- a_{L} \cdot {\hat{b}}_{2}},

such that

{\hat{c}}_{2} = - {\hat{a}}_{L}

. In this frame, we can define the parabolic trajectories derived in Ref. 9. We need these approximations as an initial guess for the required velocities to maintain the deadbanding strategy. Throughout deadbanding, we still integrate the full relative equations of motion in Eq. (61) and use the parabolic guesses from Eq. (63) as initial conditions.

Fig. 6

Final frame for dynamics, rotated from B-frame so that the lateral differential force points opposite the second axis.

We define the coordinates of the new $C$ -frame as ${[r_{S / D}]}_{C} = {[η, σ, ζ]}_{C}^{T}$ to derive the parabolic trajectories. We scale these coordinates by the deadbanding radius $R_{d}$ according to Ref. 9. Similarly, the velocities in this frame are scaled by $\sqrt{R_{d} a_{L}}$ , time by $\sqrt{R_{d} / a_{L}}$ , and accelerations by $a_{L}$ where $a_{L} \equiv ‖ a_{L} ‖$ is the magnitude of the lateral disturbance acceleration.⁹ The second derivatives of these new coordinates are given by

Eq. (68)

{[{{}^{I}a}_{S / D}]}_{C} \equiv {[\begin{matrix} \ddot{η} \\ \ddot{σ} \\ \ddot{ζ} \end{matrix}]}_{C} = {[\begin{matrix} 0 \\ - 1 \\ a_{A} / a_{L} \end{matrix}]}_{C} .

Noting that the terms on the right-hand side of this equation are all constants, Flinois et al.⁹ developed equations for parabolic flight within the deadbanding radius. The idea behind this is that the global optimal trajectory would be one that travels back and forth along the diameter of the deadbanding radius. This is in the $C$ -frame where the lateral disturbance acceleration always points downward. Figure 7 demonstrates example trajectories within the deadbanding tolerances in the plane lateral to the LOS. Although one can maximize the time of flight for an individual trajectory, the ending location of the starshade centroid within the deadbanding radius is not guaranteed to be in an optimal location for a long subsequent flight. The strategy, then, is to move the starshade as quickly as possible toward the “well” of the deadbanding radius, point $[0, - 1]$ in the $η - σ$ plane in normalized units, so that it can travel as close as possible along the vertical diameter (i.e., the longest possible flight). The necessary equations for parabolic flight that target an endpoint at the deadbanding “well” require an intercept point tangent to the deadbanding radius as shown in Fig. 7. This intercept point, as a function of the initial location on the deadbanding radius ${[η_{0}, σ_{0}]}^{T}$ , is

Eq. (69)

{[η_{i}, σ_{i}]}^{T} = {\begin{cases} [η_{0}, σ_{0}] & σ_{0} > 0.5 \\ \frac{1}{\sqrt{2}} [sgn (η_{0}) \sqrt{1 + σ_{0}}, \sqrt{1 - σ_{0}}] & σ_{0} \leq 0.5 \end{cases},

where the sgn function returns 1 or

- 1

for positive and negative arguments, respectively. The velocity components in this frame are

{[{{}^{I}v}_{S / D}]}_{C} \equiv {[\dot{η}, \dot{σ}, \dot{ζ}]}_{C}

; the initial velocity for a parabola with the initial location

{[η_{0}, σ_{0}]}^{T}

, intercept point at

{[η_{i}, σ_{i}]}^{T}

and endpoint

{[0, - 1]}^{T}

are

Eq. (70)

{\dot{η}}_{0} = - \frac{sgn (η_{0})}{\sqrt{2}} \sqrt{σ_{i} (1 - σ_{i})},

Eq. (71)

{\dot{σ}}_{0} = - \frac{η_{i}}{σ_{i}} {\dot{η}}_{0} + \frac{η_{i} - η_{0}}{{\dot{η}}_{0}},

with time of flight

Δ t = \frac{η_{0}}{{\dot{η}}_{0}}

in scaled units. If the initial point is exactly at the bottom of the deadbanding radius, the velocities then become

Eq. (72)

{\dot{σ}}_{0} = {\begin{cases} (3 σ_{i} - 1) \sqrt{\frac{1 + σ_{i}}{2 σ_{i}}} & σ_{0} < 0.5 \\ (1 - σ_{i}) \sqrt{\frac{1 + σ_{i}}{2 σ_{i}}} & σ_{0} \geq 0.5 \end{cases},

with

{\dot{η}}_{0} = 0

. These are the optimal injection velocities for maximizing drift time.

Fig. 7

Example of a deadbanding simulation. Starshade starts at the nominal position $D$ with an exaggerated initial velocity. A thruster firing is triggered when crossing the inner threshold. The second drift ends at the “well" of the deadbanding radius while the third drift is time optimal.

6.2.

Algorithm for Starshade Deadbanding Trajectories

To integrate the equations of motion for the relative starshade dynamics, we need initial conditions. Given the choice of starting locations, it is best to start the starshade on the maximal drift trajectory. At the beginning of a deadbanding trajectory, regardless of time and target star, the starshade must then be positioned at the well of the deadbanding radius. We can combine the initial parabolic velocities in Eq. (72) with the desired starshade velocity in Eq. (44) as follows:

Eq. (73)

{{}^{I}v}_{S / O} = {{}^{I}v}_{S / D} + {{}^{I}v}_{D / O},

or

Eq. (74)

{[{{}^{I}v}_{S / O}]}_{I} = {{}^{I}C}^{C} [\begin{matrix} 0 \\ {\dot{σ}}_{0} \sqrt{R_{d} a_{L}} \\ 0 \end{matrix}] + [\begin{matrix} {\dot{X}}_{T} + \dot{ϕ} \cos ϕ \cos θ - \dot{θ} \sin ϕ \sin θ \\ {\dot{Y}}_{T} + \dot{ϕ} \cos ϕ \sin θ + \dot{θ} \sin ϕ \cos θ \\ {\dot{Z}}_{T} - \dot{ϕ} \sin ϕ \end{matrix}],

where

{{}^{I}C}^{C}

is the direction cosine matrix associated with rotation from the

C

to

I

frames. This initial state (if we position the starshade at the deadbanding “well”) leads to maximum drift times in the idealized model. We use these initial conditions to begin the deadbanding simulations. They can also be used as slewing boundary conditions as discussed in Sec. 4.4.

We first integrate the equations of motion stated in Eq. (61) using these idealized initial conditions and estimate a time of flight based on the parabolic trajectory. Note that within these equations, unlike the ideal model in Sec. 6.1, we track the changing forces due to the offset trajectory. Another thing to note is that the relative dynamics are framed in the inertial frame, so the initial conditions have to be rotated from the $C$ frame. We incorporate similar lateral thresholds to those described in Ref. 9 although with different values: an inner threshold $R_{d, inner}$ of 0.9 m and an outer threshold $R_{d, outer}$ of 0.95 m. These are implemented to allow small overshoots of the trajectories due to the fact that the initial velocity conditions are calculated using the ideal equations. Event functions within the initial value problem solver are triggered when the lateral position of the starshade crosses these thresholds. Inner threshold crossings are recorded while outer threshold crossings terminate the trajectory integration. If the outer threshold is crossed in the $+ σ$ direction of the $C$ -frame, an immediate burn is required; otherwise, the next trajectory is resolved at the exiting of the inner threshold.

We define a burn as fully impulsive—a discrete and instantaneous change in velocity at the crossing of a threshold. We note the position ${[η_{1}, σ_{1}]}^{T}$ of the current crossing point and use the ideal equations, again, to define an intercept point and the initial velocities ${[{\dot{η}}_{1}, {\dot{σ}}_{1}]}^{T}$ required for a parabolic flight ending at the well of the deadbanding radius. To calculate a new parabola, the lateral acceleration $a_{L}$ is recalculated at the current time. The $Δ v$ in the lateral direction is then the difference between the new, parabolic velocities and the integrated velocities from the last drift. To counteract axial drift along the LOS, we assume a brake-damping strategy, in which we cancel the $ζ$ -direction velocity at every lateral burn. This strategy is not meant to cross the entire, much larger axial control region ( $\pm 250 km$ ) but instead to dampen the axial drift and lessen the impact on observations.⁹ The total $Δ v$ for a single crossing of the lateral radius is

Eq. (75)

Δ v_{n + 1} = \sqrt{{(- {\dot{ζ}}_{n, c})}^{2} + {({\dot{η}}_{n + 1} - {\dot{η}}_{n, c})}^{2} + {({\dot{σ}}_{n + 1} - {\dot{σ}}_{n, c})}^{2}},

where

{[{\dot{η}}_{n, c}, {\dot{σ}}_{n, c}, {\dot{ζ}}_{n, c}]}_{C}^{T}

are the integrated velocities of the

n

’th trajectory in the

C

-frame at the threshold crossing and

{[{\dot{η}}_{n + 1}, {\dot{σ}}_{n + 1}, {\dot{ζ}}_{n + 1}]}_{C}^{T}

are the parabolic guesses for the next trajectory toward the “well.” We catalogue these burns, drift times, and the number of crossings over a set observation period during an individual deadbanding simulation.

7. Global Trends and Simulation Results

Deadbanding simulations are now conducted for a target list of stars in the manner described in Sec. 6.2. The full results are intended to capture formation flying costs as a function of the ecliptic coordinates of the target, time, and the halo orbit phasing parameter. Metrics are first developed to describe the different costs associated with starshade formation flying. Resulting heatmaps of these metrics are then shown as a function of the parameters previously described. We then discuss how the lateral component of the disturbance acceleration on the starshade drives these cost metrics. Separate deadbanding simulations are also run without any axial control to demonstrate how shorter observations could be conducted while still meeting axial distance tolerances. Finally, we discuss the effect of varying the halo orbit phasing parameter on the formation flying costs.

7.1.

Station-Keeping Metrics for Observation Scheduling

We conducted station-keeping simulations for a grid of stars distributed uniformly in ecliptic longitude and latitude. Ecliptic longitudes were chosen in steps of 10 deg; ecliptic latitudes were chosen in steps of 10 deg from $- 80 \deg$ to 80 deg. Each simulation was allotted a set observation duration or simulated observation time $t_{obs}$ , during which we calculated deadbanding trajectories using the algorithm in Sec. 6.2. Integration times for observations can range from $< 1 h$ to multiple weeks.¹⁷ A 6-h $t_{obs}$ was, therefore, selected to represent short integration times. To capture the longer integration times, we conducted deadbanding simulations for every star at different days throughout a 1-year mission. We catalogued a number of different metrics from these simulations and present the results of these simulations in Fig. 8. Other starshade parameters are presented in Table 3.

Fig. 8

Maps of different station-keeping metrics, averaged over a simulation time of 6 h at three different mission times. The averaged metrics include total $Δ v$ , axial drift distance, number of thruster firings, drift times, fuel mass used per day, and the portion of the simulation time spent firing thrusters.

Table 3

Starshade deadbanding simulation parameters.

Parameter	Value
$R_{d}$	1 m
$R_{d, inner}$	0.9 m
$R_{d, outer}$	0.95 m
$m_{0}$	10,930 kg
$R_{S}$	36 m
$s$	76,600 km
$I_{sp}$	308 s
$F_{T}$	$2 \times 22 N$
$t_{obs}$	6 h
$S_{I}$	1 pc

All metrics are presented as functions of the target star ecliptic coordinates (represented as individual heat maps) at specific times in the mission. First, we catalogued the number of threshold crossings $N$ which corresponded to the required number of thruster firings throughout $t_{obs}$ . We then define the operator ${⟨ ⟩}_{obs}$ to symbolize the average value over a simulated observation. For a given metric $M$ and individual $M_{i}$ calculated during each trajectory $j$ between burns, the average of the metric is calculated using the operator as

Eq. (76)

{⟨ M ⟩}_{obs} = \frac{1}{N} \sum_{j}^{N} M_{j} .

The first two metrics we present in Fig. 8 are ${⟨ Δ v ⟩}_{obs}$ in the top row followed by the $ζ_{\max}$ in the second row. These are the average $Δ v$ used during the burn (both in the axial and lateral directions) and the final axial distance from $D$ along the LOS, respectively. Next, we catalogued the number of threshold crossings $N$ followed by the average drift time between thruster firings ${⟨ Δ t_{D} ⟩}_{obs}$ . The above metrics are independent of the type of chemical propulsion and thrusters assumed for station-keeping, they only vary with mass and area (due to SRP) of the starshade. For the last two metrics, we used parameters shown in Table 3 based on the HabEx starshade design.¹⁰^,¹¹ We define the fuel mass used during a single trajectory $j$ as

Eq. (77)

Δ m_{j} = m_{0} (1 - e^{- \frac{Δ v_{j}}{g_{0} I_{sp}}}),

from Tsiolkovsky’s ideal rocket equation where

I_{sp}

is the specific impulse of the thruster,

m_{0}

is the initial spacecraft wet mass, and

g_{0}

is a constant for the acceleration due to gravity.³³ We create a metric for the fuel usage per observation time, given as

Eq. (78)

{⟨ \dot{m} ⟩}_{obs} = \frac{1}{t_{obs}} \sum_{j}^{N} Δ m_{j},

where

N

is the number of thruster firings and

t_{obs}

is the observation time. We can also estimate the thruster firing time during each lateral correction using the equation for specific impulse

Eq. (79)

F_{T} = g_{0} I_{sp} \dot{m},

where

\dot{m}

is the mass flow rate and

F_{T}

is the maximum thrust of the engine in a single direction. We assume the thruster value in Table 3 which is a sixth of the maximum thrust listed in Ref. 11. If we assume the mass flow rate to be a constant for these short bursts, we can rewrite the equation as

Eq. (80)

Δ t_{T, j} = \frac{g_{0} I_{sp} Δ m_{j}}{F_{T}}

to estimate the duration of each thruster firing

j

. We then define a new metric

Eq. (81)

{⟨ f_{p} ⟩}_{obs} = \frac{1}{t_{obs}} \sum_{j}^{N} Δ t_{T, j},

which represents the fraction of time

f_{p}

(relative to the simulation time), in which the thruster is firing and therefore interrupting an observation.

The resulting formation flying costs are shown in Fig. 8, with colorbars chosen so that lighter colors represent lower costs (or more favorable conditions). We compared our results to those in Ref. 9, namely the worst case drift times at the worst possible gravity gradient. They estimated a worst case drift time for HabEx, with our same starshade separation distance, of $\sim 12 \min$ using a different telescope quasi-halo orbits. Our worst case drift ${⟨ t_{D} ⟩}_{obs}$ was $\sim 10 \min$ . Our results, therefore, match those in Ref. 9 with a 16% error despite the different halo orbit sizes and perturbation models. We also used a less strict $R_{d}$ value of 0.9 m rather than 0.7 m. Our largest lateral disturbance acceleration is on the order of $38 μ m / s^{2}$ compared to their $31.2 μ m / s^{2}$ due to our larger halo orbit.

With the lateral and axial control schemes described in Sec. 6.2, we see large variances of costs depending on the ecliptic coordinates of the intended target as well as the mission time. Mass used per day could vary up to 18 kg per day if a star is observed at the wrong time of year. The drift time per burn ${⟨ t_{D} ⟩}_{obs}$ averages between 10 min and 2.5 h. Although the fraction of observation time spent firing a thruster (and interrupting the observation) varies from a fraction of a percent to about 1.5%, there can be more than 30 interruptions per 6-h observation time. Observation scheduling, therefore, must take these time-varying costs to avoid conducting an observation of a particular star when the starshade is in an unfavorable, or high cost, configuration.

7.2.

Contribution Due to Lateral Disturbance Acceleration

The maximum axial distances in the second row of Fig. 8 show that our axial control law mitigates axial drifts very effectively. Recall that the axial positioning tolerance for the starshade is 250 km in either direction. The resulting 1 km variance in axial distance indicates that perhaps a less conservative control law should be employed in the future to relax fuel costs. We conducted separate deadbanding simulations without axial control—eliminating the axial $Δ v$ burn that was previously conducted at every lateral crossing of the starshade—to isolate the lateral control strategy and its effectiveness. This helps us simulate the costs of shorter observations, as these could successfully perform formation flying without a single axial maneuver and still remain within the axial tolerances. The results are shown in Fig. 9 using the same metrics as in Fig. 8. There are notable differences from the previous simulations in fuel costs. The ${⟨ Δ v ⟩}_{obs}$ values in the first row are now only in the lateral direction and are reduced by almost an order of magnitude. The fuel mass used per day and ${⟨ f_{P} ⟩}_{obs}$ are both also reduced when only conducting lateral burns. The maximum axial drift during the simulated 6-h observation, at most, is 15 km. With a simple extrapolation, we can estimate that the starshade will reach the axial tolerance of 250 km in about 4 days during the worst case conditions.

Fig. 9

Maps of different station-keeping metrics, averaged over a simulation time of 6 h at three different mission times. No axial burns are conducted here, so the starshade is free to drift in the axial direction.

The metrics in both Figs. 8 and 9 all exhibit similar patterns or regions of extrema. The ${⟨ Δ v ⟩}_{obs}$ plot was seemingly inverted when switching to a lateral-only control law. The resultant heatmap now matches the same pattern of extrema seen in the $N$ , ${⟨ Δ t_{D} ⟩}_{obs}$ , ${⟨ \dot{m} ⟩}_{obs}$ , and ${⟨ f_{P} ⟩}_{obs}$ heatmaps. Switching to a lateral-only control law did not significantly change the pattern of extrema nor the magnitudes of the $N$ and ${⟨ Δ t_{D} ⟩}_{obs}$ metrics. These cursory observations suggest that the pattern of extrema is driven by a common cause. We determine that cause to be the relevant component of the disturbance acceleration. The lateral component of the disturbance acceleration $a_{L}$ , for instance, drives the $N$ and ${⟨ Δ t_{D} ⟩}_{obs}$ metrics. The axial component $a_{A}$ drives $ζ_{\max}$ . The fuel costs have both an axial and lateral component; each ensures adherence to axial and lateral tolerances independently.

There are three distinct regions of extrema in the metric heatmaps: two poles and a banded region arcing over the ecliptic coordinates. Target stars within these regions are in favorable configurations for observation generally because of fewer interruptions and lower fuel costs. This pole-banded stucture is explained by several factors depicted in Fig. 10. In the first column, we plot the lateral component of the disturbance acceleration $a_{L}$ on the starshade desired position as a function of the target star ecliptic coordinates and the time throughout the mission. The location of the favorable metric regions in Fig. 8 match the regions of smaller lateral accelerations in Fig. 10. Small lateral accelerations allow the starshade to drift longer between the lateral tolerances, decrease the number of threshold crossings and reduce the amount of fuel needed for lateral corrections. The timing of the maximum $a_{L}$ coincides with the telescope being at its southern-most point—as well as the closest point to the Sun, Earth, and Moon—of its halo orbit, where gravitational forces are strongest.

Fig. 10

Lateral component of the lateral disturbance acceleration on the desired starshade position as a function of stellar location and time in column (a). Columns (b) and (c) show lateral force with combined keepout regions for cases 1 and 2, respectively. Different mission times are shown in each row.

We also note that the pole-banded structure of the metric maps and lateral differential accelerations all evolve over time. The formation flying geometry combines with the direction of the gravitational forces to produce this effect. The location of the pole regions, for example, is very near the location of the Sun, Earth, and Moon in the sky relative to the starshade. The Sun, Earth, and Moon locations are indirectly shown in the second and third columns of Fig. 10. Here we reproduce the first column with an added filter based, respectively, on cases 1 and 2 keepout regions discussed in Sec. 5. When the starshade is aligned with stars at the minimum poles, the differential acceleration $Δ a_{D}$ mostly aligns in the axial direction along the LOS due to the direction of gravitational forces as shown in Fig. 11(a). Lateral disturbances are, therefore, reduced. Because axial disturbances are increased, the starshade drifts further in the axial direction without any added control as shown in the second row of Fig. 9. Though the two poles align in the direction of $\sum f_{D}$ , which is dominated by gravitational forces of the Sun and Earth, the lateral differential acceleration is still nonzero due to the direction of ${{}^{I}a}_{T / O}$ shown in Fig. 11(a). The opposite pole is created when the telescope views a star opposite, or anti-parallel to, the LOS in Fig. 11(a): $\sum f_{D}$ is still aligned with the LOS but in the opposite direction. The banded structure is caused by similar reasons: in those regions, the specific forces on the starshade align nearly antiparallel with the telescope acceleration vector in the lateral direction to the LOS. This results in a small lateral component of the disturbance acceleration as shown in Fig. 11(b). We also note that, from Fig. 10, the magnitude of the lateral disturbance acceleration is periodic and shares the periodicity of the halo orbit—about 180 days.

Fig. 11

Two formation flying scenarios that lead to minimum lateral disturbance acceleration, or $a_{L}$ . The apparent Sun, Earth, and Moon locations in the sky are labeled. (a) The disturbance acceleration aligns mostly in the axial direction of the LOS resulting in pole minima and (b) the two components of the disturbance acceleration mostly cancel out in the lateral direction resulting in banded minima. Diagrams are meant to show exaggerated configurations.

7.3.

Scheduling Observations to Optimize Drift Time

The geometry and time evolution of this configuration drives the lateral disturbance acceleration on the starshade and directly drives many metrics shown in Figs. 8 and 9. We can simulate and tabulate these metrics as functions of target star coordinates and mission time to predict the location of favorable observation in the future and plan missions more effectively. We show the effects of optimal scheduling in Table 4 for a few target stars. We ran deadbanding simulations for each star with the same parameters in Table 3 used to generate the results in Figs. 8 and 9. Simulations were conducted at different times throughout a year, taken in intervals of 10 days. Metrics associated with fuel cost— ${⟨ Δ v ⟩}_{obs}$ , $⟨ \dot{m} ⟩_{obs}$ , and ${⟨ f_{P} ⟩}_{obs}$ —are split into lateral and axial components $L$ and $A$ , denoted by subscripts. We filtered the results using telescope keepout constraints, applying both case 1 and case 2 keepout scenarios. Stars showing only one column had the same results regardless of keepout conditions. We then found the optimal time $t_{opt}$ for scheduling an observation by choosing the date with the highest average drift time ${⟨ t_{D} ⟩}_{obs}$ . We also found the least optimal date regarding ${⟨ t_{D} ⟩}_{obs}$ and recorded the difference between the two in Table 4. The percent change from using the optimal over worst date is also recorded.

Table 4

Station-keeping simulation case study results with optimal time for observation scheduling using both case 1 (left) and case 2 (right) keepout conditions. Stars with one column had the same results for both keepout conditions. We also tabulate the difference in using the optimal versus worst observation times for each metric. The operator ⟨⟩obs represents averages over a 6-h observation for the given mission times.

	Beta pic	GJ 832	51 Eridani		GJ 179		47 UMa	HD 219143
$λ$ (deg)	82.54	308.62	67.31		72.44		149.07	23.74
$β$ (deg)	$- 74.42$	$- 32.47$	$- 24.31$		$- 15.93$		31.06	54.55
$S_{I}$ (pc)	19.75	4.97	29.40		12.36		13.80	6.55
$t_{opt}$ (d)	330	200	300	180	300	180	40	140
$t_{worst}$ (d)	270	170	190		190		10	180	200
Optimal versus worst	Case 1 + 2	Case 1 + 2	Case 1	Case 2	Case 1	Case 2	Case 1 + 2	Case 1	Case 2
${⟨ Δ t_{D} ⟩}_{obs}$ (min)	$+ 5 (+ 25 %)$	$+ 17 (+ 151 %)$	$+ 5 (+ 41 %)$	$+ 4 (+ 29 %)$	$+ 6 (+ 53 %)$	$+ 4 (+ 35 %)$	$+ 25 (+ 225 %)$	$+ 30 (+ 292 %)$	$+ 28 (+ 229 %)$
$N$	$- 4 (- 20 %)$	$- 19 (- 59 %)$	$- 9 (- 30 %)$	$- 7 (- 23 %)$	$- 10 (- 33 %)$	$- 8 (- 27 %)$	$- 22 (- 69 %)$	$- 26 (- 74 %)$	$- 20 (- 69 %)$
$ζ_{\max}$ (km)	$+ 0.2 (+ 1073 %)$	$+ 0.3 (+ 2451 %)$	$+ 0.1 (+ 306 %)$	$+ 0.2 (+ 452 %)$	$+ 0.2 (+ 520 %)$	$+ 0.2 (+ 540 %)$	$+ 0.3 (+ 949 %)$	$+ 0.2 (+ 112 %)$	$+ 0.2 (+ 379 %)$
${⟨ Δ v_{L} ⟩}_{obs}$ (mm/s)	$- 3 (- 22 %)$	$- 14 (- 67 %)$	$- 6 (- 29 %)$	$- 5 (- 24 %)$	$- 7 (- 35 %)$	$- 5 (- 27 %)$	$- 16 (- 76 %)$	$- 17 (- 75 %)$	$- 13 (- 70 %)$
${⟨ Δ v_{A} ⟩}_{obs}$ (mm/s)	$+ 16 (+ 1056 %)$	$+ 25 (- 2321 %)$	$- 17 (- 510 %)$	$+ 15 (+ 452 %)$	$- 22 (- 706 %)$	$+ 17 (+ 541 %)$	$+ 24 (- 976 %)$	$+ 38 (- 309 %)$	$+ 31 (- 565 %)$
${⟨ {\dot{m}}_{L} ⟩}_{obs}$ (kg/d)	$- 1 (- 38 %)$	$- 8 (- 87 %)$	$- 4 (- 51 %)$	$- 4 (- 42 %)$	$- 5 (- 57 %)$	$- 4 (- 47 %)$	$- 9 (- 92 %)$	$- 11 (- 94 %)$	$- 7 (- 91 %)$
${⟨ {\dot{m}}_{A} ⟩}_{obs}$ (kg/d)	$+ 4 (+ 825 %)$	$+ 5 (- 1002 %)$	$- 6 (- 387 %)$	$+ 5 (+ 323 %)$	$- 7 (- 504 %)$	$+ 5 (+ 370 %)$	$+ 4 (- 374 %)$	$+ 10 (- 154 %)$	$+ 6 (- 244 %)$
${⟨ f_{p, L} ⟩}_{obs}$ (%)	$- 0.1 (- 38 %)$	$- 0.7 (- 87 %)$	$- 0.3 (- 51 %)$	$- 0.3 (- 42 %)$	$- 0.4 (- 57 %)$	$- 0.3 (- 47 %)$	$- 0.7 (- 92 %)$	$- 0.9 (- 94 %)$	$- 0.6 (- 91 %)$
${⟨ f_{p, A} ⟩}_{obs}$ (%)	$+ 0.3 (+ 825 %)$	$+ 0.4 (- 1002 %)$	$- 0.4 (- 387 %)$	$+ 0.4 (+ 323 %)$	$- 0.5 (- 504 %)$	$+ 0.4 (+ 370 %)$	$+ 0.3 (- 374 %)$	$+ 0.8 (- 154 %)$	$+ 0.4 (- 244 %)$

The improvements due to optimal scheduling vary per star coordinates, but overall each star experiences an increase in drift time and decrease in number of thruster firings. The northern, mid-latitude stars experience the greatest improvements with HD 219143 seeing an average 30 more min of drift time and 26 fewer interruptions per observation. Stars spaced closely in ecliptic longitude also experience optimal observation configurations on similar dates, as evidenced by 51 Eri, and GJ 179. Although the lateral components of the fuel metrics decreased with the optimal drift time date, the fuel usage in the axial direction increased. Favorable drift times are caused by small lateral accelerations which inversely lead to large axial accelerations. For short target observations—those on the order of $t_{obs} = 6 h$ —the starshade would not exceed the axial distance tolerance and we can essentially ignore the axial fuel costs to simulate a lateral-only control law. For longer target observations, those that may take days, we can use these axial costs as a upper bound and plan less conservative axial control laws in the future. The fuel savings associated with lateral corrections, however, could result in up to 11 kg fewer spent per day. These simulations, metrics, and optimizations can be conducted for any target list to create optimal schedules for direct imaging missions that save fuel and increase the chances of successful exoplanet characterizations. Though we define the optimal observation date for each star based on drift times, one can develop more intricate cost functions that incorporate all metrics.

7.4.

Phasing the Halo Orbit

We consider one more parameter for the formation flying geometry: the phasing of the halo orbit discussed in Sec. 2.3. This phasing parameter allows us to alter the telescope-starshade geometry at mission start and its evolution throughout the mission. As demonstrated in Sec. 7.1, the time evolution of the starshade-star configuration determines the locations of minimum lateral disturbance acceleration. The halo phasing, therefore, allows us to control which—and when—stars are in favorable observation configurations. This approach is especially powerful in dealing with keepout constraints highlighted in Sec. 5. The time-dependent keepout regions filter out target coordinates at different mission times based on the relative positioning of the Sun, Earth, and Moon. Changing the halo phasing alters these relative positions so that certain stars are obscured at different times. We calculate the movement of keepout regions over a 1-year span and plot the portion of time each target is visible to the telescope in Fig. 12 for fixed halo phases and the two keepout conditions for Earth and Moon constraints. Stars at the highest latitudes are constantly obscured by the maximum keepout angle of the Sun. Stars at middle latitudes are visible for longer because they fit inside the keepout constraint annulus shown in Fig. 5. Central latitude visibility varies due to the motion and relative configuration of the Sun, Earth, and Moon for the case 2 keepout scenario; observation percentage is nearly constant for the Earth and Moon in case 1. Case 2, therefore, prevents more stars from being observed throughout the mission. We can control the amount of visibility of central and middle latitude targets in case 2 by varying the halo orbit phasing. Given high-priority targets, perhaps with known exoplanets or high completeness values, results from Fig. 12 can be used to determine favorable mission phasings for optimal mission design.

Fig. 12

Portion of time (over a 1-year span) during which individual target stars are out of keepout regions and therefore visible to the telescope. Case 1 (a) and case 2 (b) keepout conditions are applied, with stricter Earth and Moon minimum angles in the latter. Plotted with different halo phasing angles in each row.

We also use the halo orbit phasing to optimize the lateral disturbance acceleration on the starshade and combine it with the keepout results. We calculated the maximum lateral disturbance acceleration felt by the starshade when observing a grid of target stars in ecliptic coordinates throughout a year of mission time, applying the keepout filters when necessary. Figure 13 shows these results as we vary the halo orbit phasing. We also show the starting locations of the telescope corresponding to each of the halo orbit phasings (shown in the $R$ -frame).

Fig. 13

(a) Halo orbit in the Sun–Earth rotating frame with four different phasing times for the telescope. (b) Maximum lateral force on starshade over a 1-year mission lifetime as a function of target star ecliptic coordinates. Includes keepout changes over the mission, both case 1 (a) and 2 (b) scenarios. Halo orbit phasing time varies per row.

For a fixed halo orbit phase, we see larger lateral disturbance accelerations for higher latitude stars over the mission. For mid-latitude stars, regions of low lateral acceleration occur at specific longitudes. The halo orbit phasing dictates the longitudinal location of these minima and could provide reduction in lateral forcing by a factor of up to 8.8. These minima coincide with regions of longer visibility according to Fig. 12, which reinforces the importance of aligning the halo orbit phasing with high-priority targets.

7.5.

Scheduling Observations with Halo Orbit Phasing

We conducted a set of simulations to illustrate the use of these halo orbit phasings. In the previous section, we showed how the halo orbit phasings affect lateral disturbance accelerations on the starshade during observations. Here we conducted station-keeping simulations for a select number of stars and tabulated the resultant costs in Table 5. For each star, we ran deadbanding simulations on a coarse grid of different halo orbit phasings $Δ t_{P}$ (from 0 to 180 days in steps of 10 days) and a 1-year mission life span (with a step size of 5 days) for each phasing. We also applied keepout constraints on the results based on the two keepout cases as we did in Sec. 7.1. We found the optimal combination of halo orbit phasing $Δ t_{P}$ and mission time $t$ leading to maximum average drift times ${⟨ Δ t_{D} ⟩}_{obs}$ . We also found the worst case combination of $Δ t_{P}$ and $t$ for drift times. Table 5 shows the difference between metrics at those time combinations as well as the percent change from choosing the optimal over the worst case scenario. The results are more emphatic than those in Table 4: average drift times can be increased by more than 1 h in some cases. Number of thruster firings can also be reduced by 31 instances over the 6-h $t_{obs}$ . In general, axial fuel usage was high while lateral fuel usage was low for the optimal scenarios versus the worst case scenarios.

Table 5

Station-keeping simulation case study results with optimal halo orbit phasing and mission time using both case 1 (left) and case 2 (right) keepout conditions. Stars with one column had the same results for both keepout conditions. We also tabulate the difference in using the optimal versus worst halo orbit phasing and mission time combinations for each metric. The operators ⟨⟩obs represent averages over a 6-h observation for a given mission time.

	Beta pic	GJ 832	51 Eridani	GJ 179	47 UMa	HD 219143
$Δ t_{P, opt}$ (d)	140	20	140	150	20	140
$t_{opt}$ (d)	230	210	180	190	40	130
$Δ t_{P, worst}$ (d)	30	0	150	150	80	120	130
$t_{worst}$ (d)	310	170	210	210	280	240	230
Optimal versus worst	Case 1 + 2	Case 1 + 2	Case 1 + 2	Case 1 + 2	Case 1 + 2	Case 1	Case 2
${⟨ Δ t_{D} ⟩}_{obs}$ (min)	$+ 43 (+ 301 %)$	$+ 46 (+ 405 %)$	$+ 48 (+ 450 %)$	$+ 51 (+ 490 %)$	$+ 68 (+ 644 %)$	$+ 49 (+ 472 %)$	$+ 48 (+ 448 %)$
$N$	$- 20 (- 77 %)$	$- 26 (- 81 %)$	$- 28 (- 82 %)$	$- 29 (- 83 %)$	$- 31 (- 89 %)$	$- 29 (- 83 %)$	$- 28 (- 82 %)$
$ζ_{\max}$ (km)	$+ 1 (+ 525 %)$	$+ 1 (+ 4578 %)$	$+ 1 (+ 6857 %)$	$+ 1 (+ 2577 %)$	$+ 1 (+ 1028 %)$	$+ 0 (+ 387 %)$	$+ 0 (+ 310 %)$
${⟨ Δ v_{L} ⟩}_{obs}$ (mm/s)	$- 3 (- 50 %)$	$- 17 (- 82 %)$	$- 19 (- 83 %)$	$- 19 (- 83 %)$	$- 20 (- 87 %)$	$- 21 (- 90 %)$	$- 20 (- 89 %)$
${⟨ Δ v_{A} ⟩}_{obs}$ (mm/s)	$+ 57 (+ 419 %)$	$+ 39 (- 3581 %)$	$+ 38 (- 5093 %)$	$+ 55 (- 2793 %)$	$+ 43 (- 803 %)$	$+ 42 (- 462 %)$	$+ 44 (- 409 %)$
${⟨ {\dot{m}}_{L} ⟩}_{obs}$ (kg/d)	$- 2 (- 88 %)$	$- 9 (- 97 %)$	$- 11 (- 97 %)$	$- 11 (- 97 %)$	$- 11 (- 98 %)$	$- 11 (- 98 %)$	$- 11 (- 98 %)$
${⟨ {\dot{m}}_{A} ⟩}_{obs}$ (kg/d)	$+ 1 (+ 20 %)$	$+ 4 (- 753 %)$	$+ 4 (- 981 %)$	$+ 6 (- 562 %)$	$+ 5 (- 180 %)$	$+ 8 (- 162 %)$	$+ 8 (- 154 %)$
${⟨ f_{p, L} ⟩}_{obs}$ (%)	$- 0.2 (- 88 %)$	$- 0.7 (- 97 %)$	$- 0.9 (- 97 %)$	$- 0.9 (- 97 %)$	$- 0.9 (- 98 %)$	$- 0.9 (- 98 %)$	$- 0.8 (- 98 %)$
$⟨ f_{p, A} ⟩_{obs}$ (%)	$+ 0.1 (+ 20 %)$	$+ 0.3 (- 753 %)$	$+ 0.3 (- 981 %)$	$+ 0.4 (- 562 %)$	$+ 0.4 (- 180 %)$	$+ 0.6 (- 162 %)$	$+ 0.6 (- 154 %)$

With the same simulation results used in Table 5, we also studied how the selection of an optimal halo phasing affects the metrics over an entire mission. We define the operator ${⟨ ⟩}_{t}$ to denote the average value of a metric throughout the mission life span. For a given metric $M$ and average values ${⟨ M_{k} ⟩}_{obs}$ for each specific mission time $k$ , we calculate the double-averaged metric using the given operator as

Eq. (82)

{⟨ {⟨ M ⟩}_{obs} ⟩}_{t} = \frac{1}{K} \sum_{k = 0}^{K} {⟨ M_{k} ⟩}_{obs},

where

K

is the number of mission times specified. We take the resulting average drift times

{⟨ Δ t_{D} ⟩}_{obs}

over the mission time

t

to find

{⟨ {⟨ Δ t_{D} ⟩}_{obs} ⟩}_{t}

. We then define the optimal halo orbit phasing

Δ t_{P, opt}

as the corresponding halo orbit phasing that maximizes

{⟨ {⟨ Δ t_{D} ⟩}_{obs} ⟩}_{t}

. With this definition, we attempt to find a phasing that globally optimizes all possible observations of the star. We also calculate the least optimal halo orbit phasing and recorded the differences in Table 6. We see similar trends to those shown through scheduling optimization in Sec. 7.1: choosing the optimal halo orbit phasing position at the start of the mission improves almost every metric over the mission lifespan. We see larger improvements again for mid-latitude stars, though the difference is less pronounced than for the scheduling optimization. Average drift times can be increased by up to 21 min in the case of 47 UMa with 14 fewer interruptions through this strategy. As with the scheduling optimization, this strategy can be performed for any target list of stars. It is important to note that the selection of a specific halo orbit phasing, though it may optimize observation conditions for some targets, may create unfavorable conditions for other targets. Mission planners should then determine a halo orbit phasing at mission start that will produce the most favorable observations throughout an entire mission for high-priority targets.

Table 6

Station-keeping simulation case study results with optimal halo orbit phasing using both case 1 (left) and case 2 (right) keepout conditions. Stars with one column had the same results for both keepout conditions. We also tabulate the difference in using the optimal versus worst halo orbit phasing times for each metric. The operators ⟨⟩obs signify averages over a 6-h observation for a given mission time and ⟨⟩t averages over all mission times.

	Beta Pic	GJ 832	51 Eridani	GJ 179		47 UMa	HD 219143
$Δ t_{P, opt}$ (d)	90	20	130	130		30	150
$Δ t_{P, worst}$ (d)	0	140	0	20	10	120	60
Optimal versus worst	Case 1 + 2	Case 1 + 2	Case 1 + 2	Case 1	Case 2	Case 1 + 2	Case 1 + 2
${⟨ {⟨ Δ t_{D} ⟩}_{obs} ⟩}_{t}$ (min)	$+ 12 (+ 60 %)$	$+ 11 (+ 75 %)$	$+ 12 (+ 83 %)$	$+ 12 (+ 78 %)$	$+ 12 (+ 78 %)$	$+ 21 (+ 150 %)$	$+ 11 (+ 67 %)$
$⟨ N ⟩_{t}$	$- 6 (- 33 %)$	$- 7 (- 28 %)$	$- 9 (- 34 %)$	$- 8 (- 32 %)$	$- 8 (- 31 %)$	$- 14 (- 53 %)$	$- 4 (- 18 %)$
${⟨ ζ_{\max} ⟩}_{t}$ (km)	$+ 0.2 (+ 316 %)$	$+ 0.0 (+ 43 %)$	$+ 0.0 (+ 24 %)$	$- 0.0 (- 4 %)$	$- 0.0 (- 7 %)$	$+ 0.1 (+ 191 %)$	$+ 0.1 (+ 492 %)$
${⟨ {⟨ Δ v_{L} ⟩}_{obs} ⟩}_{t}$ (mm/s)	$- 4 (- 36 %)$	$- 5 (- 30 %)$	$- 5 (- 31 %)$	$- 5 (- 33 %)$	$- 5 (- 32 %)$	$- 10 (- 56 %)$	$- 3 (- 19 %)$
${⟨ {⟨ Δ v_{A} ⟩}_{obs} ⟩}_{t}$ (mm/s)	$+ 19 (+ 285 %)$	$+ 13 (- 721 %)$	$+ 13 (- 626 %)$	$+ 19 (- 205 %)$	$+ 17 (- 242 %)$	$+ 21 (- 367 %)$	$+ 7 (+ 262 %)$
$⟨ ⟨ {\dot{m}}_{L} ⟩_{obs} ⟩_{t}$ (kg/d)	$- 2 (- 55 %)$	$- 2 (- 43 %)$	$- 3 (- 48 %)$	$- 3 (- 49 %)$	$- 3 (- 48 %)$	$- 5 (- 75 %)$	$- 1 (- 16 %)$
${⟨ {⟨ {\dot{m}}_{A} ⟩}_{obs} ⟩}_{t}$ (kg/d)	$+ 2 (+ 131 %)$	$+ 2 (- 396 %)$	$+ 3 (- 390 %)$	$+ 4 (- 148 %)$	$+ 4 (- 167 %)$	$+ 4 (- 206 %)$	$- 0 (- 13 %)$
${⟨ {⟨ f_{p, L} ⟩}_{obs} ⟩}_{t}$ (%)	$- 0.1 (- 55 %)$	$- 0.2 (- 43 %)$	$- 0.2 (- 48 %)$	$- 0.2 (- 49 %)$	$- 0.2 (- 48 %)$	$- 0.4 (- 75 %)$	$- 0.1 (- 16 %)$
${⟨ {⟨ f_{p, A} ⟩}_{obs} ⟩}_{t}$ (%)	$+ 0.2 (+ 131 %)$	$+ 0.2 (- 396 %)$	$+ 0.2 (- 390 %)$	$+ 0.4 (- 148 %)$	$+ 0.3 (- 167 %)$	$+ 0.3 (- 206 %)$	$- 0.0 (- 13 %)$

8. Conclusions

In this paper, we presented a comprehensive analytical model for starshade formation flying as a function of known quantities. This model provides a clear framework for the fundamental geometric constraints of the telescope-starshade configuration A clear derivation of differential forces also helped bridge the gap between the starshade orbit geometry and JPL S5 work on station-keeping. Running simulations with publicly available Python packages, though at somewhat lower fidelity, helps bring station-keeping results to full end-to-end mission simulators like EXOSIMS. We present a number of metrics that, after generating them offline and storing them in interpolate-able look-up tables, can be used for both target selection and future scheduling of observations within these mission simulations for the first time. Finally, we highlight the importance of halo orbit phasing for optimizing observation scheduling. When selecting target stars known to harbor exoplanets, it is vital to allot as much time as possible for detection and spectral characterization. Our methods maximize the amount of time these targets are visible to the telescope and in fuel-conserving geometries, ensuring high-quality observations. These configurations are pivotal for mission design of observations with starshades. Integration times may also take multiple weeks to fully characterize the exoplanet; fuel costs can easily accumulate over these long observations. For these reasons, it is important to know formation flying costs beforehand. We can simulate costs to fuel usage and observation interruptions with EXOSIMS and use them to predict optimal times for target observation. We can also vary halo orbit phasings so that high-priority stars are visible for a bigger fraction of the mission due to keepout constraints. If we can ensure minimal fuel usage—with as much uninterrupted observation time as possible over a mission lifetime—we can increase the number of observations that can be performed as well as our chances of successful exoplanet detections and characterizations.

Acknowledgments

This work was supported by NASA JPL SURP grant RSA No. 1618976. We would also like to thank Thibault Flinois for helpful guidance regarding the JPL S5 Team’s deadbanding strategies.

References

1.

D. J. DeMarais and M. R. Walter, “Astrobiology: exploring the origins, evolution, and distribution of life in the universe,” Annu. Rev. Ecol. Evol. Syst., 30 397 –420 (1999). https://doi.org/10.1146/annurev.ecolsys.30.1.397 1543-592X Google Scholar

2.

T. D. Robinson et al., “Earth as an extrasolar planet: Earth model validation using epoxi Earth observations,” Astrobiology, 11 (5), 398 –408 (2011). https://doi.org/10.1089/ast.2011.0642 ASTRC4 1531-1074 Google Scholar

3.

M. C. Noecker, “Alignment of a terrestrial planet finder starshade at 20-100 megameters,” Proc. SPIE, 6693 669306 (2007). https://doi.org/10.1117/12.736053 PSISDG 0277-786X Google Scholar

4.

S. Seager et al., “The Exo-S probe class starshade mission,” Proc. SPIE, 9605 96050W (2015). https://doi.org/10.1117/12.2190378 PSISDG 0277-786X Google Scholar

5.

S. B. Shaklan et al., “Error budgets for the Exoplanet Starshade (Exo-S) probe-class mission study,” Proc. SPIE, 9605 96050Z (2015). https://doi.org/10.1117/12.2190384 PSISDG 0277-786X Google Scholar

6.

L. M. Palacios, A. Harness and N. J. Kasdin, “Hardware demonstration of starshade formation flying sensing and control algorithms,” Proc. SPIE, 11117 111170N (2019). https://doi.org/10.1117/12.2528437 PSISDG 0277-786X Google Scholar

7.

M. Bottom et al., “Starshade formation flying. I: Optical sensing,” J. Astron. Telesc. Instrum. Syst., 6 (1), 015003 (2020). https://doi.org/10.1117/1.JATIS.6.1.015003 Google Scholar

8.

Y. Kim et al., “Optical demonstration of a starshade at flight Fresnel numbers,” Proc. SPIE, 10400 104001A (2017). https://doi.org/10.1117/12.2273287 PSISDG 0277-786X Google Scholar

9.

T. L. B. Flinois et al., “Starshade formation flying. II: Formation control,” J. Astron. Telesc. Instrum. Syst., 6 029001 (2020). https://doi.org/10.1117/1.JATIS.6.2.029001 Google Scholar

10.

G. M. Kuan et al., “Overview of the 4m baseline architecture concept of the habitable exoplanet imaging mission (HabEx) study,” Proc. SPIE, 10698 106980Q (2018). https://doi.org/10.1117/12.2314241 PSISDG 0277-786X Google Scholar

11.

B. S. Gaudi et al., “The habitable exoplanet observatory (HabEx) mission concept study final report,” (2020). Google Scholar

12.

D. Sirbu et al., “Stationkeeping for an occulter-based exoplanetary imaging mission,” in Guidance and Control 2010—Adv. Astron. Sci., 33rd Annu. AAS Rocky Mount. Guidance and Control Conf., 129 –148 (2010). Google Scholar

13.

G. Soto et al., “Parameterizing the search space of starshade fuel costs for optimal observation schedules,” J. Guidance, Control Dyn., 42 2671 –2676 (2019). https://doi.org/10.2514/1.G003747 Google Scholar

14.

T. Flinois et al., S5: Starshade Technology to TRL5 Milestone 4 Final Report: Lateral Formation Sensing and Control, Jet Propulsion Laboratory Publications(2018). Google Scholar

15.

P. Virtanen et al., “SciPy 1.0: Fundamental algorithms for scientific computing in Python,” Nat. Methods, 17 261 –272 (2020). https://doi.org/10.1038/s41592-019-0686-2 Google Scholar

16.

D. Savransky and D. Garrett, “WFIRST-AFTA coronagraph science yield modeling with EXOSIMS,” J. Astron. Telesc. Instrum. Syst., 2 011006 (2016). https://doi.org/10.1117/1.JATIS.2.1.011006 Google Scholar

17.

R. M. Morgan et al., “Standard exoplanet yield evaluation for the LUVOIR and HabEx concept studies,” Proc. SPIE, 11117 1111701 (2019). https://doi.org/10.1117/12.2530668 PSISDG 0277-786X Google Scholar

18.

R. A. Brown, “Single-visit photometric and obscurational completeness,” Astrophys. J., 624 1010 –1024 (2005). https://doi.org/10.1086/429124 ASJOAB 0004-637X Google Scholar

19.

D. Garrett and D. Savransky, “Analytical formulation of the single-visit completeness joint probability density function,” Astrophys. J., 828 20 (2016). https://doi.org/10.3847/0004-637X/828/1/20 ASJOAB 0004-637X Google Scholar

20.

D. Savransky, N. J. Kasdin and E. Cady, “Analyzing the designs of planet-finding missions,” Publ. Astron. Soc. Pac., 122 (890), 401 –419 (2010). https://doi.org/10.1086/652181 PASPAU 0004-6280 Google Scholar

21.

E. Kolemen and N. J. Kasdin, “Optimization of an occulter-based extrasolar-planet-imaging mission,” J. Guidance Control Dyn., 35 (1), 172 –185 (2012). https://doi.org/10.2514/1.53479 Google Scholar

22.

W. S. Koon et al., Dynamical Systems, the Three-Body Problem and Space Mission Design, 1.2 ed.Marsden Books(2011). Google Scholar

23.

Astropy Collaboration, “The Astropy Project: building an open-science project and status of the v2.0 Core Package,” Astron. J., 156 123 (2018). https://doi.org/10.3847/1538-3881/aabc4f ANJOAA 0004-6256 Google Scholar

24.

J. D. Giorgini et al., “JPL’s on-line solar system data service,” in AAS/Div. Planet. Sci. Meeting , (1996). Google Scholar

25.

R. W. Farquhar, “The utilization of halo orbits in advanced lunar operations,” NASA Tech. Note, 1 1 –101 (1971). NASCA3 0077-3131 Google Scholar

26.

K. C. Howell, “Three-dimensional, periodic, ‘halo’ orbits,” Celest. Mech., 32 53 –71 (1984). https://doi.org/10.1007/BF01358403 CLMCAV 0008-8714 Google Scholar

27.

K. Howell and H. Pernicka, “Numerical determination of Lissajous trajectories in the restricted three-body problem,” Celest. Mech., 41 107 –124 (1988). https://doi.org/10.1007/BF01238756 CLMCAV 0008-8714 Google Scholar

28.

D. C. Folta and C. Webster, “Transfer trajectory options for servicing Sun-Earth-Moon libration point missions,” in Proc. 29th AAS/AIAA Space Flight Mech. Meeting, (2019). Google Scholar

29.

C. McInnes, Solar Sailing: Technology, Dynamics and Mission Applications, Springer-Praxis Books(1999). Google Scholar

30.

B. Dachwald et al., “Parametric model and optimal control of solar sails with optical degradation,” J. Guidance Control Dyn., 29 (5), 1170 –1178 (2006). https://doi.org/10.2514/1.20313 Google Scholar

31.

T. Glassman et al., “Creating optimal observing schedules for a starshade planet-finding mission,” in IEEE Aerosp. Conf. Proc., (2011). Google Scholar

32.

D. Sirbu, C. Vad Karsten and N. J. Kasdin, “Dynamical performance for sciencemode stationkeeping with an external occulter,” Proc. SPIE, 7731 773152 (2010). https://doi.org/10.1117/12.856361 PSISDG 0277-786X Google Scholar

33.

M. J. L. Turner, Rocket and Spacecraft Propulsion: Principles, Practice and New Developments, 3rd ed.Springer-Praxis Books(2009). Google Scholar

Biography

Gabriel J. Soto received his PhD in aerospace engineering from Cornell University in December 2020. His research has focused on fuel cost heuristics, formation flying, and trajectory design of spacecraft—notably starshades—in astrophysics and exoplanet direct imaging missions. He is now a postdoctoral research associate at the University of Wisconsin-Madison studying the optimization of integrated nuclear and renewable energy system design and operation.

Dmitry Savransky is an assistant professor at Sibley School of Mechanical and Aerospace Engineering of Cornell University, where he runs the Space Imaging and Optical Systems Laboratory. His research focuses on the application of estimation and computer vision techniques to control and autonomously operate advanced optical systems and on the design and optimization of astrophysics-focused space missions.

Rhonda Morgan received her BS degree in electrical engineering from Caltech and her MS and PhD degrees in optical sciences from the University of Arizona. She is a senior optical engineer at the Jet Propulsion Laboratory of California Institute of Technology. Her interests include exoplanets, high-contrast imaging, and advanced technologies for space telescopes. She is the HabEx technologist and lead of the Exoplanet Program Office Yield Standards Team.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Gabriel J. Soto, Dmitry Savransky, and Rhonda Morgan "Analytical model for starshade formation flying with applications to exoplanet direct imaging observation scheduling," Journal of Astronomical Telescopes, Instruments, and Systems 7(2), 021209 (11 January 2021). https://doi.org/10.1117/1.JATIS.7.2.021209

Received: 20 July 2020; Accepted: 16 December 2020; Published: 11 January 2021

Access the abstract

JOURNAL ARTICLE
32 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY