This post is a largely technical note that continues an exploration of the Ideal Rocket Equation in a perturbed Keplerian gravitational field. It is a continuation of an earlier post entitled "A Generalised Ideal Rocket Equation" in the Maths & Physics section of the Orbiter Forum. These notes serve as an on-line 'aide-memoire' for me - and may be of interest to others.
Some background
The purpose of this technical note is to include the effects of off-plane thrusts and gravitational perturbations on the motion of a spacecraft executing a thrust manoeuvre - principally low thrust orbit insertion and escape burns but also, inter alia, extended plane change manoeuvres. This technical note is mainly aimed at long duration, low thrust manoeuvres - under conditions where standard Orbiter tools aren't really designed to work well.
The earlier post, "A Generalised Ideal Rocket Equation" assumed that the gravitational field was a perfect Keplerian gravitational field in which, in the absence of thrust from the spacecraft's engines, the spacecraft's motion would be described accurately by circular, elliptical, parabolic or hyperbolic motion. This assumption of Keplerian gravity - and the assumption that manoeuvres were either retrograde or prograde - meant that we could also simplify the equations of motion by assuming that all motion was in a two-dimensional orbital plane. But if we have 'off plane' thrusts (i.e. 'plane change manoeuvres') then we need to extend the mathematical apparatus to allow for these.
Moreover, close to a gravitating body, these simplifying assumptions are highly accurate, but as one moves away from the gravitating body, perturbations from other gravitating bodies (e.g., the Sun, the Moon and the planets) become more significant. Generally speaking, high thrust manoeuvres are executed close to orbital periapsis, which by virtue of the Oberth effect, usually means close to the principal gravitating body. But as available thrust becomes less, these manoeuvres start and/or end further away from that body which implies that perturbations become more significant. As a rule of thumb, the equations of motion set out in the earlier post seem to work well for thrust accelerations down to about 0.01g (c.f. the Delta Glider which produces thrust accelerations of around 2.0g). Below this, where when is applying milli-g or micro-g thrusts, the equations of motion are no longer appropriate and for long burns perturbations become increasingly important.
In this note:
A quick recap of the maths from the earlier post
The earlier post set out the maths of the basic integration scheme in some detail, but the maths was dispersed across a number of posts in the ensuing thread. So, for convenience, here in summary form are the equations of motion that describe the motion of a spacecraft executing either a retrograde or prograde burn in a pure Keplerian gravitational field.
The first three of these equations encode how three of the orbital elements - the semi-major axis [MATH]a[/MATH], the orbital eccentricity [MATH]e[/MATH], and the true anomaly [MATH]\nu[/MATH] vary as thrust is applied. (In the absence of a thrust from the spacecraft, the semi-major axis and the orbital eccentricity are constant; and the true anomaly varies in accordance with Kepler's 'equal areas in equal times' rule.) The fourth equation describes how far the spacecraft has rotated around the central gravitating body during the manoeuvre.
And the last two equations describe the radial and transversed components of the prograde or retrograde thrust. If [MATH]\nu_e>0[/MATH] the thrust is prograde; and if [MATH]nu_e<0[/MATH] the thrust is retrograde.
The boundary conditions to achieve a circular orbit are:
This system of equations does not have an analytical solution, unfortunately, so to solve this set of equations we need to use a numerical integrator. The original post set out some C code with which to perform that integration.
And now moving on to three dimensions
The preceding equations work well for prograde retrograde thrusts in the strong sphere of influence of a planet or moon - so long as the thrust is applied exclusively in the orbital plane of the spacecraft. However, some modifications to the equations of motion needed to be introduced if we want them to be more generally applicable to 'out of plane' burns - or if we wish to consider general gravitational perturbations. Specifically, we need to introduce provide some additional equations that tell us how the spacecraft's orbital plane changes during the execution of the manoeuvre.
So, whereas the first three equations of motion above provide information about the three orbital elements that specify the shape of the Keplerian orbital trajectory and the location of the spacecraft on that trajectory, an additional three equations of motion are needed to specify the orientation of the orbital plane of the trajectory on that orbital plane. Traditionally, this additional information is represented by three angle variables - the inclination, the longitude of the ascending node and the argument of periapsis.
In our case, we are also going to use three angle variables - but we are going to use a Roll-Pitch-Yaw convention in which we perform a 'roll' about the [MATH]r[/MATH]-axis of angle [MATH]\alpha(t)[/MATH], then perform a 'pitch' about the original [MATH]\theta[/MATH]-axis of angle [MATH]\beta(t)[/MATH] and, finally, perform a 'yaw' about the original [MATH]h[/MATH]-axis (i.e., normal to the trajectory plane) of angle [MATH]\gamma(t)[/MATH]. Given this, we end up with the following equations of motion for these angle variables:
where, in this case, [MATH]B_h[/MATH] is the off-plane component of applied thrust (or, more generally, gravitational perturbation). These three equations replace the single equation from the preceding two-dimensional case:
although we can see an echo of that equation in the equation for [MATH]\gamma'(t)[/MATH]. In total, then to full describe the behaviour of the trajectory of the spacecraft during a trajectory, we must solve the following six first-order differential equations:
such that:
Messy, and complicated - but true.
If we manage to solve this set of equations, then from the solution for [MATH]\alpha(t)[/MATH], [MATH]\beta(t)[/MATH] and [MATH]\gamma(t)[/MATH], we can construct a rotation matrix [MATH]\boldsymbol{U}(t)[/MATH] from the following matrix elements:
Having constructed this rotation matrix, we can calculate how the coordinates of the reference frame change with time:
And from this, we can calculate the planet-centric state vectors of the spacecraft as:
Alternatively, we can calculate the orbital elements as:
The semi-major-axis, [MATH]a[/MATH]
The orbital eccentricity, [MATH]e[/MATH]
The true anomaly, [MATH]\nu[/MATH]
(such that such if [MATH]\frac{e_{\theta}(t)}{e_{r}(t)+1}>0[/MATH] then [MATH]\nu\to2\,\pi-\nu[/MATH]).
To calculate the other three orbital elements, the inclination, the longitude of the ascending node and the argument of periapsis, we need to calculate the components of the angular momentum vector [MATH]\boldsymbol{h}(t)[/MATH] and the eccentricity vector [MATH]\boldsymbol{e}(t)[/MATH] in terms of a reference x-y-z coordinate system of a reference frame - usually either the ecliptic or equatorial reference frames. The x-y-z components of these vectors are calculated as:
and
So, now we can calculate the orbital elements in the usual way:
The inclination, [MATH]\iota[/MATH]
The longitude of the ascending node, [MATH]\Omega[/MATH]
such that if [MATH]\hat{h}_x(t)<0[/MATH] then [MATH]\Omega\to2\,\pi-\Omega[/MATH].
The argument of periapsis, [MATH]\omega[/MATH]
such that if [MATH]e_z(t)<0[/MATH] then [MATH]\omega\to 2\,\pi-\omega[/MATH].
A perturbed Keplerian gravitational field
So far, we have not worried terribly much about gravitational perturbations arising from, say, the Moon or other 'nearby' gravitating bodies. How do we include these perturbations? The formalism is very simple: we simply replace the thrust vector [MATH]\boldsymbol{B}[/MATH] with one that includes the gravitational perturbation, [MATH]\Delta\boldsymbol{P}[/MATH] such that
Some background
The purpose of this technical note is to include the effects of off-plane thrusts and gravitational perturbations on the motion of a spacecraft executing a thrust manoeuvre - principally low thrust orbit insertion and escape burns but also, inter alia, extended plane change manoeuvres. This technical note is mainly aimed at long duration, low thrust manoeuvres - under conditions where standard Orbiter tools aren't really designed to work well.
The earlier post, "A Generalised Ideal Rocket Equation" assumed that the gravitational field was a perfect Keplerian gravitational field in which, in the absence of thrust from the spacecraft's engines, the spacecraft's motion would be described accurately by circular, elliptical, parabolic or hyperbolic motion. This assumption of Keplerian gravity - and the assumption that manoeuvres were either retrograde or prograde - meant that we could also simplify the equations of motion by assuming that all motion was in a two-dimensional orbital plane. But if we have 'off plane' thrusts (i.e. 'plane change manoeuvres') then we need to extend the mathematical apparatus to allow for these.
Moreover, close to a gravitating body, these simplifying assumptions are highly accurate, but as one moves away from the gravitating body, perturbations from other gravitating bodies (e.g., the Sun, the Moon and the planets) become more significant. Generally speaking, high thrust manoeuvres are executed close to orbital periapsis, which by virtue of the Oberth effect, usually means close to the principal gravitating body. But as available thrust becomes less, these manoeuvres start and/or end further away from that body which implies that perturbations become more significant. As a rule of thumb, the equations of motion set out in the earlier post seem to work well for thrust accelerations down to about 0.01g (c.f. the Delta Glider which produces thrust accelerations of around 2.0g). Below this, where when is applying milli-g or micro-g thrusts, the equations of motion are no longer appropriate and for long burns perturbations become increasingly important.
In this note:
- I extend the maths from two dimensions to three dimensions;
- I allow for the possibility that thrusts are applied 'out of plane'; and
- I include a framework for incorporating general perturbations.
A quick recap of the maths from the earlier post
The earlier post set out the maths of the basic integration scheme in some detail, but the maths was dispersed across a number of posts in the ensuing thread. So, for convenience, here in summary form are the equations of motion that describe the motion of a spacecraft executing either a retrograde or prograde burn in a pure Keplerian gravitational field.
[MATH]h'(t) =\frac{h(t)^{2}}{\mu\left(1+e_{r}(t)\right)}\,B_{\theta}[/MATH][MATH]e_{r}'(t)-\theta'(t)\,e_{\theta}(t) =2\,\frac{h(t)}{\mu}\,B_{\theta}[/MATH][MATH]e_{\theta}'(t)+\theta'(t)\,e_{r}(t) =-\frac{h(t)}{\mu}\,B_{r}+\frac{h(t)}{\mu}\,\frac{e_{\theta}(t)}{1+e_{r}(t)}\,B_{\theta}[/MATH]
and[MATH]\theta'(t)=\frac{\mu^{2}\left(1+e_{r}(t)\right){}^{2}}{h(t)^{3}}[/MATH]
with[MATH]B_r=-\frac{\nu_{e}\gamma}{m_{0}-\gamma\,t}\,\frac{e_{\theta}(t)}{\sqrt{\left(1+e_{r}(t)\right){}^{2}+e_{\theta}(t){}^{2}}}[/MATH][MATH]B_\theta=+\frac{\nu_{e}\gamma}{m_{0}-\gamma\,t}\,\frac{1+e_{r}(t)}{\sqrt{\left(1+e_{r}(t)\right){}^{2}+e_{\theta}(t){}^{2}}}[/MATH]
and- [MATH]h(t)[/MATH] is the specific angular momentum of the spacecraft;
- [MATH]e_r(t)[/MATH] and [MATH]e_\theta(t)[/MATH] are the radial and transverse components of the eccentricity vector;
- [MATH]\theta(t)[/MATH] is the angle that the spacecraft has rotated in its orbital plane; and
- [MATH]B_r[/MATH] and [MATH]B_\theta[/MATH] are the radial and transverse components of the thrust vector.
The first three of these equations encode how three of the orbital elements - the semi-major axis [MATH]a[/MATH], the orbital eccentricity [MATH]e[/MATH], and the true anomaly [MATH]\nu[/MATH] vary as thrust is applied. (In the absence of a thrust from the spacecraft, the semi-major axis and the orbital eccentricity are constant; and the true anomaly varies in accordance with Kepler's 'equal areas in equal times' rule.) The fourth equation describes how far the spacecraft has rotated around the central gravitating body during the manoeuvre.
And the last two equations describe the radial and transversed components of the prograde or retrograde thrust. If [MATH]\nu_e>0[/MATH] the thrust is prograde; and if [MATH]nu_e<0[/MATH] the thrust is retrograde.
The boundary conditions to achieve a circular orbit are:
[MATH]e_r(t^*)=0[/MATH][MATH]e_\theta(t^*)=0[/MATH]
where [MATH]t^*[/MATH] is the time when orbital insertion manoeuvre. This system of equations does not have an analytical solution, unfortunately, so to solve this set of equations we need to use a numerical integrator. The original post set out some C code with which to perform that integration.
And now moving on to three dimensions
The preceding equations work well for prograde retrograde thrusts in the strong sphere of influence of a planet or moon - so long as the thrust is applied exclusively in the orbital plane of the spacecraft. However, some modifications to the equations of motion needed to be introduced if we want them to be more generally applicable to 'out of plane' burns - or if we wish to consider general gravitational perturbations. Specifically, we need to introduce provide some additional equations that tell us how the spacecraft's orbital plane changes during the execution of the manoeuvre.
So, whereas the first three equations of motion above provide information about the three orbital elements that specify the shape of the Keplerian orbital trajectory and the location of the spacecraft on that trajectory, an additional three equations of motion are needed to specify the orientation of the orbital plane of the trajectory on that orbital plane. Traditionally, this additional information is represented by three angle variables - the inclination, the longitude of the ascending node and the argument of periapsis.
In our case, we are also going to use three angle variables - but we are going to use a Roll-Pitch-Yaw convention in which we perform a 'roll' about the [MATH]r[/MATH]-axis of angle [MATH]\alpha(t)[/MATH], then perform a 'pitch' about the original [MATH]\theta[/MATH]-axis of angle [MATH]\beta(t)[/MATH] and, finally, perform a 'yaw' about the original [MATH]h[/MATH]-axis (i.e., normal to the trajectory plane) of angle [MATH]\gamma(t)[/MATH]. Given this, we end up with the following equations of motion for these angle variables:
[MATH]\alpha'(t)=-\frac{h(t)}{\mu\left(1+e_{r}(t)\right)}\,B_{h}\,\sec\beta(t)\,\cos\gamma(t)[/MATH][MATH]\beta'(t)=+\frac{h(t)}{\mu\left(1+e_{r}(t)\right)}\,B_{h}\,\sin\gamma(t)[/MATH][MATH]\gamma'(t)=-\frac{\mu^{2}\left(1+e_{r}(t)\right){}^{2}}{h(t)^{3}}-\frac{h(t)}{\mu\left(1+e_{r}(t)\right)}\,B_{h}\,\tan\beta(t)\,\cos\gamma(t)[/MATH]
where, in this case, [MATH]B_h[/MATH] is the off-plane component of applied thrust (or, more generally, gravitational perturbation). These three equations replace the single equation from the preceding two-dimensional case:
[MATH]\theta'(t)=\frac{\mu^{2}\left(1+e_{r}(t)\right){}^{2}}{h(t)^{3}}[/MATH]
although we can see an echo of that equation in the equation for [MATH]\gamma'(t)[/MATH]. In total, then to full describe the behaviour of the trajectory of the spacecraft during a trajectory, we must solve the following six first-order differential equations:
[MATH]h'(t) =\frac{h(t)^{2}}{\mu\left(1+e_{r}(t)\right)}\,B_{\theta}[/MATH][MATH]e_{r}'(t)-\theta'(t)\,e_{\theta}(t) =2\,\frac{h(t)}{\mu}\,B_{\theta}[/MATH][MATH]e_{\theta}'(t)+\theta'(t)\,e_{r}(t) =-\frac{h(t)}{\mu}\,B_{r}+\frac{h(t)}{\mu}\,\frac{e_{\theta}(t)}{1+e_{r}(t)}\,B_{\theta}[/MATH][MATH]\alpha'(t)=-\frac{h(t)}{\mu\left(1+e_{r}(t)\right)}\,B_{h}\,\sec\beta(t)\,\cos\gamma(t)[/MATH][MATH]\beta'(t)=+\frac{h(t)}{\mu\left(1+e_{r}(t)\right)}\,B_{h}\,\sin\gamma(t)[/MATH][MATH]\gamma'(t)=-\frac{\mu^{2}\left(1+e_{r}(t)\right){}^{2}}{h(t)^{3}}-\frac{h(t)}{\mu\left(1+e_{r}(t)\right)}\,B_{h}\,\tan\beta(t)\,\cos\gamma(t)[/MATH]
such that:
[MATH] \frac{\nu_e\,\gamma}{m_0-\gamma\,t}=\sqrt{B_r^2+B_\theta^2+B_h^2}[/MATH]
Messy, and complicated - but true.
If we manage to solve this set of equations, then from the solution for [MATH]\alpha(t)[/MATH], [MATH]\beta(t)[/MATH] and [MATH]\gamma(t)[/MATH], we can construct a rotation matrix [MATH]\boldsymbol{U}(t)[/MATH] from the following matrix elements:
[MATH]U_{1,1}=\cos\beta(t)\,\cos\gamma(t)[/MATH][MATH]U_{1,2}=\sin\alpha(t)\,\sin\beta(t)\,\cos\gamma(t)-\cos\alpha(t)\,\sin\gamma(t)[/MATH][MATH]U_{1,3}=\cos\alpha(t)\,\sin\beta(t)\,\cos\gamma(t)+\sin\alpha(t)\,\sin\gamma(t)[/MATH][MATH]U_{2,1}=\cos\beta(t)\,\sin\gamma(t)[/MATH][MATH]U_{2,2}=\sin\alpha(t)\,\sin\beta(t)\,\sin\gamma(t)+\cos\alpha(t)\,\cos\gamma(t)[/MATH][MATH]U_{2,3}=\cos\alpha(t)\,\sin\beta(t)\,\sin\gamma(t)-\sin\alpha(t)\,\cos\gamma(t)[/MATH][MATH]U_{3,1}=-\sin\beta(t)[/MATH][MATH]U_{3,2}=\sin\alpha(t)\,\cos\beta(t)[/MATH][MATH]U_{3,3}=\cos\alpha(t)\,\cos\beta(t)[/MATH]
Having constructed this rotation matrix, we can calculate how the coordinates of the reference frame change with time:
[MATH]\mathbf{\hat{\boldsymbol{r}}}(t)=\boldsymbol{U}.\mathbf{\hat{\boldsymbol{r}}}_{0}[/MATH][MATH]\mathbf{\hat{\boldsymbol{\theta}}}(t)=\boldsymbol{U}.\mathbf{\hat{\boldsymbol{\theta}}}_{0}[/MATH][MATH]\mathbf{\hat{\boldsymbol{h}}}(t)=\boldsymbol{U}.\mathbf{\hat{\boldsymbol{h}}}_{0}[/MATH]
And from this, we can calculate the planet-centric state vectors of the spacecraft as:
[MATH]\boldsymbol{r}(t)=\frac{h(t)^{2}}{\mu\left(1+e_{r}(t)\right)}\,\mathbf{\hat{\boldsymbol{r}}}(t)[/MATH][MATH]\boldsymbol{v}(t)=-\frac{\mu\,e_{\theta}(t)}{h(t)}\,\mathbf{\hat{\boldsymbol{r}}}(t)+\frac{\mu\,\left(1+e_{r}(t)\right)}{h(t)}\,\hat{\boldsymbol{\theta}}(t)[/MATH]
Alternatively, we can calculate the orbital elements as:
The semi-major-axis, [MATH]a[/MATH]
[MATH]a(t)=\frac{h(t)^{2}}{\mu\,(1-e_r(t)^2-e_\theta(t)^2)}[/MATH]
The orbital eccentricity, [MATH]e[/MATH]
[MATH]e(t)=\sqrt{e_r(t)^2+e_\theta(t)^2}[/MATH]
The true anomaly, [MATH]\nu[/MATH]
[MATH]\cos\nu(t)=\frac{e_{r}(t)}{\sqrt{e_{r}(t)^{2}+e_{\theta}(t)^{2}}}[/MATH]
(such that such if [MATH]\frac{e_{\theta}(t)}{e_{r}(t)+1}>0[/MATH] then [MATH]\nu\to2\,\pi-\nu[/MATH]).
To calculate the other three orbital elements, the inclination, the longitude of the ascending node and the argument of periapsis, we need to calculate the components of the angular momentum vector [MATH]\boldsymbol{h}(t)[/MATH] and the eccentricity vector [MATH]\boldsymbol{e}(t)[/MATH] in terms of a reference x-y-z coordinate system of a reference frame - usually either the ecliptic or equatorial reference frames. The x-y-z components of these vectors are calculated as:
[MATH]\hat{h}_x(t) = \hat{\boldsymbol{x}}.\boldsymbol{U}.\hat{\boldsymbol{h}}(t)[/MATH][MATH]\hat{h}_y(t) = \hat{\boldsymbol{y}}.\boldsymbol{U}.\hat{\boldsymbol{h}}(t)[/MATH][MATH]\hat{h}_z(t) = \hat{\boldsymbol{z}}.\boldsymbol{U}.\hat{\boldsymbol{h}}(t)[/MATH]
and
[MATH]e_x(t) = e_r(t)\,\hat{\boldsymbol{r}}(t).\hat{\boldsymbol{x}}+e_\theta(t)\,\hat{\boldsymbol{\theta}}(t).\hat{\boldsymbol{x}}[/MATH][MATH]e_y(t) = e_r(t)\,\hat{\boldsymbol{r}}(t).\hat{\boldsymbol{y}}+e_\theta(t)\,\hat{\boldsymbol{\theta}}(t).\hat{\boldsymbol{y}}[/MATH][MATH]e_z(t) = e_r(t)\,\hat{\boldsymbol{r}}(t).\hat{\boldsymbol{z}}+e_\theta(t)\,\hat{\boldsymbol{\theta}}(t).\hat{\boldsymbol{z}}[/MATH]
So, now we can calculate the orbital elements in the usual way:
The inclination, [MATH]\iota[/MATH]
[MATH]\cos\iota(t) = h_z(t)[/MATH]
The longitude of the ascending node, [MATH]\Omega[/MATH]
[MATH]\cos\Omega(t) = -\frac{\hat{h}_y(t)}{\sqrt{\hat{h}_x(t)^{2}+\hat{h}_y(t)^2}}[/MATH]
such that if [MATH]\hat{h}_x(t)<0[/MATH] then [MATH]\Omega\to2\,\pi-\Omega[/MATH].
The argument of periapsis, [MATH]\omega[/MATH]
[MATH]\cos\omega=\frac{\hat{h}_{x}(t)\,e_{y}(t)-\hat{h}_{y}(t)\,e_{x}(t)}{e(t)}[/MATH]
such that if [MATH]e_z(t)<0[/MATH] then [MATH]\omega\to 2\,\pi-\omega[/MATH].
A perturbed Keplerian gravitational field
So far, we have not worried terribly much about gravitational perturbations arising from, say, the Moon or other 'nearby' gravitating bodies. How do we include these perturbations? The formalism is very simple: we simply replace the thrust vector [MATH]\boldsymbol{B}[/MATH] with one that includes the gravitational perturbation, [MATH]\Delta\boldsymbol{P}[/MATH] such that
[MATH]B_r \to B_r + \Delta P_r[/MATH][MATH]B_\theta \to B_\theta + \Delta P_\theta[/MATH][MATH]B_h \to B_h + \Delta P_h[/MATH]
and then we proceed as before while ensuring that:[MATH] \frac{\nu_e\,\gamma}{m_0-\gamma\,t}=\sqrt{B_r^2+B_\theta^2+B_h^2}[/MATH]
Last edited:
