Deriving the changes in Keplerian Elements induced by small impulses

I'm trying to make derive a concise table of the effects on each Keplerian orbital element for a small impulse in the Radial, Tangential, or Perpendicular direction at periapsis or apoapsis, good to at least 2nd order in $\delta v/v_0$.

You can think of this as doing the longhand math for the question "What orbit will this cubesat get if I shoot it out of the ISS up, forward, or sideways by 7.7 m/sec?"

For a circular orbit with zero inclination:

$$v_R, \ v_T, \ v_P \ = \ 0, \ v_0, \ 0$$

$$r(t) \ = \ a$$

I've deduced the following empirical guesses from number crunching - integrating the ODEs to obtain perturbed orbits, then extracting the orbital elements:

$$\frac{\delta a}{a} \sim \left( \frac{\delta v_R}{v_0}\right)^2 + 2 \frac{\delta v_T}{v_0} + \left( \frac{\delta v_P}{v_0}\right)^2$$

$$\delta \epsilon \sim \frac{\delta v_R}{v_0} + 2 \frac{\delta v_R}{v_0} + \left( \frac{\delta v_P}{v_0}\right)^2$$

$$\delta i \sim \frac{\delta v_P}{v_0}$$

where $a$, $\epsilon$, and $i$ are the semi-major axis, eccentricity, and inclination, respectively.

I can run a bunch of cases including elliptical orbits of varying eccentricity and use the same empirical guessing to make a full set of coefficients including their eccentricity dependence, but I feel bad because I always try to do the math first if possible.

But I'm stumped how to begin those derivations. Ideally the answer would give an example how to do these impulse perturbation expansions with good-old pen and paper, along with hints how to do the rest. A link to the final results would be helpful, but a 'how-to' discussion of the strategy for doing the expansion is what I need most.

To repeat:, I'm constraining the problem to an elliptical orbit with zero initial inclination, small impulses, only at periapsis or apoapsis, in only the radial, tangential, or perpendicular directions.

• When you want to check your work, there are some numbers in my answer to this question: space.stackexchange.com/questions/12011/… – Organic Marble Feb 21 '17 at 4:33
• @OrganicMarble I've just added the 2nd paragraph (one sentence). Lots of goodies in your link, thanks! Actually this answer also has some useful math. – uhoh Feb 21 '17 at 4:46
• @Innovine Before I do numerical simulations I always try to do the analytical side as far as I can from first principles. It helps one appreciate the sim better and helps as an independent sanity check on at least parts of the calculation. I also find "historical math" fascinating. Orbit perturbation calculations were done when pens were feathers and lamps burned oil. I'm curious why you'd call trying to get a handle on perturbations a 'curious thing'. – uhoh Feb 22 '17 at 5:07
• This sort of derivation is definitely not my jam, so I'll just mention one maybe-obvious thing: by small-angle approximation, the true solution could have some sin() terms that your empirical derivations lack. – Russell Borogove Mar 15 '17 at 16:06
• @uhoh, I wholeheartedly agree here! Far too often we rely on simulations and don't realize when we've set incorrect initial conditions, or when the simulation breaks down, simply because we trust it blindly. On another note, you won't have the complete set of usual Keplerian orbital elements for what you're looking for if you're considering an equatorial orbit. This leads to an undefined RAAN and to the argument of periapsis being the sum of the pre-computed RAAN and the argument of periapsis. You may want to check the equinoctial OEs which don't go singular. – ChrisR Mar 15 '17 at 17:44

If we assume a perfect two-body problem, absent perturbations from external bodies or non-spherical gravity sources (i.e., perfect conic orbits with no precession or variation), your constraints regarding inclination are actually unnecessary, as we may, without loss of generality, examine this problem in a perifocal reference frame.

Perifocal coordinates

Let the origin of the system be the focus of the orbit, and let the orthogonal unit vectors $\mathbf{\hat{p}}$ and $\mathbf{\hat{q}}$ lie in the orbital plane, such that the periapsis of the orbit is located on a line extending from the origin in the $\mathbf{\hat{p}}$ direction. The unit vector $\mathbf{\hat{w}} = \mathbf{\hat{p}} \times \mathbf{\hat{q}}$ completes the coordinate system. (The orbital angular momentum vector lies in the $\mathbf{\hat{w}}$ direction).

If we further constrain ourselves to only look at periapsis and -- provided $e < 1$ -- apoapsis, our position vector is then strictly in the $\pm \mathbf{\hat{p}}$ direction, and our velocity vector is strictly in the $\pm \mathbf{\hat{q}}$ direction.

Keplerian elements in relation to perifocal coordinates

The Keplerian elements $i$ (inclination), $\Omega$ (argument of ascending node), and $\omega$ (argument of periapsis) are simply Euler angles used to translate between a perifocal coordinate system and an equatorial coordinate system. By focusing on a perifocal coordinate system, we can temporarily dispense with these, leaving only $e$ (eccentricity), $a$ (semi-major axis), and $\nu$ (true anomaly).

With that established, we need to develop two more scalar terms that remain constant in an unperturbed orbit (note: $\mu$ is the gravitational parameter of the central body):

• Specific energy: $$\varepsilon = \frac{v^2}{2} - \frac{\mu}{\|\mathbf{r}\|}$$
• Specific angular momentum: $$h = \|\mathbf{r} \times \mathbf{v}\|$$

To facilitate the upcoming derivations, we will expand these in terms of components (including ones known to be zero): $$\varepsilon = \frac{v^2_p + v^2_q + v^2_w}{2} - \frac{\mu}{\sqrt{r^2_p + r^2_q + r^2_w}}$$ $$h = \sqrt{(r_q v_w - r_w v_q)^2 + (r_w v_p - r_p v_w)^2 + (r_p v_q - r_q v_p)^2}$$

We now note that $$a = -\frac{\mu}{2\varepsilon}$$ and $$e = \sqrt{1 + \frac{2 \varepsilon h^2}{\mu^2}}$$

Finally, to your question

To frame this explicitly, you are looking for how the elements are affected by velocity perturbations. This means we want some expressions for these as a function of velocity, all other terms (in particular, the position) being constant.

In general, you want to accomplish this through a multivariate Taylor expansion.

Let $Q : \mathbb{R}^3 \rightarrow \mathbb{R}$ be an arbitrary real-valued continuous and at-least-twice differentiable (since you want second order accuracy) function of its arguments. Then to second order accuracy, you have $$Q(\mathbf{x}) = Q(\mathbf{a}) + \Big[(\mathbf{x}-\mathbf{a}) \cdot \mathbf{\nabla}Q(\mathbf{a})\Big] + \frac{1}{2}\Big[(\mathbf{x}-\mathbf{a}) \cdot \mathbf{H}(\mathbf{a}) \cdot (\mathbf{x}-\mathbf{a})\Big] + \epsilon(\|\mathbf{x}-\mathbf{a}\|^3)$$ where $\mathbf{\nabla}Q$ is the gradient of $Q$, i.e., $\mathbf{\nabla}Q = [Q_1 \; Q_2 \; Q_3]^T$, (using subscript numerals to denote partial differentiation with respect to that argument), and $\mathbf{H}$ is the Hessian matrix: $$H = \begin{bmatrix}Q_{11} & Q_{12} & Q_{13}\\ Q_{21} & Q_{22} & Q_{23}\\ Q_{31} & Q_{32} & Q_{33} \end{bmatrix}$$

Writing this explicitly gets cumbersome quickly and would not be easily readable, especially since it will involve multiple invocations of the chain rule at higher orders in multiple dimensions.

Instead, we will demonstrate one specific case: the semimajor axis (one of the simplest elements to calculate), $a$, to first-order accuracy only.

Let a prefix $\delta$ denote an infinitesimal perturbation to a quantity.

We now have, $$\delta a = \frac{\partial a}{\partial \varepsilon}\left(\frac{\partial \varepsilon}{\partial v_p}\delta v_p + \frac{\partial \varepsilon}{\partial v_q}\delta v_q + \frac{\partial \varepsilon}{\partial v_w} \delta v_w\right) + \epsilon(\|\delta \mathbf{v}\|^2)$$ Expanding this, we have $$\delta a = \frac{\mu}{2\varepsilon^2} \Big(v_p \delta v_p + v_q \delta v_q + v_w \delta v_w \Big) + \epsilon(\|\delta \mathbf{v}\|^2)$$

What can we draw from this? Note that the parenthesized term is identically equal to $\mathbf{v} \cdot \mathbf{\delta v}$. That is to say, to first-order accuracy, the semimajor axis is affected by the component of the perturbation in the prograde/retrograde direction only. This is the case regardless of where in the orbit this happens, not just at one of the apsides.

If we want to explore a second order expansion, it gets ugly quickly. $$\delta a = \frac{\partial a}{\partial \varepsilon}\left(\frac{\partial \varepsilon}{\partial v_p}\delta v_p + \frac{\partial \varepsilon}{\partial v_q}\delta v_q + \frac{\partial \varepsilon}{\partial v_w} \delta v_w\right) + \frac{1}{2}\left\{ \frac{\partial a}{\partial \varepsilon} \left( \frac{\partial^2 \varepsilon}{\partial v^2_p}(\delta v_p)^2 + \frac{\partial^2 \varepsilon}{\partial v_p \partial v_q}(\delta v_p \delta v_q) + \frac{\partial^2 \varepsilon}{\partial v_p \partial v_w}(\delta v_p \delta v_w) + \frac{\partial^2 \varepsilon}{\partial v_q \partial v_p}(\delta v_q \delta v_p) + \frac{\partial^2 \varepsilon}{\partial v^2_q}(\delta v_q)^2 + \frac{\partial^2 \varepsilon}{\partial v_q \partial v_w}(\delta v_q \delta v_w) + \frac{\partial^2 \varepsilon}{\partial v_w \partial v_p}(\delta v_w \delta v_p) + \frac{\partial^2 \varepsilon}{\partial v_w \partial v_q}(\delta v_w \delta v_q) + \frac{\partial^2 \varepsilon}{\partial v^2_w}(\delta v_w)^2 \right) + \frac{\partial^2 a}{\partial \varepsilon^2} \left[ \left(\frac{\partial \varepsilon}{\partial v_p}\right)^2(\delta v_p)^2 + \left(\frac{\partial \varepsilon}{\partial v_p}\right)\left(\frac{\partial \varepsilon}{\partial v_q}\right)(\delta v_p \delta v_q) + \left(\frac{\partial \varepsilon}{\partial v_p}\right)\left(\frac{\partial \varepsilon}{\partial v_w}\right)(\delta v_p \delta v_w) + \left(\frac{\partial \varepsilon}{\partial v_q}\right)\left(\frac{\partial \varepsilon}{\partial v_p}\right)(\delta v_q \delta v_p) + \left(\frac{\partial \varepsilon}{\partial v_q}\right)^2(\delta v_q)^2 + \left(\frac{\partial \varepsilon}{\partial v_q}\right)\left(\frac{\partial \varepsilon}{\partial v_w}\right)(\delta v_q \delta v_w) + \left(\frac{\partial \varepsilon}{\partial v_w}\right)\left(\frac{\partial \varepsilon}{\partial v_p}\right)(\delta v_w \delta v_p) + \left(\frac{\partial \varepsilon}{\partial v_w}\right)\left(\frac{\partial \varepsilon}{\partial v_q}\right)(\delta v_w \delta v_q) + \left(\frac{\partial \varepsilon}{\partial v_w}\right)^2(\delta v_w)^2 \right]\right\} + \epsilon(\|\delta \mathbf{v}\|^3)$$

Hoo boy!

We now make the following substitutions: $$\frac{\partial a}{\partial \varepsilon} = \frac{\mu}{2\varepsilon^2}$$ $$\frac{\partial^2 a}{\partial \varepsilon^2} = -\frac{\mu}{4\varepsilon^3}$$ $$\frac{\partial \varepsilon}{\partial v_i} = v_i,\;i=p,q,w$$ $$\frac{\partial^2 \varepsilon}{\partial v_i \partial v_j} = \delta_{ij} = \left\{\begin{array}{l}1, & i=j\\0, & i \neq j\end{array}\right., \; i,j = p,q,w$$

Now, finally

$$\delta a = \frac{\mu}{2\varepsilon^2} \Big(v_p \delta v_p + v_q \delta v_q + v_w \delta v_w \Big) + + \frac{1}{2}\left\{ \frac{\mu}{2 \varepsilon^2} \Big[ (\delta v_p)^2 + (\delta v_q)^2 + (\delta v_w)^2 \Big] - \frac{\mu}{4 \varepsilon^3} \Big[ v_p^2(\delta v_p)^2 + v_q^2(\delta v_q)^2 + v_w^2(\delta v_w)^2 + 2 v_p v_q(\delta v_p \delta v_q) + 2 v_p v_w(\delta v_p \delta v_w) + 2 v_q v_w(\delta v_q \delta v_w) \Big]\right\} +\epsilon(\|\delta \mathbf{v}\|^3)$$

Writing this more compactly: $$\delta a = \frac{\mu}{2\varepsilon^2} \left( \mathbf{v} \cdot \delta \mathbf{v} + \frac{\|\delta \mathbf{v}\|^2}{2} \right) - \frac{\mu}{4 \varepsilon^3} \Big( \mathbf{v} \cdot \delta \mathbf{v} \Big)^2 +\epsilon(\|\delta \mathbf{v}\|^3)$$

Important consequences: at second-order accuracy, the previous assertion regarding prograde/retrograde is not applicable.

Final thoughts

Though there are a few hiccups involved in applying this to the three orientation elements due to some degeneracy in the rotation formalism defined by Euler angles, this answer should provide all the tools you need to calculate what you are asking. It is a simple -- though tedious -- matter of simply following the Taylor expansion for each of the remaining elements and applying coordinate transformations where appropriate.

This should also provide a healthy illustration for why nobody really does second-order expansions. The equations get nasty very quickly, with very little benefit to offer from a numerical analysis standpoint.

• I see what you mean! I've never used perifocal coordinates. While they seem to have a clear definition here and in Wikipedia, I don't understand it - I keep thinking there's some trick I'm missing because the wording is so precisely chosen. – uhoh Mar 15 '17 at 22:10
• $\mathbf{\hat p}$ points from origin to periapsis, but what about $\mathbf{\hat q}$? It is in the orbital plane, and points in a direction that the body will be after it passes through periapsis, but is it orthogonal to $\mathbf{\hat p}$ in the real world? Are $\mathbf{\hat p, \ \hat q, \ \hat w }$ like $\mathbf{\hat x, \ \hat y, \ \hat z }$ if I were to draw them with x pointing to periapsis and z perpendicular to the plane? Or is there something going on I don't understand? – uhoh Mar 15 '17 at 22:10
• Perifocal coordinates define $\mathbf{\hat{p}}$ in the direction of periapsis, $\mathbf{\hat{w}}$ in the direction of the angular momentum vector, and $\mathbf{\hat{q}} = \mathbf{\hat{w}} \times \mathbf{\hat{p}}$ completes the right-handed orthonormal triad. – Tristan Mar 15 '17 at 22:31
• The Taylor expansion is derived around basically any arbitrary point in the orbit. The values of $\mathbf{v}$ and $\mathbf{r}$ will determine the effects. The perifocal coordinate system makes it easy to see what terms go to zero at periapsis or apoapsis. – Tristan Mar 15 '17 at 22:33
• Yes. You use precisely those numbers. I decline to pit specific numbers in the answer to avoid giving the impression that this approach is only applicable to earth orbit. The approach is general and works for any body that may be treated as a point mass. – Tristan Mar 19 '17 at 12:49

An approach could be based on analyzing energy and angular momentum changes. For example, energy can be expressed as:

$$\mathcal{E} = \frac{v^2}{2} - \frac{\mu}{r}.$$

Equivalently:

$$\mathcal{E} = \frac{\boldsymbol{v}^\mathrm{T}\boldsymbol{v}}{2} - \frac{\mu}{r}.$$

That means that the energy after an impulsive maneuver $\boldsymbol{\Delta v}$ is given by

$$\mathcal{E}_{+} = \frac{\left(\boldsymbol{v}+\boldsymbol{\Delta v}\right)^\mathrm{T}\left(\boldsymbol{v}+\boldsymbol{\Delta v}\right)}{2} - \frac{\mu}{r}.$$

Subtracting the two quantities we can obtain the total energy change

$$\mathcal{E}_+-\mathcal{E} = \Delta \mathcal{E} = \frac{\left(\boldsymbol{v}+\boldsymbol{\Delta v}\right)^\mathrm{T}\left(\boldsymbol{v}+\boldsymbol{\Delta v}\right)}{2} -\frac{\boldsymbol{v}^\mathrm{T}\boldsymbol{v}}{2}.$$

Notice that the $\mu/r$ term is eliminated because an impulsive maneuver does not change the position. After some simplification we arrive to an expression for the change in orbital energy due to an impulsive maneuver:

$$\Delta \mathcal{E} = \boldsymbol{v}^\mathrm{T}\boldsymbol{\Delta v} + \frac{\left\|\boldsymbol{\Delta v}\right\|^2}{2}.$$

Now, how can we decompose this change in radial/tangential/normal (RTN) components? First, we define the principal axes of the RTN frame:

$$\boldsymbol{i}_r = \frac{\boldsymbol{r}} {\left\|\boldsymbol{r}\right\|},\quad \boldsymbol{i}_t = \frac{\boldsymbol{i}_r\times\boldsymbol{i}_n} {\left\|\boldsymbol{i}_r\times\boldsymbol{i}_n\right\|},\quad \boldsymbol{i}_n = \frac{\boldsymbol{r}\times\boldsymbol{v}} {\left\|\boldsymbol{r}\times\boldsymbol{v}\right\|}.$$

and use them to define a rotation matrix from the RTN frame to the inertial frame:

$$\boldsymbol{Q} = \left[ \boldsymbol{i}_r\quad \boldsymbol{i}_t\quad \boldsymbol{i}_n \right]$$

With this infrastructure we can now express the change in energy due to a $\boldsymbol{\delta v}$ maneuver; I use $\delta$ instead of $\Delta$ to denote that the maneuver is expressed in the RTN frame:

$$\boldsymbol{\delta v} = \left[ \Delta v_r\quad \Delta v_t \quad \Delta v_n \right].$$

Our expression is:

$$\Delta \mathcal{E} = \boldsymbol{v}^\mathrm{T} \boldsymbol{Q}\boldsymbol{\delta v} + \frac{\left\|\boldsymbol{\delta v}\right\|^2}{2}.$$

We can use that expression to express the change in energy due to individual changes in the impulsive maneuver:

$$\Delta \mathcal{E} = \Delta \mathcal{E}_r + \Delta \mathcal{E}_t + \Delta \mathcal{E}_n$$

where, for the subindex $\gamma\in\left\{r, t, n\right\}$, each component is given by:

$$\Delta \mathcal{E}_\gamma = \boldsymbol{v}^\mathrm{T}\boldsymbol{Q} \boldsymbol{\delta v}_\gamma + \frac{\left\|\boldsymbol{\delta v}_\gamma\right\|^2}{2}.$$

The vector $\boldsymbol{\delta v}_\gamma$ has all components but $\gamma$ set to zero. You can work out the matrix product by hand if you want to.

Does this work?

Take an arbitrary orbit (canonical elements; $\mu=1$)

$$\boldsymbol{r} = \left[1, 2, 3\right],\quad \boldsymbol{v} = \left[-0.3, -0.2, -0.1\right]$$

its orbital energy is $\mathcal{E} = -1.9726124191242439\times10^{-1}$. Applying an impulsive maneuver

$$\boldsymbol{\delta v} = [1, -1, 2] \times 10^{-3}$$

leads to the new energy

$$\mathcal{E}_+ = \frac{\left\|\boldsymbol{v} + \boldsymbol{Q}\boldsymbol{\delta v}\right\|}{2} - \frac{\mu}{r} = -1.9778736462262000\times10^{-1}.$$

The change is $\Delta\mathcal{E} = -5.2612271019561528\times10^{-4}$.

Our formula correctly decomposes this change as:

\begin{align} \Delta \mathcal{E}_r &= -2.6676124191242441\times 10^{-4}\\ \Delta \mathcal{E}_t &= -2.6136146828319081\times 10^{-4}\\ \Delta \mathcal{E}_n &= +1.9999999999999321\times 10^{-6} \end{align}

which add up to $-5.2612271019561528\times10^{-4}$.

We can use those changes in energy to directly map out the changes in semimajor axis. A similar approach can be used to map out changes in angular momentum, and thus decompose eccentricity and inclination.

• This is a helpful and useful summary, and I will go through it carefully in due course. But I'm starting with a $\Delta v$ of 7.7 km/s in a 400 km orbit around the Earth, not canonical units. See this comment. – uhoh Mar 18 '17 at 6:04
• So I see that I can use real world values for velocity and $GM_E$ for $\mu$ above and so likely I can here as well. I really like this form of the derivation, and I'll go through this carefully now. There is an option in stackexchange to award a bonus to an answer ex post facto so if I get everything working I'll award a second +100 in appreciation of your time and help, and concise form of the derivation. Thanks! – uhoh Mar 21 '17 at 3:42
• This is quite an elegant and concise form. But I don't understand how to use $\Delta\mathcal{E}$ to address a change inclination. $\mathcal{E}$ is a scalar that comes from two other scalars; $v^2$ and $1/r$. Getting eccentricity and the semi-major axis is straightforward, but how identify a change in the orientation of the orbital plane itself? – uhoh Mar 26 '17 at 12:33