diff --git a/LaTeX/fig/kepler.png b/LaTeX/fig/kepler.png index 0da3b93..b6eb866 100644 Binary files a/LaTeX/fig/kepler.png and b/LaTeX/fig/kepler.png differ diff --git a/LaTeX/flowcharts/kepler.svg b/LaTeX/flowcharts/kepler.svg index 0ed6dec..15bb4a1 100644 --- a/LaTeX/flowcharts/kepler.svg +++ b/LaTeX/flowcharts/kepler.svg @@ -7,7 +7,7 @@ viewBox="0 0 279.4 215.9" version="1.1" id="svg5" - inkscape:version="1.1.1 (3bf5ae0d25, 2021-09-20, custom)" + inkscape:version="1.1.2 (0a00cf5339, 2022-02-04, custom)" sodipodi:docname="kepler.svg" inkscape:export-filename="/home/connor/projects/thesis/LaTeX/fig/kepler.png" inkscape:export-xdpi="300" @@ -34,12 +34,12 @@ inkscape:snap-grids="false" inkscape:snap-object-midpoints="false" inkscape:zoom="1.4361764" - inkscape:cx="589.76042" + inkscape:cx="589.41228" inkscape:cy="333.17634" - inkscape:window-width="747" - inkscape:window-height="1024" - inkscape:window-x="1157" - inkscape:window-y="40" + inkscape:window-width="1912" + inkscape:window-height="1040" + inkscape:window-x="0" + inkscape:window-y="32" inkscape:window-maximized="1" inkscape:current-layer="layer1" /> asin(E) + + orbit path + aux circle diff --git a/LaTeX/trajectory_design.tex b/LaTeX/trajectory_design.tex index b83f54a..aee2060 100644 --- a/LaTeX/trajectory_design.tex +++ b/LaTeX/trajectory_design.tex @@ -81,69 +81,68 @@ \ddot{\vec{r}} = - \frac{\mu}{r^2} \hat{r} \end{equation} - \subsubsection{Kepler's Laws and Equations} + \subsection{Kepler's Laws} - 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). + 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). - 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}: + 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''. + \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''. - Specifically the path of the orbit follows the trajectory equation: + Specifically the path of the orbit follows the trajectory equation: - \begin{equation} - r = \frac{\sfrac{h^2}{\mu}}{1 + e \cos(\theta)} - \end{equation} + \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. + 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. - \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: + \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} + \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$. + 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. Specifically, the relationship is: + \item The square of the orbital period is proportional to the cube of the semi-major + axis of the orbit, regardless of eccentricity. Specifically, the relationship is: - \begin{equation} - T = 2 \pi \sqrt{\frac{a^3}{\mu}} - \end{equation} + \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} + Where $T$ is the period and $a$ is the semi-major axis. + \end{enumerate} - \subsection{Analytical Solutions to Kepler's Equations} + \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}. Since 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}. + orbiting body around a primary body as a function of time, which we'll derive now for + the elliptical case\cite{vallado2001fundamentals}. Since 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 @@ -151,13 +150,15 @@ \caption{Geometric Representation of Auxiliary Circle}\label{aux_circ} \end{figure} - In order to find the area swept by the spacecraft, $k$, we can take advantage of the fact - that that area is the triangle $k_1$ subtracted from the elliptical segment $PCB$: + 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} + \noindent Where the area of the triangle $k_1$ can be found easily using geometric formulae: \begin{align} @@ -165,8 +166,11 @@ &= \frac{ab}{2} \left(e \sin E - \cos E \sin E \right) \end{align} - 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 + \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} @@ -185,10 +189,10 @@ \end{equation} Which we can then substitute back into the equation for the swept area as a function of - time (Equation~\ref{swept}): + time (Equation~\ref{swept}) for period of time since the spacecraft left periapsis: \begin{equation} - \frac{\Delta t}{T} = \frac{E - e \sin E}{2 \pi} + \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: @@ -200,18 +204,81 @@ 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. + algorithm for solving this equation in order to use this equation to propagate a + spacecraft. - \subsubsection{LaGuerre-Conway Algorithm}\label{laguerre} + \subsection{LaGuerre-Conway Algorithm}\label{laguerre} - For this application, I used an algorithm known as the LaGuerre-Conway - algorithm\cite{laguerre_conway}, which was presented in 1986 as a faster and more - robust algorithm for directly solving Kepler's equation and has been in use in many - applications since. This algorithm is known for its convergence robustness and also - its speed of convergence when compared to higher order Newton methods. + 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. - This thesis will omit a step-through of the algorithm itself, but psuedo-code for - the algorithm will be discussed briefly in Section~\ref{conway_pseudocode}. + 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 relatively easily derived by examining the polynomial equation with + $m$ roots: + + \begin{equation} + g(x) = (x - x_1) (x - x_2) ... ( x - x_m) + \end{equation} + + \noindent + 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} + + \noindent + 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} + + \noindent + 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} + + \noindent + 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. More relevantly, 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 very useful for propagating spacecraft states. \section{Interplanetary Considerations}\label{interplanetary}