thesis/LaTeX/trajectory_design.tex

\chapter{Interplanetary Trajectory Design} \label{traj_dyn}

	In order to optimize trajectories in interplanetary space, it is necessary to understand and
	utilize the dynamical systems that will be acting on the spacecraft throughout the trajectory.
	This section will explore the models available for understanding the natural motion of a
	spacecraft in interplanetary space, beginning with the Two Body Problem and a method of patched
	conics for combining the gravitational effects of multiple primary bodies.

	\section{System Dynamics}

		\subsection{The Two-Body Problem}

			The motion of a spacecraft in space is governed by a complex dynamical environment.
			However, in the process of designing a trajectory, we often have to compute the
			path of the spacecraft many hundreds, thousands, or even millions of times. Utilizing
			very high-fidelity force models that account for aerodynamic pressure, solar radiation
			pressure, multi-body effects, and other forces may be too time intensive for a
			particular application. Initial surveys of the solution space often don't require such
			complex models in order to gain valuable preliminary insight.

			A common approach (and the one utilized in this implementation) is to first use a
			lower-fidelity dynamical model that captures only the gravitational force due to the
			primary body around which the spacecraft is orbiting. This approach can provide an
			excellent low-to-medium fidelity model that is useful as an underlying model in an
			algorithm for quickly categorizing a search space for initial mission feasibility
			explorations.

			In order to explore the Two Body Problem, we must first examine the full set of
			assumptions associated with the force model\cite{vallado2001fundamentals}. Firstly, we
			are only concerned with the gravitational influence between the nominative two bodies:
			the spacecraft and the planetary body around which it is orbiting. Secondly, both of
			these bodies are modeled as point masses with constant mass. This removes the need to
			account for non-uniform densities and asymmetry. Finally, for convenience in notation at
			the end, we'll also assume that the mass of the spacecraft ($m_2$) is much much smaller
			than the mass of the planetary body ($m_1$) and enough so as to be considered
			negligible.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.65\textwidth]{LaTeX/fig/2bp}
				\caption{Figure representing the positions of the bodies relative to each other and
				the center of mass in the two body problem}
				\label{2bp_fig}
			\end{figure}

			Under these assumptions, the force acting on the body due to the law of universal
			gravitation is:
			\begin{align}
				F_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\
				F_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|}
			\end{align}

			And by Newton's second law (force is the product of mass and acceleration), we can
			derive the following differential equations for $r_1$ and $r_2$:
			\begin{align}
				m_2 \ddot{\vec{r}}_2 &= - \frac{G m_1 m_2}{r^2} \frac{\vec{r}}{\left| r \right|} \\
				m_1 \ddot{\vec{r}}_1 &= \frac{G m_2 m_1}{r^2} \frac{\vec{r}}{\left| r \right|}
			\end{align}

			Where $\vec{r}$ is the position of the spacecraft relative to the primary body,
			$\vec{r}_1$ is the position of the primary body relative to the origin of the inertial
			frame, and $\vec{r}_2$ is the position of the spacecraft relative to the center of the
			inertial frame. $G$ is the universal gravitational parameter, $m_1$ is the mass of the
			planetary body, and $m_2$ is the mass of the spacecraft. From these equations, we can
			then determine the acceleration of the spacecraft relative to the planet:
			\begin{equation}
				\ddot{\vec{r}} = \ddot{\vec{r}}_2 - \ddot{\vec{r}}_1 =
				- \frac{G \left( m_1 + m_2 \right)}{r^2} \frac{\vec{r}}{\left| r \right|}
			\end{equation}

			Further assuming that the mass of the spacecraft is significantly smaller than the mass
			of the primary body ($m_1 >> m_2$) we can simplify the problem by removing the
			negligible $m_2$ term. We can also introduce, for convenience, a gravitational parameter
			$\mu$ which represents the gravity constant for the system about the center of motion
			($\mu = G (m_1 + m_2) \approx G m_1$). Doing so and simplifying produces:
			\begin{equation}
				\ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r}
			\end{equation}

			Since the spacecraft is acting only under the gravitational influence of the planet and
			no other forces, we can define the total specific mechanical energy as
			\cite{vallado2001fundamentals}:
			\begin{equation} \label{energy}
				\xi = \frac{v^2}{2} - \frac{\mu}{r}
			\end{equation}
			Where the first term represents the kinetic energy of the spacecraft and the second term
			represents the gravitational potential energy.

			Now that we've fully qualified the forces acting within the Two Body Problem, we can concern
			ourselves with more practical applications of it as a force model. It should be noted,
			firstly, that the spacecraft's position and velocity (given an initial position and velocity
			and of course the $\mu$ value of the primary body) is actually analytically solvable for all
			future points in time. This can be easily observed by noting that there are three
			one-dimensional equations (one for each component of the three-dimensional space) and
			three unknowns (the three components of the second derivative of the position).

		\subsection{Kepler's Laws}

			In the early 1600s, Johannes Kepler produced just such a solution, by taking advantages of
			what is also known as ``Kepler's Laws'' which are\cite{murray1999solar}:

			\begin{enumerate}
				\item Each planet's orbit is an ellipse with the Sun at one of the foci. This can be
					expanded to any orbit by re-wording as ``all orbital paths follow a conic section
					(circle, ellipse, parabola, or hyperbola) with a primary mass at one of the foci''.

					The conic trajectory equation explains this observation and offers a description
					of the path as:
					\begin{equation}
						r = \frac{\sfrac{h^2}{\mu}}{1 + e \cos(\theta)}
					\end{equation}
					where $h$ is the angular momentum of the satellite, $e$ is the
					eccentricity of the orbit, and $\theta$ is the true anomaly, or simply
					the angular distance the satellite has traversed along the orbit path from
					periapsis.

				\item The area swept out by the imaginary line connecting the primary and secondary
					bodies increases linearly with respect to time. This implies that the magnitude of the
					orbital speed is not constant. For the moment, we'll just take this
					value to be a constant:
					\begin{equation}\label{swept}
						\frac{\Delta t}{T} = \frac{k}{\pi a b}
					\end{equation}
					where $k$ is the constant value, $a$ and $b$ are the semi-major and
					semi-minor axis of the conic section, and $T$ is the period. In the
					following section, we'll derive the value for $k$.

				\item The square of the orbital period is proportional to the cube of the semi-major
					axis of the orbit, regardless of eccentricity. For an elliptical orbit this
					observation connects to the following known expression for the orbit period:
					\begin{equation}
						T = 2 \pi \sqrt{\frac{a^3}{\mu}}
					\end{equation}
					where $T$ is the period and $a$ is the semi-major axis.
			\end{enumerate}

		\subsection{Kepler's Equation}

			Kepler was able to produce an equation to represent the angular displacement of an
			orbiting body around a primary body as a function of time, which we'll derive now for
			the elliptical case\cite{vallado2001fundamentals}. Because the total area of an ellipse
			is the product of $\pi$, the semi-major axis, and the semi-minor axis ($\pi a b$), we
			can relate (by Kepler's second law) the area swept out by an orbit as a function of
			time, as we did in Equation~\ref{swept}. This leaves just one unknown variable $k$,
			which we can determine through use of the geometric auxiliary circle, which is a circle
			with radius equal to the ellipse's semi-major axis and center directly between the two
			foci, as in Figure~\ref{aux_circ}.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.8\textwidth]{LaTeX/fig/kepler}
				\caption{Geometric representation of auxiliary circle}\label{aux_circ}
			\end{figure}

			In order to find the area swept by the spacecraft\cite{vallado2001fundamentals}, $k$, we
			can take advantage of the fact that that area is the triangle $k_1$ subtracted from the
			elliptical segment $PCB$:
			\begin{equation}\label{areas_eq}
				k = area(seg_{PCB}) - area(k_1)
			\end{equation}
			Where the area of the triangle $k_1$ can be found easily using geometric formulae:
			\begin{align}
				area(k_1) &= \frac{1}{2} \left( ae - a \cos E \right) \left( \frac{b}{a} a \sin E \right) \\
									&= \frac{ab}{2} \left(e \sin E - \cos E \sin E \right)
			\end{align}

			\noindent
			Where $E$, notably, is the eccentric anomaly, or the angular distance from the
			periapsis to the vertical projection of the spacecraft on the auxiliary circle. Now we
			can find the area for the elliptical segment $PCB$ by first finding the circular segment
			$POB'$, subtracting the triangle $COB'$, then applying the fact that an ellipse is
			merely a vertical scaling of a circle by the amount $\frac{b}{a}$.
			\begin{align}
				area(PCB) &= \frac{b}{a} \left( area(POB') - area(COB') \right) \\
				&= \frac{b}{a} \left( \frac{a^2 E}{2} - \frac{1}{2} \left( a \cos E \right)
				\left( a \sin E \right) \right) \\
				&= \frac{abE}{2} - \frac{ab}{2} \left( \cos E \sin E \right) \\
				&= \frac{ab}{2} \left( E - \cos E \sin E \right)
			\end{align}

			By substituting the two areas back into Equation~\ref{areas_eq} we can get the $k$ area
			swept out by the spacecraft:
			\begin{equation}
				k = \frac{ab}{2} \left( E - e \sin E \right)
			\end{equation}
			Which we can then substitute back into the equation for the swept area as a function of
			time (Equation~\ref{swept}) for period of time since the spacecraft left periapsis:
			\begin{equation}
				\frac{\Delta t}{T} = \frac{t_2 - t_{peri}}{T} = \frac{E - e \sin E}{2 \pi}
			\end{equation}
			Which is, effectively, Kepler's equation. It is commonly known by a different form:
			\begin{equation}
				M = \sqrt{\frac{\mu}{a^3}} \Delta t = E - e \sin E
			\end{equation}
			where we've defined the mean anomaly as $M$ and used the fact that $T =
			\sqrt{\frac{a^3}{\mu}}$. This provides us a useful relationship between eccentric anomaly
			($E$) which can be related to spacecraft position, and time, but we still need a useful
			algorithm for solving this equation in order to use this equation to propagate a
			spacecraft.

	\subsection{LaGuerre-Conway Algorithm}\label{laguerre}

		For this thesis, the algorithm used to solve Kepler's equation was the general numeric
		root-finding scheme first developed by LaGuerre in the 1800s and first applied to Kepler's
		equation by Bruce Conway in 1985\cite{laguerre_conway}. In his paper, Conway makes a
		compelling argument for utilizing the less common LaGuerre method over higher order Newton
		or Newton-Raphson methods. The Newton-Raphson methods, while found to generally have quite
		impressive convergence rates (generally successfully solving Kepler's equation correctly
		within 5 iterations), were prone to failures in convergence given certain specific initial
		conditions. Therefore LaGuerre's algorithm is proposed as an alternative.

		The algorithm can be derived by examining the polynomial equation with $m$ roots:
		\begin{equation}
			g(x) = (x - x_1) (x - x_2) ... ( x - x_m)
		\end{equation}
		We can then generate some useful convenience functions as:
		\begin{align}
			\ln|g(x)| &= \ln|(x - x_1)| + \ln|(x - x_2)| + ... + \ln|( x - x_m)| \\
			\frac{d\ln|g(x)|}{dx} &= \frac{1}{x - x_1} + \frac{1}{x - x_2} + ... + \frac{1}{x -
			x_m} = G_1(x)
		\end{align}
		and
		\begin{align}
			\frac{-d^2\ln|g(x)|}{dx^2} &= \frac{1}{(x - x_1)^2} + \frac{1}{(x - x_2)^2} + ... +
			\frac{1}{(x - x_m)^2} = G_2(x)
		\end{align}

		Now we define the targeted root as $x_1$ and make the approximation that all of the
		other roots are equidistant from the targeted root, which means:
		\begin{equation}
			x - x_i = b, i=2,3,...,m
		\end{equation}
		We can then rewrite $G_1$ and $G_2$ as:
		\begin{align}
			G_1 &= \frac{1}{a} + \frac{n-1}{b} \\
			G_2 &= \frac{1}{a^2} + \frac{n-1}{b^2}
		\end{align}
		Which may be solved for $a$ in terms of $G_1$, $G_2$:
		\begin{equation}
			a = \frac{n}{G_1 \pm \sqrt{(n-1)(nG_2 - G_1^2)}}
		\end{equation}
		With corresponding iteration function:
		\begin{equation}
			x_{i+1} = x_i - \frac{n g(x_i)}{g'(x_i) \pm \sqrt{(n-1)^2 f'(x_i)^2 - n (n-1) f(x_i)
			f''(x_i)}}
		\end{equation}

		This iteration scheme can be shown to be globally convergent, regardless of the initial
		guess. Conway also showed that the application of this method to Kepler's equation was shown
		to converge with similar speed to many of the best common higher order Newton-Raphson
		solvers. However, LaGuerre's method was also found to be incredibly robust, converging to
		the correct value for every one of Conway's 500,000 tests. Because of this robustness, it is
		useful for solving Kepler's equation.

	\section{Interplanetary Trajectories}\label{interplanetary}

		In interplanetary travel, the primary body most responsible for gravitational forces might
		be a number of different bodies, dependent on the phase of the mission. In fact, at some
		points along the trajectory, there may not be a ``primary'' body, but instead a number of
		different forces of roughly equal magnitude vying for ``primary'' status.

		In the ideal case, every relevant body would be considered as an ``n-body'' perturbation
		during the entire trajectory. For some approaches, this method is sufficient and preferred.
		However, for most uses, a more efficient model is necessary. The method of patched conics
		can be applied in this case to simplify the model.

		Interplanetary travel does not merely complicate trajectory optimization. The increased
		complexity of the search space also opens up new opportunities for mission designers. The
		primary strategy investigated by this thesis will be the gravity assist, a technique for
		utilizing the gravitational energy of a planet to modify the direction of solar velocity.

		Finally, the concept of multiple gravity assists and the techniques used to visualize the
		space in which we might accomplish stringing together multiple flybys will be analyzed. A
		number of tools have been developed to assist mission designers in manually visualizing the
		search space, but some of these tools can also be leveraged by the automated optimization
		algorithm.

		\subsection{Patched Conics}

			The first hurdle to deal with in interplanetary space is the problem of reconciling
			Two-Body dynamics with the presence of multiple and varying planetary bodies. The most
			common method for approaching this is the method of patched conics
			\cite{bate2020fundamentals}. In this model, we break the interplanetary trajectory up
			into a series of smaller sub-trajectories. During each of these sub-trajectories, a
			single primary is considered to be responsible for the trajectory of the orbit, via the
			Two-Body problem.

			The transition point can be calculated a variety of ways. The most typical method is to
			calculate the gravitational force due to the two bodies separately, via the Two-Body
			models. Whichever primary is a larger influence on the motion of the spacecraft is
			considered to be the primary at that moment. In other words, the spacecraft, at that
			epoch, is within the Sphere of Influence of that primary. Generally for missions in this
			Solar System, the spacecraft is either within the Sphere of Influence of a planetary
			body or the Sun. However, there are points in the Solar System where the gravitational
			influence of two planetary bodies are roughly equivalent to each other and to the
			influence of the Sun.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.8\textwidth]{LaTeX/fig/patched_conics}
				\caption{Patched Conics Example Figure}
				\label{patched_conics_fig}
			\end{figure}

			This effectively breaks the trajectory into a series of arcs each governed by a distinct
			Two-Body problem patched together by distinct transition points. These transition points
			occur along the spheres of influence of the planets nearest to the spacecraft. A
			conceptual example of this process, labeled the method of patched conics, appears in
			Figure~\ref{patched_conics_fig}.

			Therefore, we must understand how to convert our spacecraft's state from the Sun frame
			to the planetary frame as it crosses this boundary. An elliptical orbit about the sun
			will have enough orbital energy to represent a hyperbolic orbit around the planet. So we
			first need to determine the velocity of the spacecraft relative to the planet as it
			crosses the SOI, which we can determine by subtraction \cite{vallado2001fundamentals}:
			\begin{equation}
				\vec{v}_{sc/p} = \vec{v}_{sc/sun} - \vec{v}_{planet/sun}
			\end{equation}

			Since the orbit around the planet is hyperbolic, in order to characterize the hyperbola
			we must determine the velocity of the spacecraft when it has infinite distance relative
			to the planet. Since this never occurs, a further approximation is made that the
			velocity of the spacecraft (relative to the planet) as it crosses the SOI can be modeled
			as the $\vec{v}_\infty$ of that hyperbolic arc.

			As an example, we may wish to determine the velocity relative to the planet that the
			spacecraft has at the periapsis of its hyperbolic trajectory during the flyby. This
			could be useful, perhaps, for sizing the $\Delta V$ required during the insertion stage
			of the mission if the spacecraft is intended to be captured into an elliptical orbit
			around its target planet. For a given incoming hyperbolic $\vec{v}_\infty$, we can first
			determine the specific mechanical energy of the hyperbola at infinite distance by using
			Equation~\ref{energy}:
			\begin{equation}
				\xi = \frac{v^2}{2} - \frac{\mu}{r} = \frac{v_\infty^2}{2}
			\end{equation}

			We can then leverage the conservation of energy to determine the velocity at a
			particular point, $r_{ins}$:
			\begin{align}
				\xi_{ins} &= \frac{v_{ins}^2}{2} - \frac{\mu}{r_{ins}} \\
				\xi_{ins} &= \xi_\infty = \frac{v_\infty^2}{2} \\
				v_{ins} &= \sqrt{\frac{2\mu}{r_{ins}} + v_\infty^2}
			\end{align}

		\subsection{Launch Considerations}

			For a satellite of given size, a certain amount of orbital energy can be imparted to the
			satellite by the launch vehicle. In practice, this value, for a particular mission, is
			actually determined as a parameter of the mission trajectory to be optimized. The excess
			velocity at infinity of the hyperbolic orbit of the spacecraft that leaves the Earth can
			be used to derive the launch energy. This is usually qualified as the quantity $C_3$,
			which is actually double the kinetic orbital energy with respect to the Sun, or simply
			the square of the excess hyperbolic velocity at infinity\cite{wie1998space}.

			This algorithm will assume that the initial trajectory at the beginning of the mission
			will be some hyperbolic orbit with velocity enough to leave the Earth. That initial
			$v_\infty$ will be used as a tunable parameter in the optimization routine. This allows
			the mission designer to include the launch $C_3$ in the cost function and, hopefully,
			determine the mission trajectory that includes the least initial launch energy. This can
			then be fed back into a mass-$C_3$ curve for prospective launch providers to determine
			what the maximum mass any launch provider is capable of imparting that specific $C_3$
			to.

			A similar approach is taken at the end of the trajectory. This algorithm doesn't attempt
			to exactly match the velocity of the planet. Instead, the excess hyperbolic velocity is
			also treated as a parameter that can be minimized by the cost function. If a trajectory
			is to then end in insertion, a portion of the mass budget can then be used for an
			impulsive thrust engine, which can provide a final insertion burn. This approach also
			allows flexibility for missions that might end in a flyby rather than insertion.

		\subsection{Gravity Assist Maneuvers}

			As previously mentioned, there are methods for utilizing the orbital energy of the other
			planets in the Solar System. This is achieved via a technique known as a Gravity Assist,
			or a Gravity Flyby. During a gravity assist, the spacecraft enters into the
			gravitational sphere of influence of the planet and, because of its excess velocity,
			proceeds to exit the sphere of influence. Relative to the planet, the speed of the
			spacecraft increases as it approaches, then decreases as it departs. From the
			perspective of the planet, the velocity of the spacecraft is unchanged. However, the
			planet is also orbiting the Sun.

			From the perspective of a Sun-centered frame, though, this is effectively an elastic
			collision. The overall momentum remains the same, with the spacecraft either gaining or
			losing some in the process (dependent on the directions of travel). The planet also
			loses or gains momentum enough to maintain the overall system momentum, but this amount
			is negligible compared to the total momentum of the planet. The overall effect is that
			the spacecraft arrives at the planet from one direction and, because of the influence of
			the planet, leaves in a different direction\cite{negri2020historical}.

			This can be visualized in Figure~\ref{grav_assist_fig}, which shows the bend in the
			spacecraft's velocity due to the hyperbolic arc as it passes the planet. This turns the
			direction of the spacecraft's velocity relative to the planet, which has an overall
			effect on kinetic energy that can be seen by adding the two vectors to the velocity of
			the planet relative to the sun. By passing in front of the planet or behind it (relative
			to its velocity), energy can be removed or added to the spacecraft by the maneuver.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.8\textwidth]{LaTeX/fig/flyby}
				\caption{Velocity changes during a gravity assist}
				\label{grav_assist_fig}
			\end{figure}

			This effect can be used strategically. The ``bend'' due to the flyby is actually
			tunable via the exact placement of the fly-by in the b-frame, or the frame centered at
			the planet, from the perspective of the spacecraft at $v_\infty$. By modifying the
			turning angle of this bend. In doing so, one can effectively achieve a (restricted) free
			impulsive thrust event.

		\subsection{Flyby Periapsis Altitude}

			Now that we understand gravity assists, the natural question is then how to leverage
			them for achieving certain velocity changes\cite{cho2017b}. But first, we must consider
			mathematically the effect that a gravity flyby can have on the velocity of a spacecraft
			as it orbits the Sun. Specifically, we can determine the turning angle of the bend
			mentioned in the previous section, given an excess hyperbolic velocity entering the
			planet's sphere of influence ($\vec{v}_{\infty, in}$) and a target excess hyperbolic
			velocity as the spacecraft leaves the sphere of influence ($\vec{v}_{\infty, out}$):
			\begin{equation}\label{turning_angle_eq}
				\delta = \arccos \left( \frac{\vec{v}_{\infty,in} \cdot
				\vec{v}_{\infty,out}}{|\vec{v}_{\infty,in}| |\vec{v}_{\infty,out}|} \right)
			\end{equation}

			From this turning angle, we can also determine, importantly, the periapsis of the flyby
			that we must target in order to achieve the required turning angle. The periapsis of the
			flyby, however, can provide a useful check on what turning angles are possible for a
			given flyby, since the periapsis:
			\begin{equation}\label{periapsis_eq}
				r_p = \frac{\mu}{v_\infty^2} \left[ \frac{1}{\sin\left(\frac{\delta}{2}\right)} - 1 \right]
			\end{equation}
			cannot be lower than some safe value that accounts for the radius of the planet and
			perhaps its atmosphere if applicable.

		\subsection{Multiple Gravity Assist Techniques}

			Now that we can leverage gravity flybys for their change in velocity direction, the
			final remaining question is that of stringing together flybys. How, for instance, do we
			know which planets can provide feasible flyby opportunities given a known hyperbolic
			energy leaving the previous planet?

			\subsubsection{Lambert's Problem}

				The answer comes from the application of the solution to the problem, posed by
				Johann Lambert in the 18th century\cite{blanchard1969unified}, of how to determine,
				given an initial position, a final position (the ephemeris of the two planets), and
				a time of flight between the two positions, what velocity was necessary to connect
				the two states.

				There are many algorithms developed to solve Lambert's problem, but the universal
				variable method is used here for its robustness in finding trajectories regardless
				of geometry. This method is concerned with the determination of the variables $y$
				and $A$ by a method of iterating $\psi$, which represent the square root of the
				distance traveled between the two points. These variables can then be used to build
				$f$ and $g$ functions, which can completely constrain the initial and final states.
				This problem can be solved by any root-finding method, with bisection being used
				here for its robustness given any initial guess \cite{battin1984elegant}.

				Firstly, some geometric considerations must be accounted for. For any initial
				position, $\vec{r}_1$, and final position, $\vec{r}_2$, and time of flight $\Delta
				t$, there are actually two separate transfer orbits that can connect the two points
				with paths that traverse less than one full orbit. Therefore, there are
				actually then two trajectories that can connect the points
				\cite{vallado2001fundamentals}. The first of the two will have a $\Delta \theta$ of
				less than 180 degrees, which we classify as a Type I trajectory, and the second will
				have a $\Delta \theta$ of greater than 180 degrees, which we call a Type II
				trajectory. They will also differ in their direction of motion (clockwise or
				counter-clockwise about the focus). This can be seen in Figure~\ref{type1type2},
				where both of the Lambert's solutions are presented for sample points in an orbit
				around the Sun.

				\begin{figure}[H]
					\centering
					\includegraphics[width=0.8\textwidth]{LaTeX/fig/lamberts}
					\caption{Visualization of the possible solutions to Lambert's Problem}
					\label{type1type2}
				\end{figure}

				The iteration used in this thesis will start by first calculating the change in true
				anomaly, $\Delta \theta$, as well as the cosine of this value, which can be found
				by:
				\begin{align}
					\cos (\Delta \theta) &= \frac{\vec{r}_1 \cdot \vec{r}_2}{|\vec{r}_1| |\vec{r}_2|} \\
					\Delta \theta &= \arctan(y_2/x_2) - \arctan(y_1/x_1)
				\end{align}

				The direction of motion is then chosen such that counter-clockwise orbits are
				considered, as travelling in the same direction as the planets is generally more
				efficient. Next, the variable $A$ is defined:
				\begin{equation}
					A = DM \sqrt{|r_1| |r_2| (1 - \cos(\Delta \theta))}
				\end{equation}

				A is independent of $\psi$, and therefore won't need updating as the iteration
				proceeds. Then $\psi$ is initialized to any number within its bounds
				($[-4\pi,4\pi^2]$), arbitrarily set to 0, representing a parabolic arc as a starting
				point.

				From here, the iteration loop can begin. Specifically, time of flight is calculated
				at each step and compared to the expected value. The iteration proceeds until the
				time of flight matches the expected value to within a provided tolerance. In order
				to calculate the time of flight at each step, we must first calculate some useful
				coefficients:
				\begin{equation}\label{loop_start}
					c_2 = \begin{cases}
						\frac{1-\cos(\sqrt{\psi})}{\psi} \quad &\text{if} \, \psi > 10^{-6} \\
						\frac{1-\cosh(\sqrt{-\psi})}{\psi} \quad &\text{if} \, \psi < -10^{-6} \\
						1/2 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6}
					\end{cases}
				\end{equation}
				\begin{equation}
					c_3 = \begin{cases}
						\frac{\sqrt{\psi} - \sin \sqrt{\psi}}{\psi^{3/2}} \quad &\text{if} \, \psi > 10^{-6} \\
						\frac{\sinh\sqrt{-\psi} - \sqrt{-\psi}}{(-\psi)^{3/2}} \quad &\text{if} \, \psi < -10^{-6} \\
						1/6 \quad &\text{if} \, 10^{-6} > \psi > -10^{-6}
					\end{cases}
				\end{equation}
				Where the conditions of this piecewise function represent the elliptical,
				hyperbolic, and parabolic cases, respectively. Once we have these, we can calculate
				another variable, $y$:
				\begin{equation}
					y = |r_1| + |r_2| + \frac{A (c_3 \psi - 1)}{\sqrt{c_2}}
				\end{equation}

				We can then finally calculate the variable $\chi$, and from that, the time of
				flight:
				\begin{equation}
					\chi = \sqrt{\frac{y}{c_2}}
				\end{equation}
				\begin{equation}
					\Delta t = \frac{c_3 \chi^3 + A \sqrt{y}}{\sqrt{c_2}}
				\end{equation}

				Based on the value of this time of flight and how it compares to the expected value,
				the bounds on $\psi$ are adjusted, a new $\psi$ is calculated at the midpoint
				between the bounds, and the iteration begins again at Equation~\ref{loop_start}. If
				the time of flight is sufficiently close to the expected value, the algorithm is
				allowed to complete.

				The resulting $f$ and $g$ functions (and the derivative of $g$, $\dot{g}$) can then
				be calculated:
				\begin{align}
					f &= 1 - \frac{y}{|r_1|} \\
					g &= A \sqrt{\frac{y}{\mu}} \\
					\dot{g} &= 1 - \frac{y}{|r_2|}
				\end{align}
				And from these, we can calculate the velocities of the transfer points as:
				\begin{align}
					\vec{v}_1 &= \frac{\vec{r}_1 - f \vec{r}_2}{g} \\
					\vec{v}_2 &= \frac{\dot{g} \vec{r}_2 - \vec{r}_1}{g}
				\end{align}
				Fully describing the connecting path with the specified flight time.

			\subsubsection{Planetary Ephemeris}

				Applying Lambert's problem to interplanetary travel requires knowing the positions
				of the planets in the inertial reference frame at a specific epoch. Fortunately,
				many packages have been developed for this purpose. The most commonly used for this
				is the SPICE package, developed by NASA in the 1980's. This software package, which
				has ports that are widely available in a number of languages, including Julia,
				contains many useful functions for astrodynamics.

				The primary use of SPICE in this thesis, however, was to determine the planetary
				ephemeris at a known epoch. Using the NAIF0012 and DE430 kernels, ephemeris in the
				International Celestial Reference Frame could be easily determined for a given
				epoch, provided as a decimal Julian Day since the J2000 epoch.

			\subsubsection{Porkchop Plots}

				Armed with a routine for quickly determining the outgoing velocity necessary to
				reach a planet at a given time, as well as the ephemeris of the planets in question
				at any given time, one can produce a grid of departure and arrival times between two
				planetary encounters. Within this grid, one can then plot a variety of useful
				values.

				The solution to Lambert's equation provides both the velocity vectors at departure
				and the velocity vectors at arrival. Often, these will be overlayed on the gridded
				time plots, as normalized values, or sometimes converted to characteristic energy
				$C_3$. This ``porkchop plot'' allows for a quick and concise view of what orbital
				energies are required to reach a planet at a given time from a given location, as
				well as an idea of what outgoing velocities one can expect.

				Using porkchop plots such as the one in Figure~\ref{porkchop}, mission designers can
				quickly visualize which natural trajectories are possible between planets. Using the
				fact that incoming and outgoing $v_\infty$ magnitudes must be the same for a flyby,
				a mission designer can even begin to work out what combinations of flybys might be
				possible for a given timeline, spacecraft state, and planet selection.

				\begin{figure}[H]
					\centering
					\includegraphics[width=\textwidth]{LaTeX/fig/porkchop}
					\caption{A sample porkchop plot of an Earth-Mars transfer}
					\label{porkchop}
				\end{figure}

	\section{Modeling Low Thrust Control} \label{low_thrust}

		In this section, we'll discuss the intricacies of continuous low-thrust trajectories in
		particular. There are many methods for optimizing such profiles and we'll briefly discuss
		the difference between a direct and indirect optimization of a low-thrust trajectory as well
		as introduce the concept of a control law and the notation used in this thesis for modelling
		low-thrust trajectories more simply.

		\subsection{Engine Model}

			The primary advantage of continuous thrust methods over their impulsive counterparts is
			in their fuel-efficiency in generating changes in velocity. Put specifically, all
			thrusters are capable of translating a mass flow (the rate of mass ejection from the
			thruster during operation) to a thrust imparted to the craft. Low thrust techniques
			suffer from limitations in the amount of thrust they can produce, but benefit from high
			efficiency by means of achieving that thrust by means of very low mass ejection rates.

			This efficiency is often captured in a single variable called specific impulse, often
			denoted as $I_{sp}$. We can derive the specific impulse by starting with the rocket
			thrust equation\cite{sutton2016rocket}:
			\begin{equation}
				F = \dot{m} v_e + \Delta p A_e
			\end{equation}
			Where $F$ is the thrust imparted, $\dot{m}$ is the fuel mass rate, $v_e$ is the exhaust
			velocity of the fuel, $\Delta p$ is the change in pressure across the exhaust opening,
			and $A_e$ is the area of the exhaust opening. We can then define a new variable
			$v_{eq}$, such that the thrust equation becomes:
			\begin{align}
				v_{eq} &= v_e + \frac{\Delta p A_e}{\dot{m}} \\
				F &= \dot{m} v_{eq} \label{isp_1}
			\end{align}
			And we can then take the integral of this value with respect to time to find the total
			impulse, dividing by the weight of the fuel to derive the specific impulse:
			\begin{align}
				I &= \int F dt = \int \dot{m} v_{eq} dt = m_e v_{eq} \\
				I_{sp} &= \frac{I}{m_e g_0} = \frac{m_e v_{eq}}{m_e g_0} = \frac{v_{eq}}{g_0}
			\end{align}
			Plugging Equation~\ref{isp_1} into the previous equation we can derive the following
			formula for $I_{sp}$:
			\begin{equation} \label{isp_real}
				I_{sp} = \frac{F}{\dot{m} g_0}
			\end{equation}
			Which is generally taken to be a value with units of seconds and effectively represents
			the efficiency with which a thruster converts mass to thrust.

		\subsection{Sims-Flanagan Transcription}

			In this thesis the following approach is used for modeling low-thrust paths: the
			Sims-Flanagan transcription (SFT)\cite{sims1999preliminary}. The SFT allows for
			flexibility in the trade-off between fidelity and performance, which makes it very
			useful for this sort of preliminary analysis.

			First the continuous arc is subdivided into a number ($N$) of individual consistent
			timesteps of length $\frac{tof}{N}$ where the $tof$ represents the total length of time
			for that particular mission phase. The control thrust is then applied as an impulsive
			maneuver at the center of each of these time steps. This approach can be seen visualized
			in Figure~\ref{sft_fig}.

			\begin{figure}[H]
				\centering
				\includegraphics[width=0.6\textwidth]{LaTeX/fig/sft}
				\caption{Example of an orbit raising using the Sims-Flanagan Transcription with 7
				Sub-Trajectories}
				\label{sft_fig}
			\end{figure}

			Using the SFT, it is relatively straightforward to propagate a state (in the context of
			the Two-Body Problem) that utilizes a continuous low-thrust control, without the need
			for computationally expensive numerical integration algorithms, by simply solving
			Kepler's equation (using the LaGuerre-Conway algorithm introduced in
			Section~\ref{laguerre}) $N+1$ times. First, the state is propagated to the middle of the
			first arc. Then a discontinuity is allowed in the velocity at that point and the state
			is propagated again to the middle of the next arc. That process is repeated $N-1$ times,
			and then finally, the last half-arc is propagated after applying the final velocity
			change.

			This greatly reduces the computation complexity, which is particularly useful for cases
			in which low-thrust trajectories need to be calculated many millions of times, as is the
			case in this thesis. The fidelity of the model can also be easily fine-tuned. By simply
			increasing the number of sub-arcs, one can rapidly approach a fidelity equal to a
			continuous low-thrust trajectory within the Two-Body Problem, with only
			linearly-increasing computation time\cite{sims1999preliminary}.

		\subsection{Low-Thrust Control Vector Description}

			In determining a low-thrust arc, a number of variables must be accounted for and,
			ideally, optimized. Generally speaking, this means that a control law must be determined
			for the thruster. This involves determining the appropriate direction, magnitude, and
			presence of a thrust at each point along its continuous orbit.

			\subsubsection{Angle of Thrust}

			Firstly, we can examine the direction at which to point the thrusters while they are on.
			The methods for determining this direction varies greatly depending on the particular
			control law chosen for that mission. Often, this process involves first determining a
			useful frame to think about the kinematics of the spacecraft. In this case, we'll use a
			frame often used in these low-thrust control laws: the spacecraft-centered $\hat{R}
			\hat{\theta} \hat{H}$ frame. In this frame, the $\hat{R}$ direction is the radial
			direction from the center of the primary to the center of the spacecraft. The $\hat{H}$
			hat is perpendicular to this, in the direction of orbital momentum (out-of-plane) and
			the $\hat{\theta}$ direction completes the right-handed orthonormal triad.

			This frame is useful because, for a given orbit, especially a nearly circular one, the
			$\hat{\theta}$ direction is nearly aligned with the velocity direction for that orbit at
			that moment. This allows us to define a set of two angles, which we'll call $\alpha$ and
			$\beta$, to represent the in and out of plane pointing direction of the thrusters. This
			convention is useful because, in a near-circular path, a $(0,0)$ set represents a thrust
			force more or less directly in line with the direction of the velocity, a commonly
			useful thrusting direction for most effectively increasing (or decreasing if negative)
			the angular momentum and orbital energy of the trajectory.

			Using these conventions, we can then redefine our thrust vector in terms of the $\alpha$
			and $\beta$ angles in the chosen frame:
			\begin{align}
				F_r &= F \cos(\beta) \sin (\alpha) \\
				F_\theta &= F \cos(\beta) \cos (\alpha) \\
				F_h &= F \sin(\beta)
			\end{align}

		\subsubsection{Thrust Magnitude}

			There is another variable that can be varied by the majority of electric thrusters.
			Either by controlling the input power of the thruster or the duty cycle, the thrust
			magnitude can also be varied, limited by the maximum thrust available to the thruster.
			Not all control laws allow for this fine-tuned control of the thruster.

			The approach used in this thesis does vary the magnitude of the thrust control. In
			certain cases it actually can be useful to have some fine-tuned control over the
			magnitude of the thrust. Since the optimization in this algorithm is automatic, it is
			relatively straightforward to consider the control thrust as a 3-dimensional vector in
			space limited in magnitude by the maximum thrust, which allows for that increased
			flexibility.

			For instance, we can consider the case of a simple orbit raising. Given an initial orbit
			with some eccentricity and some semi-major axis, we can define a new orbit that we'd
			like to achieve that simply has a higher semi-major axis value, regardless of the
			eccentricity of the new orbit. It is well known by analysis of the famous Hohmann
			Transfer\cite{hohmann1960attainability}, that thrusting for orbit raising is most
			effective near the periapsis of an orbit, where changes in velocity will have a higher
			impact on total orbital energy. Therefore, for a given low-thrust control law that
			allows for the presence or absence of thrusting at different efficiency cutoffs, we can
			easily come up with two different orbits, each of which achieve the same semi-major
			axis, but in two different ways at two different rates, both in time and fuel use, as
			can be seen in Figures~\ref{low_efficiency_fig} and \ref{high_efficiency_fig}.

			\begin{figure}[H]
				\centering
				\includegraphics[width=\textwidth]{LaTeX/fig/low_efficiency}
				\caption{Graphic of an orbit-raising with a low efficiency cutoff}
				\label{low_efficiency_fig}
			\end{figure}

			\begin{figure}[H]
				\centering
				\includegraphics[width=\textwidth]{LaTeX/fig/high_efficiency}
				\caption{Graphic of an orbit-raising with a high efficiency cutoff}
				\label{high_efficiency_fig}
			\end{figure}