~/blog on 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 m1m_1 and m2m_2, at positions r1\mathbf{r}_1 and r2\mathbf{r}_2 respectively. The gravitational potential between these particles is then clearly given by

U(r1,r2)=γr2r1U(\mathbf{r}_1,\mathbf{r}_2)=-\frac{\gamma}{\|\mathbf{r}_2-\mathbf{r}_1\|}

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

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

Knowing that the potential depends only on the displacement vector, it’s clear that r\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

R=m1r1+m2r2m1+m2=m1r1+m2r2M\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=m1+m2M=m_1+m_2 is called the total mass. Using R\mathbf{R} with r\mathbf{r}, we can express the kinetic energy of the system as

T=12(m1r˙12+m2r˙22)=12MR˙2+12m1m2Mr˙2T=\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 r1=Rm2m1+m2r\mathbf{r}_1=\mathbf{R}-\frac{m_2}{m_1+m_2}\mathbf{r} and r2=R+m1m1+m2r\mathbf{r}_2=\mathbf{R}+\frac{m_1}{m_1+m_2}\mathbf{r}. The quantity μ=m1m2M\mu=\frac{m_1m_2}{M} is called the reduced mass, and employing it gives the kinetic energy the further simplified form

T=12MR˙2+12μr˙2T=\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 R\mathbf{R} and r\mathbf{r} as generalized coordinates, as

L=12MR˙2+12μr˙2+γr\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 R=0\mathbf{R}=\mathbf{0}, the Lagrangian takes the form

L=12μr˙2+γr\mathcal{L}=\frac{1}{2}\mu\dot{\mathbf{r}}^2+\frac{\gamma}{\|\mathbf{r}\|}

which is clearly much simpler. While simply eliminating R\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. R\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

L=m1r1×r˙1+m2r2×r˙2=m1m2M2(m2r×r˙+m1r×r˙)=μr×r˙\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 r\mathbf{r} and r˙\dot{\mathbf{r}} lie in a plane, so the motion is planar. As such, it is perfectly valid to write r\mathbf{r} in the polar form r=rr^\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 r\mathbf{r} in polar form allows us to rewrite our Lagrangian yet again as

L=12μ(r˙2+r2φ˙2)+γr\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

LφddtLφ˙=ddtμr2φ˙=0\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 μr2φ˙\mu r^2\dot{\varphi} is constant. In fact, this is the magnitude of the angular momentum found earlier, so we denote it LL. We then turn to considering the equation of motion in rr, which is derived from the relevant Euler-Lagrange equation

LrddtLr˙=μrφ˙2γr2μr¨=0\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

μr¨=L2μr3γr2\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

μr¨=ddr(12L2μr2γr)\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 r˙\dot{r} to yield

μr˙r¨=ddt12μr˙2\mu\dot{r}\ddot{r}=\frac{\text{d}}{\text{d}t}\frac{1}{2}\mu\dot{r}^2

and

r˙ddr(12L2μr2γr)=ddt(12L2μr2γr)-\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

ddt(12μr˙2+12L2μr2γr)=0\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 12μr˙2+12L2μr2γr\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 EE. If we rearrange this expression, we readily obtain the separable differential equation

drdt=2μ(E+γr12L2μr2)\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(φ)r(\varphi) instead. Using the expression for φ˙\dot{\varphi} found previously, and substituting it into the above differential equation, we have

drdφ=drdtdtdφ=μr2L2μ(E+γr12L2μr2)\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

φ=drμr2L2μ(E+γr12L2μr2)=drr22μEL2+2μγL2r1r2\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=1ru=\frac{1}{r}. Doing so casts the integral into the form

φ=du2μEL2+2μγuL2u2\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

φ=CarccosL2uμγ11+2EL2μγ2\varphi=C-\arccos{\frac{\frac{L^2u}{\mu\gamma}-1}{\sqrt{1+\frac{2EL^2}{\mu\gamma^2}}}}

which can be solved for uu to yield

u=μγL2(1+1+2EL2μγ2cos(φC))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=1ru=\frac{1}{r} may be applied to yield the final form of the equation of the orbit:

r=μγL21+1+2EL2μγ2cos(φC)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<0E<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=μγL2c=\frac{\mu\gamma}{L^2} and ε=1+2EL2μγ2\varepsilon=\sqrt{1+\frac{2EL^2}{\mu\gamma^2}} are typically introduced, giving the orbital equation the final form

r=c1+εcos(φC)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.