# ~/blogon main cat orbital-mechanics-part-1-the-theory-of-orbits.md

I am currently residing in something of an interim period wherein I lack any real responsibilities. The last of my final exams has been written, but my summer job has yet to commence. This has led to my suddenly accumulating an ungodly number of hours spent playing Kerbal Space Program. Incidentally, I recently completed a course in classical mechanics which discussed some of the basics of orbital mechanics. This, combined with my recent KSP binge, was enough to inspire me to write a series of posts about the basics of orbital mechanics, with applications to KSP.

In this post, I will be discussing the basic theory of orbits. The principle concern here is with the solution of the two-body problem, which is to say the equations of motion of a pair of gravitationally interacting bodies. The argument followed in this article is synthesized from those presented in Goldstein’s Classical Mechanics and Taylor’s Classical Mechanics.

We begin by considering a pair of particles with masses $m_1$ and $m_2$, at positions $\mathbf{r}_1$ and $\mathbf{r}_2$ respectively. The gravitational potential between these particles is then clearly given by

$U(\mathbf{r}_1,\mathbf{r}_2)=-\frac{\gamma}{\|\mathbf{r}_2-\mathbf{r}_1\|}$

where $\gamma=Gm_1m_2$ is called the force constant, but this is only really a function of the displacement vector $\mathbf{r}=\mathbf{r}_2-\mathbf{r}_1$, so our expression for the potential is readily rewritten as

$U(\mathbf{r})=-\frac{\gamma}{\|\mathbf{r}\|}$

Knowing that the potential depends only on the displacement vector, it’s clear that $\mathbf{r}$ is a sensible generalized coordinate for the Lagrangian of this system. Another generalized coordinate is required, and the most sensible choice turns out to be the centre of mass vector

$\mathbf{R}=\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{m_1+m_2}=\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{M}$

where $M=m_1+m_2$ is called the total mass. Using $\mathbf{R}$ with $\mathbf{r}$, we can express the kinetic energy of the system as

$T=\frac{1}{2}\left(m_1\dot{\mathbf{r}}_1^2+m_2\dot{\mathbf{r}}_2^2\right)=\frac{1}{2}M\dot{\mathbf{R}}^2+\frac{1}{2}\frac{m_1m_2}{M}\dot{\mathbf{r}}^2$

using the equalities $\mathbf{r}_1=\mathbf{R}-\frac{m_2}{m_1+m_2}\mathbf{r}$ and $\mathbf{r}_2=\mathbf{R}+\frac{m_1}{m_1+m_2}\mathbf{r}$. The quantity $\mu=\frac{m_1m_2}{M}$ is called the reduced mass, and employing it gives the kinetic energy the further simplified form

$T=\frac{1}{2}M\dot{\mathbf{R}}^2+\frac{1}{2}\mu\dot{\mathbf{r}}^2$

with which we can write the Lagrangian of this system, using $\mathbf{R}$ and $\mathbf{r}$ as generalized coordinates, as

$\mathcal{L}=\frac{1}{2}M\dot{\mathbf{R}}^2+\frac{1}{2}\mu\dot{\mathbf{r}}^2+\frac{\gamma}{\|\mathbf{r}\|}$

This may seem daunting, but we can greatly simplify this problem using a change of coordinates. If we consider this problem in the centre of mass reference frame, wherein $\mathbf{R}=\mathbf{0}$, the Lagrangian takes the form

$\mathcal{L}=\frac{1}{2}\mu\dot{\mathbf{r}}^2+\frac{\gamma}{\|\mathbf{r}\|}$

which is clearly much simpler. While simply eliminating $\mathbf{R}$ may seem concerning at first blush, doing so allows us to focus on what is in some sense the interesting part of the orbit: the motion of the orbiting bodies relative to one another. $\mathbf{R}$ encodes the position of the two-body system as a whole, which is entirely uninteresting for the purposes of basic orbital mechanics. Before attempting to solve this system, we apply the conservation of angular momentum, noting that since

$\mathbf{L}=m_1\mathbf{r}_1\times\dot{\mathbf{r}}_1+m_2\mathbf{r}_2\times\dot{\mathbf{r}}_2=\frac{m_1m_2}{M^2}(m_2\mathbf{r}\times\dot{\mathbf{r}}+m_1\mathbf{r}\times\dot{\mathbf{r}})=\mu\mathbf{r}\times\dot{\mathbf{r}}$

is constant, it must be that $\mathbf{r}$ and $\dot{\mathbf{r}}$ lie in a plane, so the motion is planar. As such, it is perfectly valid to write $\mathbf{r}$ in the polar form $\mathbf{r}=r\hat{\mathbf{r}}$. It may at first appear that the decision to treat the motion as planar will be problematic when considering normal and anti-normal boosts, but it proves to be unproblematic for reasons that will be discussed later. For now, writing $\mathbf{r}$ in polar form allows us to rewrite our Lagrangian yet again as

$\mathcal{L}=\frac{1}{2}\mu(\dot{r}^2+r^2\dot{\varphi}^2)+\frac{\gamma}{r}$

from which we can finally obtain a pair of equations of motion. We first consider the equation of motion in $\varphi$, which is derived from the relevant Euler-Lagrange equation

$\frac{\partial\mathcal{L}}{\partial\varphi}-\frac{\text{d}}{\text{d}t}\frac{\partial\mathcal{L}}{\partial\dot{\varphi}}=-\frac{\text{d}}{\text{d}t}\mu r^2\dot{\varphi}=0$

which tells us that this quantity $\mu r^2\dot{\varphi}$ is constant. In fact, this is the magnitude of the angular momentum found earlier, so we denote it $L$. We then turn to considering the equation of motion in $r$, which is derived from the relevant Euler-Lagrange equation

$\frac{\partial\mathcal{L}}{\partial r}-\frac{\text{d}}{\text{d}t}\frac{\partial\mathcal{L}}{\partial\dot{r}}=\mu r\dot{\varphi}^2-\frac{\gamma}{r^2}-\mu\ddot{r}=0$

which we combine with the result of the $\varphi$ Euler-Lagrange equation and rearrange to obtain the differential equation

$\mu\ddot{r}=\frac{L^2}{\mu r^3}-\frac{\gamma}{r^2}$

This differential equation may seem alarmingly difficult, and it is in fact a non-trivial problem. To begin solving it, we note that the right hand side can be written as a derivative, giving us

$\mu\ddot{r}=-\frac{\text{d}}{\text{d}r}\left(\frac{1}{2}\frac{L^2}{\mu r^2}-\frac{\gamma}{r}\right)$

where both sides can then be multiplied by $\dot{r}$ to yield

$\mu\dot{r}\ddot{r}=\frac{\text{d}}{\text{d}t}\frac{1}{2}\mu\dot{r}^2$

and

$-\dot{r}\frac{\text{d}}{\text{d}r}\left(\frac{1}{2}\frac{L^2}{\mu r^2}-\frac{\gamma}{r}\right)=-\frac{\text{d}}{\text{d}t}\left(\frac{1}{2}\frac{L^2}{\mu r^2}-\frac{\gamma}{r}\right)$

which, when combined and rearranged, tell us that

$\frac{\text{d}}{\text{d}t}\left(\frac{1}{2}\mu\dot{r}^2+\frac{1}{2}\frac{L^2}{\mu r^2}-\frac{\gamma}{r}\right)=0$

meaning that $\frac{1}{2}\mu\dot{r}^2+\frac{1}{2}\frac{L^2}{\mu r^2}-\frac{\gamma}{r}$ is constant. This quantity is actually the total energy of our system, so we denote it $E$. If we rearrange this expression, we readily obtain the separable differential equation

$\frac{\text{d}r}{\text{d}t}=\sqrt{\frac{2}{\mu}\left(E+\frac{\gamma}{r}-\frac{1}{2}\frac{L^2}{\mu r^2}\right)}$

which, while soluble, yields an exceedingly complex integral. At any rate, we are more interested in knowing the shape of our orbit than we are in knowing where along it we will be at any particular moment in time, so we will attempt to find $r(\varphi)$ instead. Using the expression for $\dot{\varphi}$ found previously, and substituting it into the above differential equation, we have

$\frac{\text{d}r}{\text{d}\varphi}=\frac{\text{d}r}{\text{d}t}\frac{\text{d}t}{\text{d}\varphi}=\frac{\mu r^2}{L}\sqrt{\frac{2}{\mu}\left(E+\frac{\gamma}{r}-\frac{1}{2}\frac{L^2}{\mu r^2}\right)}$

which is once again a separable differential equation, this time suggesting the solution

$\varphi=\int\frac{\text{d}r}{\frac{\mu r^2}{L}\sqrt{\frac{2}{\mu}\left(E+\frac{\gamma}{r}-\frac{1}{2}\frac{L^2}{\mu r^2}\right)}}=\int\frac{\text{d}r}{r^2\sqrt{\frac{2\mu E}{L^2}+\frac{2\mu\gamma}{L^2r}-\frac{1}{r^2}}}$

which is simplified by making the substitution $u=\frac{1}{r}$. Doing so casts the integral into the form

$\varphi=-\int\frac{\text{d}u}{\sqrt{\frac{2\mu E}{L^2}+\frac{2\mu\gamma u}{L^2}-u^2}}$

which is a known form, and is readily looked up in a table of integrals. Looking up this integral, we find that

$\varphi=C-\arccos{\frac{\frac{L^2u}{\mu\gamma}-1}{\sqrt{1+\frac{2EL^2}{\mu\gamma^2}}}}$

which can be solved for $u$ to yield

$u=\frac{\mu\gamma}{L^2}\left(1+\sqrt{1+\frac{2EL^2}{\mu\gamma^2}}\cos\left(\varphi-C\right)\right)$

and finally the substitution $u=\frac{1}{r}$ may be applied to yield the final form of the equation of the orbit:

$r=\frac{\frac{\mu\gamma}{L^2}}{1+\sqrt{1+\frac{2EL^2}{\mu\gamma^2}}\cos\left(\varphi-C\right)}$

which (at least for $E<0$) is an ellipse, as expected. While this equation is entirely complete, it is somewhat unwieldy due to the large number of constants in play. In the interest of simplifying the expression, as well as making it easier to examine boosts, the constants $c=\frac{\mu\gamma}{L^2}$ and $\varepsilon=\sqrt{1+\frac{2EL^2}{\mu\gamma^2}}$ are typically introduced, giving the orbital equation the final form

$r=\frac{c}{1+\varepsilon\cos\left(\varphi-C\right)}$

I’m aware that this post may have been considerably more mathematically intensive than might be expected from a blog post which is ultimately about Kerbal Space Program, but we have successfully laid out the theoretical groundwork for the analysis of orbits. In subsequent posts, we will concern ourselves with some of the features of the orbital equation, such as parabolic and hyperbolic orbits and apoapsis and periapsis, followed by the effects of prograde, retrograde, radial in, radial out, normal, and anti-normal boosts on orbits. We will also discuss various methods of determining the equation of an orbit in Kerbal Space Program using information available in the game.