A generalised ideal rocket equation
- Thread starter Keithth G
- Start date
Enjo
Mostly harmless
- Joined
- Nov 25, 2007
- Messages
- 1,665
- Reaction score
- 13
- Points
- 38
- Location
- Germany
- Website
- www.enderspace.de
- Preferred Pronouns
- Can't you smell my T levels?
Great to hear. Of course the rocket equation soon to come is not the only reason why I'm happy 
RGClark
Mathematician
- Joined
- Jan 27, 2010
- Messages
- 1,635
- Reaction score
- 1
- Points
- 36
- Location
- Philadelphia
- Website
- exoscientist.blogspot.com
If using Schilling's calculator you have to be aware of some quirks.
First, for the option "Restartable Upper Stage?" select No, rather than the default Yes, otherwise the payload will be reduced. I think this is because the calculator keeps some propellant on reserve for the restart when you select Yes. Notably, this happens even if the rocket happens to be a SSTO, which wouldn't have an upper stage. That is, the payload is reduced when the Yes option is checked even in this case.
Also, for the Isp and thrust values always input the vacuum values, even for the first stage and/or boosters. This is because the calculator takes into account the diminution in thrust and Isp at sea level.
Another aspect is you should input for the inclination angle of the orbit the same value as the latitude of the launch site. So for Cape Canaveral it should be 28.5 in degrees. This is because having a different launch angle from the launch site latitude is equivalent to having to make an orbital angle change, which is a high delta-v cost.
BTW, so as not to derail KeithG's thread I'm moving this discussion to a separate thread.
Bob Clark
Last edited:
As promised, this note outlines an algorithm for finding the time when a retrograde burn must be started (and stopped) so as to insert into spacecraft into a circular orbit around a central gravitating body such as a planet or moon assuming constant engine thrust during the manoeuvre and while the thrust is directed along the retrograde direction. In addition, the algorithm permits ancillary quantities including, *inter alia*, the amount of fuel used, the Delta-V required to execute the burn, and the radius of the final circular orbit around the gravitating body.
The basics conceptual model upon which the algorithm is built
Although the mathematical framework described herein is general, and is applicable to craft with widely varying thrust capabilities, the basic scheme imagined here is that of a spacecraft approaching a gravitating body (planet or moon) in a hyperbolic orbit - or, possibly, a high eccentricity elliptical orbit - executing a single, retrograde burn to insert it into a circular orbit around the central gravitating body. To keep things simple, we make the assumption that the ballistic trajectory of the incoming spacecraft is Keplerian, which is to say that non-spherical gravity sources and other gravitational perturbations are not considered.
Outline of the algorithm
At the core of the algorithm is a mathematical function that will be described fully in the next section.
The core function has two inputs. Together, albeit indirectly, these inputs allow one to calculate the start and stop times for the retrograde burn given the thrust characteristics of the spacecraft. The two inputs are:
1. [MATH]R[/MATH], the orbital radius of the spacecraft at the start of the orbit insertion burn; and
2. [MATH]\delta t[/MATH], the duration of the circular orbit insertion burn once it has started.
Given these two inputs, the core function returns two numbers. These are the two (planar) components of the final orbital eccentricity vector, [MATH]e_r[/MATH] and [MATH]e_\theta[/MATH]. The orbital eccentricity vector is a vector that, from the centre of the gravitating body, points towards periapsis; and the length of which is the usual Keplerian orbital eccentricity, [MATH]e[/MATH]. The vector component [MATH]e_r[/MATH] is the magnitude of the eccentricity vector in the radial direction (i.e., from the centre of the planet to the position of the spacecraft); and the other vector component [MATH]e_\theta[/MATH] is the magnitude of the eccentricity perpendicular to the radial direction. Together, these two components yields information about the shape and orientation of the final orbit.
The conditions for achieving a circular orbit are to find input variables [MATH]R[/MATH] and [MATH]\delta t[/MATH] such that the orbital eccentricity output quantities are zero as the end of the orbit insertion burn. This search for input variables is a root-finding exercise in two variables and may be aided if a robust root-finding algorithm - e.g., perhaps, Brent's method - is available.
Having found the roots of the core function, one may then calculate the following quantities of interest:
1. [MATH]T^*[/MATH], the time when the orbital insertion burn should commence. [MATH]T^*[/MATH] may be earlier than or later than the current time. If it is earlier than the current time, then the opportunity to insert into a circular orbit is already passed;
2. [MATH]\delta t[/MATH], the duration of the burn;
2. [MATH]r_c[/MATH], the radius of the final circular orbit;
3. [MATH]\delta m[/MATH], the mass of the fuel consumed; and
4. [MATH]\Delta V[/MATH] the Delta-V required for orbital insertion (as calculated using the standard Ideal Rocket Equation)
The core numerical integration function in detail
As mentioned above, at the core function performs a numerical integration of a specific dynamical system. This numerical integration requires that a suitable numerical integrator, such as RK4, be available.
The specific dynamical system is as follows:
[MATH]h'(t) = -\frac{1}{\mu}\,\frac{\gamma\,v_{e}}{m_{0}-\gamma\,t}\,\frac{h(t)^{2}}{\sqrt{(1+e_{r}(t))^{2}+e_{\theta}(t)^{2}}}[/MATH]
[MATH]e_r'(t)-\theta'(t)\,e_{\theta}(t) = -\frac{2}{\mu}\,\frac{\gamma\,v_{e}}{m_{0}-\gamma\,t}\,\frac{1+e_{r}(t)}{\sqrt{(1+e_{r}(t))^{2}+e_{\theta}(t)^{2}}}\,h(t)[/MATH]
[MATH]e_\theta'(t)+\theta'(t)\,e_{r}(t) = -\frac{2}{\mu}\,\frac{\gamma\,v_{e}}{m_{0}-\gamma\,t}\,\frac{e_{\theta}(t)}{\sqrt{(1+e_{r}(t))^{2}+e_{\theta}(t)^{2}}}\,h(t)[/MATH]
[MATH] \theta'(t) =\frac{\mu ^2 \left(1+e_r(t)\right){}^2}{h(t)^3} [/MATH]
This system of equations is integrated from [MATH]0[/MATH] to [MATH]\delta t[/MATH] (one of the two input parameters) subject to the following initial conditions:
[MATH] h(0) = \sqrt{\mu\,a\,(1-e^{2})} [/MATH]
[MATH] e_r(0) = \frac{a\,(1-e^{2})}{R}-1 [/MATH]
[MATH] e_\theta(0) = \frac{1}{R}\,\sqrt{(1-e^{2})\,(a\,(1+e)-R)\,(R-a\,(1-e))} [/MATH]
[MATH] \theta(0) = 0 [/MATH]
Here, [MATH]R[/MATH] is the second of the two input parameters, the orbital radius of the spacecraft at the start of the orbit insertion burnThis dynamical system also contains a number of other parameters. These are:
1. [MATH]a[/MATH], the semi-major axis of the initial spacecraft's Keplerian trajectory prior to commencement of the burn;
2. [MATH]e[/MATH], the eccentricity of the initial spacecraft trajectory;
3. [MATH]\mu[/MATH], the gravitational parameter 'GM' for the central gravitating body;
4. [MATH]\gamma[/MATH], the fuel flow rate (kg/s) for the engines that will execute the retrograde manoeuvre;
5. [MATH]v_e[/MATH], the effective velocity of the exhaust gases from the engines; and
6. [MATH]m_0[/MATH], the total mass of the spacecraft at the start of the retrograde burn.
So, what do these four equations describe? Well, there are four dynamical quantities that are being tracked by integrating this dynamical system:
1. the specific angular momentum, [MATH]h(t)[/MATH] of the spacecraft;
2. the two components of the eccentricity vector, [MATH]e_r(t)[/MATH] and [MATH]e_\theta(t)[/MATH]; and
3. the change in angular position of the spacecraft in an inertial reference frame, [MATH]\theta(t)[/MATH], as it executes the orbit insertion burn.
Although these equations appear different from those outlined in the original post of this thread, they are formally equivalent. However, they make better use of invariants of a 2-body Keplerian system - namely the angular momentum and the orbital eccentricity.
Note that there are bounds on permissible values of [MATH]R[/MATH]. For a hyperbolic approach orbit, we require that [MATH]a\,(1-e) \le R[/MATH]$. And for a (high eccentricity) elliptic approach orbit, we require that [MATH]a\,(1-e) \le R \le a\,(a+e)[/MATH].
Calculating ancillary quantities
Once one has solved the problem for [MATH]R^*[/MATH] and [MATH]\delta t^*[/MATH], we can calculate other quantities that may be of interest:
a. The radius of the final circular orbit, [MATH]r_c[/MATH]:
[MATH]r_c = \frac{1}{\mu}\,h(\delta t^*)^2[/MATH]
b. The total mass of fuel consumed, [MATH]\delta m[/MATH]:
[MATH] \delta m=\gamma\,\delta t^*[/MATH]
c. The Delta-V required to execute the orbit insertion manoeuvre, [MATH]\Delta V[/MATH]:
[MATH] \Delta V = v_e\,\log\frac{m_0}{m_0-\delta m}[/MATH]
d. The time to the start of the orbit insertion manoeuvre, [MATH]\Delta T[/MATH]. First, calculate the true anomaly corresponding to both the spacecraft's current position and its position at the start of the orbit insertion burn:
[MATH] \nu_0 = -\tan^{-1}\frac{\sqrt{(1-e^{2})\,(a\,(1+e)-R_{0})\,(R_{0}-a\,(1-e))}}{a\,(1-e^{2})-R_{0}}[/MATH]and
[MATH] \nu^* = -\tan^{-1}\frac{\sqrt{(1-e^{2})\,(a\,(1+e)-R^{*})\,(R^{*}-a\,(1-e))}}{a\,(1-e^{2})-R^{*}}[/MATH]Next, calculate the mean anomaly at these two locations using Kepler's equation. And then, finally, convert the change in mean anomaly to a time difference using standard Keplerian conversions.
e. And for completeness, from the integrated solution to the equations of motion, one can calculate the orbital radius and radial velocity of the spacecraft by applying the following transformations:
[MATH] r(t) = \frac{h(t)^2}{\mu \left(1+e_r(t)\right)} [/MATH]
[MATH] r'(t) = -\frac{\mu \, e_{\theta }(t)}{h(t)}[/MATH]
A mathematical aside
It is of some interest to quickly examining the form of these equations in the absence of a thrust from the engines. In this case, the equations of motion reduce to the relatively simple form:
[MATH]h(t) = H [/MATH]
[MATH] e_r'(t)-\theta'(t)\,e_{\theta}(t) = 0[/MATH]
[MATH] e_{\theta}'(t)+\theta'(t)\,e_{r}(t) = 0 [/MATH]
[MATH] \theta'(t) = \frac{\mu^2}{H^3}\,\left(1+e_r(t)\right)^2 [/MATH]
or, if one wishes to change the 'time' variable in the problem from clock time to the true anomalay, [MATH]\nu[/MATH], we get:
[MATH]h(\nu) = H [/MATH]
[MATH] \frac{\partial^2}{\partial\nu^2}e_r(\nu)+e_r(\nu) = 0[/MATH]
[MATH] t'(\nu) = \frac{H^3}{\mu^2}\,\left(1+e_r(t)\right)^{-2}[/MATH]
These equations describe a system for which:
1. the specific angular momentum is conserved;
2. the eccentricity vector is of constant length and executes a rotation around completing one revolution in one orbital period.
3. The last equation is a disguised form of Kepler's Equation enabling a conversion from the true anomaly, [MATH]\nu[/MATH], back to ordinary clock time.
In short, although the form of the equations may be unfamiliar, this system of equations describes a standard Keplerian system. Consequently, the terms involving thrust ([MATH]\gamma \,v_e[/MATH]) on the right-hand side in the first three equations of the full system above describe the impact of a retrograde burn on the evolution of the shape and orientation of the conic trajectory of the craft's orbit.
Next steps
Given that this note is now overly long, in a subsequent note I will go through the details of applying this algorithm to a specific, low thrust orbit insertion example.
The basics conceptual model upon which the algorithm is built
Although the mathematical framework described herein is general, and is applicable to craft with widely varying thrust capabilities, the basic scheme imagined here is that of a spacecraft approaching a gravitating body (planet or moon) in a hyperbolic orbit - or, possibly, a high eccentricity elliptical orbit - executing a single, retrograde burn to insert it into a circular orbit around the central gravitating body. To keep things simple, we make the assumption that the ballistic trajectory of the incoming spacecraft is Keplerian, which is to say that non-spherical gravity sources and other gravitational perturbations are not considered.
Outline of the algorithm
At the core of the algorithm is a mathematical function that will be described fully in the next section.
The core function has two inputs. Together, albeit indirectly, these inputs allow one to calculate the start and stop times for the retrograde burn given the thrust characteristics of the spacecraft. The two inputs are:
1. [MATH]R[/MATH], the orbital radius of the spacecraft at the start of the orbit insertion burn; and
2. [MATH]\delta t[/MATH], the duration of the circular orbit insertion burn once it has started.
Given these two inputs, the core function returns two numbers. These are the two (planar) components of the final orbital eccentricity vector, [MATH]e_r[/MATH] and [MATH]e_\theta[/MATH]. The orbital eccentricity vector is a vector that, from the centre of the gravitating body, points towards periapsis; and the length of which is the usual Keplerian orbital eccentricity, [MATH]e[/MATH]. The vector component [MATH]e_r[/MATH] is the magnitude of the eccentricity vector in the radial direction (i.e., from the centre of the planet to the position of the spacecraft); and the other vector component [MATH]e_\theta[/MATH] is the magnitude of the eccentricity perpendicular to the radial direction. Together, these two components yields information about the shape and orientation of the final orbit.
The conditions for achieving a circular orbit are to find input variables [MATH]R[/MATH] and [MATH]\delta t[/MATH] such that the orbital eccentricity output quantities are zero as the end of the orbit insertion burn. This search for input variables is a root-finding exercise in two variables and may be aided if a robust root-finding algorithm - e.g., perhaps, Brent's method - is available.
Having found the roots of the core function, one may then calculate the following quantities of interest:
1. [MATH]T^*[/MATH], the time when the orbital insertion burn should commence. [MATH]T^*[/MATH] may be earlier than or later than the current time. If it is earlier than the current time, then the opportunity to insert into a circular orbit is already passed;
2. [MATH]\delta t[/MATH], the duration of the burn;
2. [MATH]r_c[/MATH], the radius of the final circular orbit;
3. [MATH]\delta m[/MATH], the mass of the fuel consumed; and
4. [MATH]\Delta V[/MATH] the Delta-V required for orbital insertion (as calculated using the standard Ideal Rocket Equation)
The core numerical integration function in detail
As mentioned above, at the core function performs a numerical integration of a specific dynamical system. This numerical integration requires that a suitable numerical integrator, such as RK4, be available.
The specific dynamical system is as follows:
[MATH]h'(t) = -\frac{1}{\mu}\,\frac{\gamma\,v_{e}}{m_{0}-\gamma\,t}\,\frac{h(t)^{2}}{\sqrt{(1+e_{r}(t))^{2}+e_{\theta}(t)^{2}}}[/MATH]
[MATH]e_r'(t)-\theta'(t)\,e_{\theta}(t) = -\frac{2}{\mu}\,\frac{\gamma\,v_{e}}{m_{0}-\gamma\,t}\,\frac{1+e_{r}(t)}{\sqrt{(1+e_{r}(t))^{2}+e_{\theta}(t)^{2}}}\,h(t)[/MATH]
[MATH]e_\theta'(t)+\theta'(t)\,e_{r}(t) = -\frac{2}{\mu}\,\frac{\gamma\,v_{e}}{m_{0}-\gamma\,t}\,\frac{e_{\theta}(t)}{\sqrt{(1+e_{r}(t))^{2}+e_{\theta}(t)^{2}}}\,h(t)[/MATH]
[MATH] \theta'(t) =\frac{\mu ^2 \left(1+e_r(t)\right){}^2}{h(t)^3} [/MATH]
This system of equations is integrated from [MATH]0[/MATH] to [MATH]\delta t[/MATH] (one of the two input parameters) subject to the following initial conditions:
[MATH] h(0) = \sqrt{\mu\,a\,(1-e^{2})} [/MATH]
[MATH] e_r(0) = \frac{a\,(1-e^{2})}{R}-1 [/MATH]
[MATH] e_\theta(0) = \frac{1}{R}\,\sqrt{(1-e^{2})\,(a\,(1+e)-R)\,(R-a\,(1-e))} [/MATH]
[MATH] \theta(0) = 0 [/MATH]
Here, [MATH]R[/MATH] is the second of the two input parameters, the orbital radius of the spacecraft at the start of the orbit insertion burnThis dynamical system also contains a number of other parameters. These are:
1. [MATH]a[/MATH], the semi-major axis of the initial spacecraft's Keplerian trajectory prior to commencement of the burn;
2. [MATH]e[/MATH], the eccentricity of the initial spacecraft trajectory;
3. [MATH]\mu[/MATH], the gravitational parameter 'GM' for the central gravitating body;
4. [MATH]\gamma[/MATH], the fuel flow rate (kg/s) for the engines that will execute the retrograde manoeuvre;
5. [MATH]v_e[/MATH], the effective velocity of the exhaust gases from the engines; and
6. [MATH]m_0[/MATH], the total mass of the spacecraft at the start of the retrograde burn.
So, what do these four equations describe? Well, there are four dynamical quantities that are being tracked by integrating this dynamical system:
1. the specific angular momentum, [MATH]h(t)[/MATH] of the spacecraft;
2. the two components of the eccentricity vector, [MATH]e_r(t)[/MATH] and [MATH]e_\theta(t)[/MATH]; and
3. the change in angular position of the spacecraft in an inertial reference frame, [MATH]\theta(t)[/MATH], as it executes the orbit insertion burn.
Although these equations appear different from those outlined in the original post of this thread, they are formally equivalent. However, they make better use of invariants of a 2-body Keplerian system - namely the angular momentum and the orbital eccentricity.
Note that there are bounds on permissible values of [MATH]R[/MATH]. For a hyperbolic approach orbit, we require that [MATH]a\,(1-e) \le R[/MATH]$. And for a (high eccentricity) elliptic approach orbit, we require that [MATH]a\,(1-e) \le R \le a\,(a+e)[/MATH].
Calculating ancillary quantities
Once one has solved the problem for [MATH]R^*[/MATH] and [MATH]\delta t^*[/MATH], we can calculate other quantities that may be of interest:
a. The radius of the final circular orbit, [MATH]r_c[/MATH]:
[MATH]r_c = \frac{1}{\mu}\,h(\delta t^*)^2[/MATH]
b. The total mass of fuel consumed, [MATH]\delta m[/MATH]:
[MATH] \delta m=\gamma\,\delta t^*[/MATH]
c. The Delta-V required to execute the orbit insertion manoeuvre, [MATH]\Delta V[/MATH]:
[MATH] \Delta V = v_e\,\log\frac{m_0}{m_0-\delta m}[/MATH]
d. The time to the start of the orbit insertion manoeuvre, [MATH]\Delta T[/MATH]. First, calculate the true anomaly corresponding to both the spacecraft's current position and its position at the start of the orbit insertion burn:
[MATH] \nu_0 = -\tan^{-1}\frac{\sqrt{(1-e^{2})\,(a\,(1+e)-R_{0})\,(R_{0}-a\,(1-e))}}{a\,(1-e^{2})-R_{0}}[/MATH]and
[MATH] \nu^* = -\tan^{-1}\frac{\sqrt{(1-e^{2})\,(a\,(1+e)-R^{*})\,(R^{*}-a\,(1-e))}}{a\,(1-e^{2})-R^{*}}[/MATH]Next, calculate the mean anomaly at these two locations using Kepler's equation. And then, finally, convert the change in mean anomaly to a time difference using standard Keplerian conversions.
e. And for completeness, from the integrated solution to the equations of motion, one can calculate the orbital radius and radial velocity of the spacecraft by applying the following transformations:
[MATH] r(t) = \frac{h(t)^2}{\mu \left(1+e_r(t)\right)} [/MATH]
[MATH] r'(t) = -\frac{\mu \, e_{\theta }(t)}{h(t)}[/MATH]
A mathematical aside
It is of some interest to quickly examining the form of these equations in the absence of a thrust from the engines. In this case, the equations of motion reduce to the relatively simple form:
[MATH]h(t) = H [/MATH]
[MATH] e_r'(t)-\theta'(t)\,e_{\theta}(t) = 0[/MATH]
[MATH] e_{\theta}'(t)+\theta'(t)\,e_{r}(t) = 0 [/MATH]
[MATH] \theta'(t) = \frac{\mu^2}{H^3}\,\left(1+e_r(t)\right)^2 [/MATH]
or, if one wishes to change the 'time' variable in the problem from clock time to the true anomalay, [MATH]\nu[/MATH], we get:
[MATH]h(\nu) = H [/MATH]
[MATH] \frac{\partial^2}{\partial\nu^2}e_r(\nu)+e_r(\nu) = 0[/MATH]
[MATH] t'(\nu) = \frac{H^3}{\mu^2}\,\left(1+e_r(t)\right)^{-2}[/MATH]
These equations describe a system for which:
1. the specific angular momentum is conserved;
2. the eccentricity vector is of constant length and executes a rotation around completing one revolution in one orbital period.
3. The last equation is a disguised form of Kepler's Equation enabling a conversion from the true anomaly, [MATH]\nu[/MATH], back to ordinary clock time.
In short, although the form of the equations may be unfamiliar, this system of equations describes a standard Keplerian system. Consequently, the terms involving thrust ([MATH]\gamma \,v_e[/MATH]) on the right-hand side in the first three equations of the full system above describe the impact of a retrograde burn on the evolution of the shape and orientation of the conic trajectory of the craft's orbit.
Next steps
Given that this note is now overly long, in a subsequent note I will go through the details of applying this algorithm to a specific, low thrust orbit insertion example.
Last edited:
As a specific example of the above, consider an unloaded - but fully-fuelled - Shuttle A approaching the Moon. The semi-major axis and eccentricity of its approach orbit are:
[MATH] a = 84324191.4 m [/MATH]
[MATH] e = 0.9599018 [/MATH]
This is a high eccentricity elliptic orbit with periapsis altitude of 1643 km. Its current orbital radius is 5000 km.
Now, imagine that (for some reason) only the Shuttle A's RCS thrusters are serviceable. Can the Shuttle A enter into a circular orbit? If so, what will the altitude of its final orbit be? How much fuel will it burn? And when should it commence its orbit insertion burn?
Well, as one might expect the preceding theory (previous post) can be used to answer the question. Let's determine a few basic parameters:
[MATH] \mu = 4.9048695\times10^{12} [/MATH] (for the Moon)
[MATH] v_e = 30,000[/MATH] m/s
[MATH] m_0 = 82,700 [/MATH] kg
[MATH] \gamma = 0.3 [/MATH] kg/s
If one now cranks the handle of the algorithm, one should be able to determine that:
So, yes, orbit insertion is possible (just). Orbit insertion will require 2,700 kg of fuel, it will take 9000 seconds; and the total Delta-V requirement for the manoeuvre is 943.6 m/s. At the end of the burn, the Shuttle A's RCS reserves will have been exhausted.
Now, obviously, I've back-solved the system so as to give a 20 km final orbital altitude and a 9000 second burn. But, equally, it is possible to use the above algorithm to confirm that these numbers are indeed correct.
[MATH] a = 84324191.4 m [/MATH]
[MATH] e = 0.9599018 [/MATH]
This is a high eccentricity elliptic orbit with periapsis altitude of 1643 km. Its current orbital radius is 5000 km.
Now, imagine that (for some reason) only the Shuttle A's RCS thrusters are serviceable. Can the Shuttle A enter into a circular orbit? If so, what will the altitude of its final orbit be? How much fuel will it burn? And when should it commence its orbit insertion burn?
Well, as one might expect the preceding theory (previous post) can be used to answer the question. Let's determine a few basic parameters:
[MATH] \mu = 4.9048695\times10^{12} [/MATH] (for the Moon)
[MATH] v_e = 30,000[/MATH] m/s
[MATH] m_0 = 82,700 [/MATH] kg
[MATH] \gamma = 0.3 [/MATH] kg/s
If one now cranks the handle of the algorithm, one should be able to determine that:
- the Shuttle A should commence its RCS orbit insertion burn when its orbital radius is 4384.3 km;
- the altitude of the final orbit is 20 km (a reduction of orbital periapsis altitude of 1623 km during the orbit insertion burn);
- the RCS thrusters need to burn for 9000 seconds (yes, 2.5 hours!); and
- in executing the orbit insertion manoeuvre, the Shuttle A will rotate 260 degrees around the Moon.
So, yes, orbit insertion is possible (just). Orbit insertion will require 2,700 kg of fuel, it will take 9000 seconds; and the total Delta-V requirement for the manoeuvre is 943.6 m/s. At the end of the burn, the Shuttle A's RCS reserves will have been exhausted.
Now, obviously, I've back-solved the system so as to give a 20 km final orbital altitude and a 9000 second burn. But, equally, it is possible to use the above algorithm to confirm that these numbers are indeed correct.
Last edited:
And for anyone interested, here (below) is a simple, stand-alone C program that performs the integration of the equations of motion for the example given above.
This implementation uses a Taylor series method to carryout the integration. Given initial conditions for the dynamical variables, the Taylor series method calculates high order derivatives (with respect to time) of those variables - in this case, up to 20th order. Those derivatives are calculated using a process called 'automatic differentiation' which is a fast, iterative technique for calculating derivatives of ODEs up to arbitrary order.
Taylor series methods have become more popular in recent years, in large measure because they are competitive with other, more familiar integration techniques but they also allow for high accuracy and considerable control over integration step-size and integration order.
In the code snippet below, I have made a quick estimate of optimal order and integration time-step for the system (given the input parameters); The optimal order is a little under 20 (but I rounded up to 20 to be on the safe side); and the optimal integration step-size is around 1000 seconds. This yields an accuracy that is close to machine number precision for this specific example.
So, to integrate the 9000 seconds (2.5 hours of the orbit insertion burn) requires just 9 integration steps. And, yes, the integration confirms that the circular orbit solution given in an earlier post in this thread is indeed one that leads to circular orbit insertion at an orbital altitude of 20 km.
Feel free to use and experiment with.
And for anyone running this program, as a check, you should find that the output of the program is:
This implementation uses a Taylor series method to carryout the integration. Given initial conditions for the dynamical variables, the Taylor series method calculates high order derivatives (with respect to time) of those variables - in this case, up to 20th order. Those derivatives are calculated using a process called 'automatic differentiation' which is a fast, iterative technique for calculating derivatives of ODEs up to arbitrary order.
Taylor series methods have become more popular in recent years, in large measure because they are competitive with other, more familiar integration techniques but they also allow for high accuracy and considerable control over integration step-size and integration order.
In the code snippet below, I have made a quick estimate of optimal order and integration time-step for the system (given the input parameters); The optimal order is a little under 20 (but I rounded up to 20 to be on the safe side); and the optimal integration step-size is around 1000 seconds. This yields an accuracy that is close to machine number precision for this specific example.
So, to integrate the 9000 seconds (2.5 hours of the orbit insertion burn) requires just 9 integration steps. And, yes, the integration confirms that the circular orbit solution given in an earlier post in this thread is indeed one that leads to circular orbit insertion at an orbital altitude of 20 km.
Feel free to use and experiment with.
Code:
//
// main.c
// IdealRocket
//
// A procedure to use 'automatic differentiation' to allow fast computation of the derivtaives
// to arbitrarily high order - first used in celestial mechanics problems
// For further details see, for example, http://www.maia.ub.edu/~angel/taylor/taylor.pdf
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int taylorIntegation( double er[21], // time derivatives of the radial component of the orbital eccentricity vector
double et[21], // time derivatives of the transbsrese component of the orbital eccentricity vector
double th[21], // time derivatives of the angular component of the spacecraft's position
double h [21], // time derivatives of the spacecrafts specific angular momentum
double f [21], // time derivatives of the spacecraft's acceleration
double ve, // the effective exhaust velocity (m/s) of the spacecraft's engines
double mu // the 'GM' parameter for the central gravitating body
) {
// declare the working arrays
double u01[20], u02[20], u05[20], u06[20], u07[20], u08[20], u09[20], u10[20], u11[20], u12[20], u13[20], u14[20], u15[20];
double v01[20], v02[20];
double w01[20], w02[20], w04[20], w05[20], w06[20], w07[20], w08[20], w09[20], w10[20];
// declare miscellaeous
double alpha;
// the integration step - building Taylor series expansions up to 20th order.
for ( int n = 0; n < 20; n++) {
u01[n] = er[n];
u02[n] = et[n];
v01[n] = f [n];
w01[n] = h [n];
// calculate f' - (the derivatives of the spacecraft's thrust function)
v02[n] = 0;
for (int j = 0; j <= n; j++) { // f(t)^2
v02[n] += v01[j] * v01[n-j];
}
v02[n] = v02[n] / ve; // f(t)^2 / ve
f[n+1] = v02[n] / (n+1); // f'(t) = f(t)^2 / ve
// calculate h' and higher order derivatives
w02[n] = 0;
for (int j = 0; j <= n; j++) { // h(t)^2
w02[n] += w01[j] * w01[n-j];
}
w04[n] = 0;
for (int j = 0; j <= n; j++) { // f(t) * h(t)^2
w04[n] += w02[j] * v01[n-j];
}
w05[n] = u01[n];
if (n == 0) { // 1 + e_r(t)
w05[n] += 1;
}
w06[n] = 0;
for (int j = 0; j <= n; j++) { // (1 + e_r(t))^2
w06[n] += w05[j] * w05[n-j];
}
w07[n] = 0;
for (int j = 0; j <= n; j++) { // e_t(t)^2
w07[n] += u02[j] * u02[n-j];
}
w08[n] = w06[n] + w07[n]; // A = (1 + e_r(t))^2 + e_t(t)^2
alpha = -0.5;
if (n == 0) {
w09[n] = pow( w08[0], alpha ); // A^-1/2
} else {
w09[n] = 0;
for (int j = 0; j <= n-1; j++) {
w09[n] += (n * alpha - j * (1 + alpha)) * w08[n-j] * w09[j];
}
w09[n] = w09[n] / n / w08[0];
}
w10[n] = 0;
for (int j = 0; j <= n; j++) {
w10[n] += w09[j] * w04[n-j];
}
w10[n] = - w10[n] / mu; // - mu^-1 * f(t) * h(t)^2 * A^-1/2
h[n+1] = w10[n] / (n+1);
// calculate theta' and higher order derivatives
alpha = -3.0;
if (n == 0) {
u05[n] = pow( w01[0], alpha ); // h(t)^-3
} else {
u05[n] = 0;
for (int j = 0; j <= n-1; j++) {
u05[n] += (n * alpha - j * (1 + alpha)) * w01[n-j] * u05[j];
}
u05[n] = u05[n] / n / w01[0];
}
u06[n] = 0;
for (int j = 0; j <= n; j++) {
u06[n] += u05[j] * w06[n-j];
}
u06[n] = u06[n] * mu * mu; // mu^2 * (1 + e_r(t))^2 * h(t)^-3
th[n+1] = u06[n] / (n+1); // th'(t) = u06 / (n+1)
// calculate e_r' and higher order derivatives
u07[n] = 0;
for (int j = 0; j <= n; j++) {
u07[n] += u06[j] * u02[n-j]; // mu^2 * (1 + e_r(t))^2 * h(t)^-3 * e_t(t)
}
u08[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t)^2 * A^-1/2 * (1 + e_r(t))
u08[n] += w10[j] * w05[n-j];
}
alpha = -1.0;
if (n == 0) {
u09[n] = pow( w01[0], alpha ); // h(t)^-1
} else {
u09[n] = 0;
for (int j = 0; j <= n-1; j++) {
u09[n] += (n * alpha - j * (1 + alpha)) * w01[n-j] * u09[j];
}
u09[n] = u09[n] / n / w01[0];
}
u10[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t) * A^-1/2 * (1 + e_r(t))
u10[n] += u09[j] * u08[n-j];
}
u11[n] = + u07[n] + 2 * u10[n]; // mu^2*(1+e_r(t))^2*h(t)^-3*e_t(t) - 2*mu^-1*f(t)*h(t)*A^-1/2*(1 + e_r(t))
er[n+1] = u11[n] / (n+1); // e_r'(t) = u11 / (n+1)
//calculate e_t' and higher order derivatives
u12[n] = 0;
for (int j = 0; j <= n; j++) { // mu^2 * (1 + e_r(t))^2 * h(t)^-3 * e_r(t)
u12[n] += u06[j] * u01[n-j];
}
u13[n] = 0;
for (int j = 0; j <= n; j++) {
u13[n] += w10[j] * u02[n-j]; // - mu^-1 * f(t) * h(t)^2 * A^-1/2 * e_t(t)
}
u14[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t) * A^-1/2 * e_t(t)
u14[n] += u13[j] * u09[n-j];
}
u15[n] = - u12[n] + 2 * u14[n]; // - mu^2*(1 + e_r(t))^2*h(t)^-3*e_r(t) - 2*mu^-1*f(t)*h(t)*A^-1/2*e_t(t)
et[n+1] = u15[n] / (n+1); // e_t'(t) = u15 / (n+1)
}
return 0;
}
int main(int argc, const char * argv[]) {
// declare constants - quantities that need to be passed to the Taylor integration routine
double mu = 4.9048695e12; // the standard gravitational parameter for the central body
double ve = 30000.0; // the effective velocity of engine exhaust gases
// devlare constants - quatities used to define the intitial conditons for the Taylor integration routine
double a = 84324191.410; // the semi-major axis of the spacecraft's current orbit
double e = 0.959901756305179; // the orbital eccentricity of the spacecraft's current orbit
double m0 = 87200.0; // the initial mass of the spacecraft
double g = 0.30000; // the fuel consumption rate
double r = 6.122297437930348e6; // the orbital radius at the start of the orbit insertion burn
// declare result arrays - normalised Taylor series derivatives up to and including 10th order
double er[21]; // the component of the eccentricity vector in the radial direction
double et[21]; // the component of the eccentricity vector in the 'theta' direction
double th[21]; // the spatial angular coordinate of the spacecraft
double h [21]; // the specific angular momentum of the spacecraft
double f [21]; // the spacraft's thrust function g * ve / (m0 - g * t)
// initalise the values to start the integration processteration process
er[0] = a * (1 - e * e) / r - 1;
et[0] = sqrt( (1 - e * e) * (a * e + a - r) * (a * e - a + r)) / r;
th[0] = 0.0;
h [0] = sqrt(mu * a * (1 - e * e));
f [0] = g * ve / m0;
// declare miscellanous quantities
double temp;
int status;
// set the time-step to 1000 seconds
double deltat = 1000.0;
// set the initial time at the start of the burn to 0.0 seconds
double tm = 0.0;
// carry out the integration step nine times - with each integration step spanning 1000 seconds
for (int count = 0; count < 9; count++) {
status = taylorIntegation(er, et, th, h, f, ve, mu);
// print out the value of the variables at the end of each 1000 second time-step using Horner's method
tm += deltat;
printf("tm = %2.2f\n", tm); // print the time after start of the burn (seconds)
temp = 0;
for (int j = 20; j >=0; j--) {
temp = (temp * deltat + er[j]);
}
printf("e_r = %2.10f\n", temp); // print the value of e_r - the radial component of the eccentricity vector.
er[0] = temp; // reset the initial value of the
temp = 0;
for (int j = 20; j >=0; j--) {
temp = (temp * deltat + et[j]);
}
printf("e_t = %2.10f\n", temp); // print the value of e_t - the other component of the eccenricity vector.
et[0] = temp; // reset the initial value of the transverse component of the eccentricity
temp = 0;
for (int j = 20; j >=0; j--) {
temp = (temp * deltat + th[j]);
}
printf("theta = %2.10f\n", temp); // print the angle of rotation since the start of the burn (radians)
th[0] = temp; // reset the initial value of the angle variable
temp = 0;
for (int j = 20; j >=0; j--) {
temp = (temp * deltat + f[j]);
}
printf("f = %2.10f\n", temp); // print the vehicle acceleration (m/s)
f[0] = temp; // reset the initial value of the spacecraft's acceleration
temp = 0;
for (int j = 20; j >=0; j--) {
temp = (temp * deltat + h[j]);
}
printf("h = %10.2f\n", temp); // print the specific angular momentum of the spacecraft
h[0] = temp; // reset the initial value of the specific angular momentum
printf("r(t) = %8.1f\n", h[0] * h[0] / mu / (1 + er[0])); // print the orbital radius (m)
printf("v_r(t)= %7.2f\n", - mu * et[0] / h[0]); // print the radial speed (m/s)
printf("v_t(t)= %7.2f\n", + mu * (1 + er[0]) / h[0]); // print the transverse speed (m/s)
printf("\n");
}
return 0;
}And for anyone running this program, as a check, you should find that the output of the program is:
Code:
tm = 1000.00
e_r = 0.0503950381
e_t = 0.7990964443
theta = 0.1674254229
f = 0.1035673188
h = 5244131904.98
r(t) = 5337859.0
v_r(t)= -747.40
v_t(t)= 982.44
tm = 2000.00
e_r = 0.0252932435
e_t = 0.6689526119
theta = 0.3712433940
f = 0.1039260970
h = 4822586804.35
r(t) = 4624710.5
v_r(t)= -680.37
v_t(t)= 1042.79
tm = 3000.00
e_r = 0.0113568631
e_t = 0.5631421578
theta = 0.6230878888
f = 0.1042873696
h = 4440068757.65
r(t) = 3974179.8
v_r(t)= -622.09
v_t(t)= 1117.23
tm = 4000.00
e_r = 0.0139696676
e_t = 0.4766587271
theta = 0.9408899815
f = 0.1046511628
h = 4098956739.84
r(t) = 3378269.2
v_r(t)= -570.38
v_t(t)= 1213.33
tm = 5000.00
e_r = 0.0389015112
e_t = 0.3999912296
theta = 1.3533248165
f = 0.1050175029
h = 3800193274.63
r(t) = 2834063.2
v_r(t)= -516.26
v_t(t)= 1340.90
tm = 6000.00
e_r = 0.0868518646
e_t = 0.3149779526
theta = 1.9047473335
f = 0.1053864169
h = 3542567746.37
r(t) = 2354173.8
v_r(t)= -436.10
v_t(t)= 1504.80
tm = 7000.00
e_r = 0.1334483968
e_t = 0.1971123852
theta = 2.6459230085
f = 0.1057579318
h = 3321015688.84
r(t) = 1983867.4
v_r(t)= -291.12
v_t(t)= 1674.01
tm = 8000.00
e_r = 0.1111519905
e_t = 0.0610401719
theta = 3.5677363946
f = 0.1061320755
h = 3124190301.40
r(t) = 1790911.1
v_r(t)= -95.83
v_t(t)= 1744.47
tm = 9000.00
e_r = 0.0000000004
e_t = 0.0000000279
theta = 4.5388623123
f = 0.1065088757
h = 2936453749.33
r(t) = 1758000.0
v_r(t)= -0.00
v_t(t)= 1670.34
Last edited:
Enjo
Mostly harmless
- Joined
- Nov 25, 2007
- Messages
- 1,665
- Reaction score
- 13
- Points
- 38
- Location
- Germany
- Website
- www.enderspace.de
- Preferred Pronouns
- Can't you smell my T levels?
Thanks. That's closer to what I need, but still, there are a lot of outputs. How to get them down just to the burn time required at a specific point in time?
Enjo
I'll expand on this a little later, but there are two controllable inputs and two outputs that need to be set to zero.
The controllable inputs in the routine are:
1. 'r', the orbital radius at which the orbit insertion burn commences; and
2. 'tm', the total burn time.
These two inputs need to be adjusted until 'e_r = 0' and 'e_t = 0' (the conditions needed to achieve a circular orbit). This is a root-finding problem in two-dimensions. If you don't have a routine to hand to solve that problem, I will write one for you - or, better yet, download one that someone else has posted.
All of the other outputs can be ignored until you have found the solution to the root-finding problem, The most relevant additional output is the radius of the final circular orbit. In general, this will be different from (and lower than) orbital periapsis prior to the start of the insertion burn.
---------- Post added at 01:31 PM ---------- Previous post was at 06:47 AM ----------
Enjo,
Some repackaged code which makes more explicit the nature of the map from input to outputs. When I find a bit more time, I'll insert this into a two-variable root-finding algorithm.
I'll expand on this a little later, but there are two controllable inputs and two outputs that need to be set to zero.
The controllable inputs in the routine are:
1. 'r', the orbital radius at which the orbit insertion burn commences; and
2. 'tm', the total burn time.
These two inputs need to be adjusted until 'e_r = 0' and 'e_t = 0' (the conditions needed to achieve a circular orbit). This is a root-finding problem in two-dimensions. If you don't have a routine to hand to solve that problem, I will write one for you - or, better yet, download one that someone else has posted.
All of the other outputs can be ignored until you have found the solution to the root-finding problem, The most relevant additional output is the radius of the final circular orbit. In general, this will be different from (and lower than) orbital periapsis prior to the start of the insertion burn.
---------- Post added at 01:31 PM ---------- Previous post was at 06:47 AM ----------
Enjo,
Some repackaged code which makes more explicit the nature of the map from input to outputs. When I find a bit more time, I'll insert this into a two-variable root-finding algorithm.
Code:
// main.c
// IdealRocket
//
// Created by Keith Gelling on 27/10/2016.
// Copyright © 2016 Keith Gelling. All rights reserved.
//
//
// A procedure to use 'automatic differentiation' to allow fast computation of the derivtaives
// to arbitrarily high order - first used in celestial mechanics problems
// For further details see, for example, http://www.maia.ub.edu/~angel/taylor/taylor.pdf
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// parameters used to define the initial orbit of the spacecraft
#define MU 4.9048695e12 // the 'GM' parameter for the central gravitating body (m^3 s^-2)
#define SMAJ 84324191.410 // the semi-major axis of the spacecraft's current orbit (m)
#define ECC 0.959901756305179 // the orbital eccentricity of the spacecraft's current orbit (n/a)
// parameters used to define the initial state of the spacecraft
#define M0 87200.0 // the spacecraft's initial mass (kg)
#define GAMMA 0.30000 // the fuel fow rate (kg/s)
#define VE 30000.0 // the effective velocity of engine exhaust gases (m/s)
// parameters used to define numerical integration order and step-size
#define DELTAT 1000.0 // the integration setp-size (s)
#define N 20 // the order of the Taylor series expansion
int taylorIntegation( double er[N+1], // time derivatives of the radial component of the orbital eccentricity vector
double et[N+1], // time derivatives of the transbsrese component of the orbital eccentricity vector
double th[N+1], // time derivatives of the angular component of the spacecraft's position
double h [N+1], // time derivatives of the spacecrafts specific angular momentum
double f [N+1] // time derivatives of the spacecraft's acceleration
) {
// declare the working arrays
double u01[N], u02[N], u05[N], u06[N], u07[N], u08[N], u09[N], u10[N], u11[N], u12[N], u13[N], u14[N], u15[N];
double v01[N], v02[N];
double w01[N], w02[N], w04[N], w05[N], w06[N], w07[N], w08[N], w09[N], w10[N];
// declare miscellaeous
double alpha;
// the integration step - building Taylor series expansions up to 20th order.
for ( int n = 0; n < N; n++) {
u01[n] = er[n];
u02[n] = et[n];
v01[n] = f [n];
w01[n] = h [n];
// calculate f'(t) - (the derivatives of the spacecraft's thrust function)
v02[n] = 0;
for (int j = 0; j <= n; j++) { // f(t)^2
v02[n] += v01[j] * v01[n-j];
}
v02[n] = v02[n] / VE; // f(t)^2 / ve
f[n+1] = v02[n] / (n+1); // f'(t) = f(t)^2 / ve
// calculate h'(t) and higher order derivatives
w02[n] = 0;
for (int j = 0; j <= n; j++) { // h(t)^2
w02[n] += w01[j] * w01[n-j];
}
w04[n] = 0;
for (int j = 0; j <= n; j++) { // f(t) * h(t)^2
w04[n] += w02[j] * v01[n-j];
}
w05[n] = u01[n];
if (n == 0) { // 1 + e_r(t)
w05[n] += 1;
}
w06[n] = 0;
for (int j = 0; j <= n; j++) { // (1 + e_r(t))^2
w06[n] += w05[j] * w05[n-j];
}
w07[n] = 0;
for (int j = 0; j <= n; j++) { // e_t(t)^2
w07[n] += u02[j] * u02[n-j];
}
w08[n] = w06[n] + w07[n]; // A = (1 + e_r(t))^2 + e_t(t)^2
alpha = -0.5;
if (n == 0) {
w09[n] = pow( w08[0], alpha ); // A^-1/2
} else {
w09[n] = 0;
for (int j = 0; j <= n-1; j++) {
w09[n] += (n * alpha - j * (1 + alpha)) * w08[n-j] * w09[j];
}
w09[n] = w09[n] / n / w08[0];
}
w10[n] = 0;
for (int j = 0; j <= n; j++) {
w10[n] += w09[j] * w04[n-j];
}
w10[n] = - w10[n] / MU; // - mu^-1 * f(t) * h(t)^2 * A^-1/2
h[n+1] = w10[n] / (n+1);
// calculate theta'(t) and higher order derivatives
alpha = -3.0;
if (n == 0) {
u05[n] = pow( w01[0], alpha ); // h(t)^-3
} else {
u05[n] = 0;
for (int j = 0; j <= n-1; j++) {
u05[n] += (n * alpha - j * (1 + alpha)) * w01[n-j] * u05[j];
}
u05[n] = u05[n] / n / w01[0];
}
u06[n] = 0;
for (int j = 0; j <= n; j++) {
u06[n] += u05[j] * w06[n-j];
}
u06[n] = u06[n] * MU * MU; // mu^2 * (1 + e_r(t))^2 * h(t)^-3
th[n+1] = u06[n] / (n+1); // th'(t) = u06 / (n+1)
// calculate e_r'(t) and higher order derivatives
u07[n] = 0;
for (int j = 0; j <= n; j++) {
u07[n] += u06[j] * u02[n-j]; // mu^2 * (1 + e_r(t))^2 * h(t)^-3 * e_t(t)
}
u08[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t)^2 * A^-1/2 * (1 + e_r(t))
u08[n] += w10[j] * w05[n-j];
}
alpha = -1.0;
if (n == 0) {
u09[n] = pow( w01[0], alpha ); // h(t)^-1
} else {
u09[n] = 0;
for (int j = 0; j <= n-1; j++) {
u09[n] += (n * alpha - j * (1 + alpha)) * w01[n-j] * u09[j];
}
u09[n] = u09[n] / n / w01[0];
}
u10[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t) * A^-1/2 * (1 + e_r(t))
u10[n] += u09[j] * u08[n-j];
}
u11[n] = + u07[n] + 2 * u10[n]; // mu^2*(1+e_r(t))^2*h(t)^-3*e_t(t) - 2*mu^-1*f(t)*h(t)*A^-1/2*(1 + e_r(t))
er[n+1] = u11[n] / (n+1); // e_r'(t) = u11 / (n+1)
//calculate e_t'(t) and higher order derivatives
u12[n] = 0;
for (int j = 0; j <= n; j++) { // mu^2 * (1 + e_r(t))^2 * h(t)^-3 * e_r(t)
u12[n] += u06[j] * u01[n-j];
}
u13[n] = 0;
for (int j = 0; j <= n; j++) {
u13[n] += w10[j] * u02[n-j]; // - mu^-1 * f(t) * h(t)^2 * A^-1/2 * e_t(t)
}
u14[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t) * A^-1/2 * e_t(t)
u14[n] += u13[j] * u09[n-j];
}
u15[n] = - u12[n] + 2 * u14[n]; // - mu^2*(1 + e_r(t))^2*h(t)^-3*e_r(t) - 2*mu^-1*f(t)*h(t)*A^-1/2*e_t(t)
et[n+1] = u15[n] / (n+1); // e_t'(t) = u15 / (n+1)
}
return 0;
}
int updateInitialValues( double er[N+1], // time derivatives of the radial component of the orbital eccentricity vector
double et[N+1], // time derivatives of the transbsrese component of the orbital eccentricity vector
double th[N+1], // time derivatives of the angular component of the spacecraft's position
double h [N+1], // time derivatives of the spacecrafts specific angular momentum
double f [N+1], // time derivatives of the spacecraft's acceleration
double dt // the time step
) {
double temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + er[j]);}; er[0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + et[j]);}; et[0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + th[j]);}; th[0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + f [j]);}; f [0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + h [j]);}; h [0] = temp;
return 0;
}
int f ( double R, // the orbital radius at the start of the orbit insertion burn
double tm, // the duration of the orbit insertion burn
double e[2] // the radial and transverse components of the eccentricity vector at the end of the orbital inserion burn
) {
// declare result arrays - normalised Taylor series derivatives up to and including 10th order
double er[N+1]; // the component of the eccentricity vector in the radial direction
double et[N+1]; // the component of the eccentricity vector in the 'theta' direction
double th[N+1]; // the spatial angular coordinate of the spacecraft
double h [N+1]; // the specific angular momentum of the spacecraft
double f [N+1]; // the spacraft's thrust function g * ve / (m0 - g * t)
// initalise the values to start the integration processteration process
er[0] = SMAJ * (1 - ECC * ECC) / R - 1;
et[0] = sqrt( (1 - ECC * ECC) * (SMAJ * ECC + SMAJ - R) * (SMAJ * ECC - SMAJ + R)) / R;
th[0] = 0.0;
h [0] = sqrt(MU * SMAJ * (1 - ECC * ECC));
f [0] = GAMMA * VE / M0;
// declare miscellanous quantities
int status;
double t;
t = tm;
// carry out the integration step nine times - with each integration step spanning 1000 seconds
while (t >= DELTAT) {
status = taylorIntegation (er, et, th, h, f);
status = updateInitialValues (er, et, th, h, f, DELTAT);
t -= DELTAT;
}
status = taylorIntegation (er, et, th, h, f);
status = updateInitialValues (er, et, th, h, f, t);
t = 0;
e[0] = er[0];
e[1] = et[0];
return 0;
}
int main(int argc, const char * argv[]) {
// INPUTS
double r = 6.122297437930348e6; // the orbital radius at the start of the orbit insertion burn
double tm = 9000.0; // the duration of the orbit insertion burn
// OUTPUTS
double e[2]; // the radial and transverse components of the eccentricity vector at the end of the orbital insertion burn
// e[0] holds the radial component; e[1] holds the transverse component.
// map inputs to outputs
int status;
status = f(r, tm, e);
// print the outputs
printf("e_r %2.10f\n", e[0]);
printf("e_t %2.10f\n", e[1]);
return 0;
}
Last edited:
Enjo
Here (below) is a revised version of the code that includes the two-variable root-finding solution. This code is valid for a wide range of spacecraft with widely varying thrust characteristics approaching periapsis on an elliptical or hyperbolic trajectory and wishing to enter a circular orbit by executing a retrograde burn. (I have not yet considered circularisation at apoapsis, but the code would need some minor changes to address that case.)
The calculation of the duration of the burn and the time (before periapsis) to execute the orbit insertion burn is iterative. This means that the calculation needs a reasonable initial guess in order to converge quickly. I propose that BTC's current estimates provide that initial guess. In particular, the two quantities of interest are:
1. the number of seconds [MATH]\delta t[/MATH] before periapsis when the orbit insertion burn starts; and
2. the duration, [MATH]t_m[/MATH] of the orbit insertion burn.
First, you need to convert [MATH]\delta t[/MATH] into an orbital radius by solving for the eccentric anomaly, if the approach orbit is elliptical:
[MATH]\delta t=\sqrt{\frac{a^3}{\mu }}\, (E-e\, \sin E)[/MATH]
and the following if the approach orbit is hyperbolic:
[MATH]\delta t=\sqrt{\frac{-a^3}{\mu }}\, (e\, \sinh E - E)[/MATH]
Second, you need to calculate the initial vale of [MATH]R[/MATH]:
[MATH]R = a (1- e\, \cos E)[/MATH]
if the approach orbit is elliptical; and
[MATH]R = a (1 - e\, \cosh E )[/MATH]
if the approach orbit is hyperbolic.
The resulting value of [MATH]R[/MATH] and [MATH]t_m[/MATH] are the initial values for the orbit insertion algorithm below.
In most cases, convergence should be rapid - but, as with all things involving non-linear root-finding, there are no guarantees. At the moment, the algorithm has few checks to test for convergence. Enhancing the 'robustness' of the algorithm may be worthwhile thing to do at some later stage.
The outputs of the algorithm are:
1. the time in seconds prior to periapsis when the orbit insertion burn should start. This is a 'refined' estimate of BTC's initial guess.
2. the duration of the orbit insertion burn in seconds. Again, this is refinement of BTCs estimate of the same.
3. the radius of the final circular orbit. This is an additional piece of information that BTC does not currently provide which, in effect, assumes that the radius of the final orbit will be the same as the orbital periapsis. For short-duration burn, this is almost true, but for low thrust, long duration orbit insertion burns there can be a significant drop in the final orbital radius. What you choose to do with this additional information is up to you, but I think it a value that should be reported by BTC.)
The algorithm also needs to link-up to relevant Orbiter parameters:
1. the value of the current semi-major axis, [MATH]a[/MATH];
2. the value of the current orbital eccentricity, [MATH]e[/MATH];
3. the value of 'GM' for the central gravitating body, [MATH]\mu[/MATH];
4. the current mass of the spacecraft, [MATH]m_0[/MATH];
5. the mass flow rate for the engines chosen to execute the retrograde burn, [MATH]\gamma[/MATH]; and
6. the effective exhaust velocity of the combustion gases of those engines [MATH]v_e[/MATH].
To obtain these values, you will of course need to interrogate Orbiter's system.
For the specific example given in the code below (which considers a Shuttle A on approach to perilune aiming to use its main engines to enter into a circular orbit. The algorithm reports that:
a. the duration of the orbit insertion burn is: 19.595 seconds;
b. the orbit insertion burn should commence 10.277 seconds before periapsis; and
c. the radius of the final circular orbit will be 3380913.4 metres.
In comparison, BTC reports that the duration of the orbit insertion burn is: 19.591 seconds; and the insertion burn should commence 9.831 seconds before periapsis. Here, the algorithm opts for a slightly longer burn time because the final orbital radius is very slightly lower than that assumed by BTC. The material difference between the two solutions is that the algorithm estimates that the burn should commence 0.45 seconds earlier than that calculated by BTC. For fast, high-thrust orbit insertion burns such as this the differences in solutions are not particularly material. But for low-thrust orbit insertion burns, the differences are much greater.
At this stage, to hook this up to BTC, is there anything else that you need?
Here (below) is a revised version of the code that includes the two-variable root-finding solution. This code is valid for a wide range of spacecraft with widely varying thrust characteristics approaching periapsis on an elliptical or hyperbolic trajectory and wishing to enter a circular orbit by executing a retrograde burn. (I have not yet considered circularisation at apoapsis, but the code would need some minor changes to address that case.)
The calculation of the duration of the burn and the time (before periapsis) to execute the orbit insertion burn is iterative. This means that the calculation needs a reasonable initial guess in order to converge quickly. I propose that BTC's current estimates provide that initial guess. In particular, the two quantities of interest are:
1. the number of seconds [MATH]\delta t[/MATH] before periapsis when the orbit insertion burn starts; and
2. the duration, [MATH]t_m[/MATH] of the orbit insertion burn.
First, you need to convert [MATH]\delta t[/MATH] into an orbital radius by solving for the eccentric anomaly, if the approach orbit is elliptical:
[MATH]\delta t=\sqrt{\frac{a^3}{\mu }}\, (E-e\, \sin E)[/MATH]
and the following if the approach orbit is hyperbolic:
[MATH]\delta t=\sqrt{\frac{-a^3}{\mu }}\, (e\, \sinh E - E)[/MATH]
Second, you need to calculate the initial vale of [MATH]R[/MATH]:
[MATH]R = a (1- e\, \cos E)[/MATH]
if the approach orbit is elliptical; and
[MATH]R = a (1 - e\, \cosh E )[/MATH]
if the approach orbit is hyperbolic.
The resulting value of [MATH]R[/MATH] and [MATH]t_m[/MATH] are the initial values for the orbit insertion algorithm below.
In most cases, convergence should be rapid - but, as with all things involving non-linear root-finding, there are no guarantees. At the moment, the algorithm has few checks to test for convergence. Enhancing the 'robustness' of the algorithm may be worthwhile thing to do at some later stage.
The outputs of the algorithm are:
1. the time in seconds prior to periapsis when the orbit insertion burn should start. This is a 'refined' estimate of BTC's initial guess.
2. the duration of the orbit insertion burn in seconds. Again, this is refinement of BTCs estimate of the same.
3. the radius of the final circular orbit. This is an additional piece of information that BTC does not currently provide which, in effect, assumes that the radius of the final orbit will be the same as the orbital periapsis. For short-duration burn, this is almost true, but for low thrust, long duration orbit insertion burns there can be a significant drop in the final orbital radius. What you choose to do with this additional information is up to you, but I think it a value that should be reported by BTC.)
The algorithm also needs to link-up to relevant Orbiter parameters:
1. the value of the current semi-major axis, [MATH]a[/MATH];
2. the value of the current orbital eccentricity, [MATH]e[/MATH];
3. the value of 'GM' for the central gravitating body, [MATH]\mu[/MATH];
4. the current mass of the spacecraft, [MATH]m_0[/MATH];
5. the mass flow rate for the engines chosen to execute the retrograde burn, [MATH]\gamma[/MATH]; and
6. the effective exhaust velocity of the combustion gases of those engines [MATH]v_e[/MATH].
To obtain these values, you will of course need to interrogate Orbiter's system.
For the specific example given in the code below (which considers a Shuttle A on approach to perilune aiming to use its main engines to enter into a circular orbit. The algorithm reports that:
a. the duration of the orbit insertion burn is: 19.595 seconds;
b. the orbit insertion burn should commence 10.277 seconds before periapsis; and
c. the radius of the final circular orbit will be 3380913.4 metres.
In comparison, BTC reports that the duration of the orbit insertion burn is: 19.591 seconds; and the insertion burn should commence 9.831 seconds before periapsis. Here, the algorithm opts for a slightly longer burn time because the final orbital radius is very slightly lower than that assumed by BTC. The material difference between the two solutions is that the algorithm estimates that the burn should commence 0.45 seconds earlier than that calculated by BTC. For fast, high-thrust orbit insertion burns such as this the differences in solutions are not particularly material. But for low-thrust orbit insertion burns, the differences are much greater.
At this stage, to hook this up to BTC, is there anything else that you need?
Code:
// main.c
// A procedure to use 'automatic differentiation' to allow fast computation of the derivtaives
// to arbitrarily high order - first used in celestial mechanics problems
// For further details see, for example, http://www.maia.ub.edu/~angel/taylor/taylor.pdf
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// parameters used to define the initial orbit of the spacecraft
#define MU 4.9048695e12 // the 'GM' parameter for the central gravitating body (m^3 s^-2)
#define SMAJ 84781532.99 // the semi-major axis of the spacecraft's current orbit (m)
#define ECC 0.960122 // the orbital eccentricity of the spacecraft's current orbit (n/a)
// parameters used to define the initial state of the spacecraft
#define M0 87200.0 // the spacecraft's initial mass (kg)
#define GAMMA 64.507 // the fuel mass fow rate (kg/s)
#define VE 33000.0 // the effective velocity of engine exhaust gases (m/s)
// parameters used to define numerical integration order and step-size
#define DELTAT 20.0 // the integration step-size (s)
#define N 20 // the order of the Taylor series expansion of the integration
int taylorIntegation( double er[N+1], // time derivatives of the radial component of the orbital eccentricity vector
double et[N+1], // time derivatives of the transbsrese component of the orbital eccentricity vector
double th[N+1], // time derivatives of the angular component of the spacecraft's position
double h [N+1], // time derivatives of the spacecrafts specific angular momentum
double f [N+1] // time derivatives of the spacecraft's acceleration
) {
// declare the working arrays
double u01[N], u02[N], u05[N], u06[N], u07[N], u08[N], u09[N], u10[N], u11[N], u12[N], u13[N], u14[N], u15[N];
double v01[N], v02[N];
double w01[N], w02[N], w04[N], w05[N], w06[N], w07[N], w08[N], w09[N], w10[N];
// declare miscellaeous
double alpha;
// the integration step - building Taylor series expansions up to 20th order.
for ( int n = 0; n < N; n++) {
u01[n] = er[n];
u02[n] = et[n];
v01[n] = f [n];
w01[n] = h [n];
// calculate f'(t) - (the derivatives of the spacecraft's thrust function)
v02[n] = 0;
for (int j = 0; j <= n; j++) { // f(t)^2
v02[n] += v01[j] * v01[n-j];
}
v02[n] = v02[n] / VE; // f(t)^2 / ve
f[n+1] = v02[n] / (n+1); // f'(t) = f(t)^2 / ve
// calculate h'(t) and higher order derivatives
w02[n] = 0;
for (int j = 0; j <= n; j++) { // h(t)^2
w02[n] += w01[j] * w01[n-j];
}
w04[n] = 0;
for (int j = 0; j <= n; j++) { // f(t) * h(t)^2
w04[n] += w02[j] * v01[n-j];
}
w05[n] = u01[n];
if (n == 0) { // 1 + e_r(t)
w05[n] += 1;
}
w06[n] = 0;
for (int j = 0; j <= n; j++) { // (1 + e_r(t))^2
w06[n] += w05[j] * w05[n-j];
}
w07[n] = 0;
for (int j = 0; j <= n; j++) { // e_t(t)^2
w07[n] += u02[j] * u02[n-j];
}
w08[n] = w06[n] + w07[n]; // A = (1 + e_r(t))^2 + e_t(t)^2
alpha = -0.5;
if (n == 0) {
w09[n] = pow( w08[0], alpha ); // A^-1/2
} else {
w09[n] = 0;
for (int j = 0; j <= n-1; j++) {
w09[n] += (n * alpha - j * (1 + alpha)) * w08[n-j] * w09[j];
}
w09[n] = w09[n] / n / w08[0];
}
w10[n] = 0;
for (int j = 0; j <= n; j++) {
w10[n] += w09[j] * w04[n-j];
}
w10[n] = - w10[n] / MU; // - mu^-1 * f(t) * h(t)^2 * A^-1/2
h[n+1] = w10[n] / (n+1);
// calculate theta'(t) and higher order derivatives
alpha = -3.0;
if (n == 0) {
u05[n] = pow( w01[0], alpha ); // h(t)^-3
} else {
u05[n] = 0;
for (int j = 0; j <= n-1; j++) {
u05[n] += (n * alpha - j * (1 + alpha)) * w01[n-j] * u05[j];
}
u05[n] = u05[n] / n / w01[0];
}
u06[n] = 0;
for (int j = 0; j <= n; j++) {
u06[n] += u05[j] * w06[n-j];
}
u06[n] = u06[n] * MU * MU; // mu^2 * (1 + e_r(t))^2 * h(t)^-3
th[n+1] = u06[n] / (n+1); // th'(t) = u06 / (n+1)
// calculate e_r'(t) and higher order derivatives
u07[n] = 0;
for (int j = 0; j <= n; j++) {
u07[n] += u06[j] * u02[n-j]; // mu^2 * (1 + e_r(t))^2 * h(t)^-3 * e_t(t)
}
u08[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t)^2 * A^-1/2 * (1 + e_r(t))
u08[n] += w10[j] * w05[n-j];
}
alpha = -1.0;
if (n == 0) {
u09[n] = pow( w01[0], alpha ); // h(t)^-1
} else {
u09[n] = 0;
for (int j = 0; j <= n-1; j++) {
u09[n] += (n * alpha - j * (1 + alpha)) * w01[n-j] * u09[j];
}
u09[n] = u09[n] / n / w01[0];
}
u10[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t) * A^-1/2 * (1 + e_r(t))
u10[n] += u09[j] * u08[n-j];
}
u11[n] = + u07[n] + 2 * u10[n]; // mu^2*(1+e_r(t))^2*h(t)^-3*e_t(t) - 2*mu^-1*f(t)*h(t)*A^-1/2*(1 + e_r(t))
er[n+1] = u11[n] / (n+1); // e_r'(t) = u11 / (n+1)
//calculate e_t'(t) and higher order derivatives
u12[n] = 0;
for (int j = 0; j <= n; j++) { // mu^2 * (1 + e_r(t))^2 * h(t)^-3 * e_r(t)
u12[n] += u06[j] * u01[n-j];
}
u13[n] = 0;
for (int j = 0; j <= n; j++) {
u13[n] += w10[j] * u02[n-j]; // - mu^-1 * f(t) * h(t)^2 * A^-1/2 * e_t(t)
}
u14[n] = 0;
for (int j = 0; j <= n; j++) { // - mu^-1 * f(t) * h(t) * A^-1/2 * e_t(t)
u14[n] += u13[j] * u09[n-j];
}
u15[n] = - u12[n] + 2 * u14[n]; // - mu^2*(1 + e_r(t))^2*h(t)^-3*e_r(t) - 2*mu^-1*f(t)*h(t)*A^-1/2*e_t(t)
et[n+1] = u15[n] / (n+1); // e_t'(t) = u15 / (n+1)
}
return 0;
}
int updateInitialValues( double er[N+1], // time derivatives of the radial component of the orbital eccentricity vector
double et[N+1], // time derivatives of the transbsrese component of the orbital eccentricity vector
double th[N+1], // time derivatives of the angular component of the spacecraft's position
double h [N+1], // time derivatives of the spacecrafts specific angular momentum
double f [N+1], // time derivatives of the spacecraft's acceleration
double dt // the time step
) {
double temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + er[j]);}; er[0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + et[j]);}; et[0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + th[j]);}; th[0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + f [j]);}; f [0] = temp;
temp = 0; for (int j = N; j >=0; j--) {temp = (temp * dt + h [j]);}; h [0] = temp;
return 0;
}
int f ( double R, // the orbital radius at the start of the orbit insertion burn
double tm, // the duration of the orbit insertion burn
double e[4], // the radial and transverse components of the eccentricity vector at the end of the orbital inserion burn
double g[3]
) {
// declare result arrays - normalised Taylor series derivatives up to and including 10th order
double er[N+1]; // the component of the eccentricity vector in the radial direction
double et[N+1]; // the component of the eccentricity vector in the 'theta' direction
double th[N+1]; // the spatial angular coordinate of the spacecraft
double h [N+1]; // the specific angular momentum of the spacecraft
double f [N+1]; // the spacecraft's thrust function g * ve / (m0 - g * t)
// initalise the values to start the integration processteration process
er[0] = SMAJ * (1 - ECC * ECC) / R - 1;
et[0] = sqrt( (1 - ECC * ECC) * (SMAJ * ECC + SMAJ - R) * (SMAJ * ECC - SMAJ + R)) / R;
th[0] = 0.0;
h [0] = sqrt(MU * SMAJ * (1 - ECC * ECC));
f [0] = GAMMA * VE / M0;
// declare miscellanous quantities
int status;
double t, temp;
t = tm;
// carry out the integration step nine times - with each integration step spanning 1000 seconds
while (t >= DELTAT) {
status = taylorIntegation (er, et, th, h, f);
status = updateInitialValues (er, et, th, h, f, DELTAT);
t -= DELTAT;
}
status = taylorIntegation (er, et, th, h, f);
status = updateInitialValues (er, et, th, h, f, t);
// calculate the components of the eccentricity vector at the end of the orbit insertion burn
e[0] = er[0];
e[1] = et[0];
// calculate the firs derivatives with respect to time of the eccentricity vector at the end of the orbit insertion burn
temp = 0; for (int j = N; j >=1; j--) {temp = (temp * t + j * er[j]);}; e[2] = temp;
temp = 0; for (int j = N; j >=1; j--) {temp = (temp * t + j * et[j]);}; e[3] = temp;
g[0] = th[0];
g[1] = f [0];
g[2] = h [0];
return 0;
}
int main(int argc, const char * argv[]) {
// INPUTS
double r = 3380938.0; // the orbital radius at the start of the orbit insertion burn
double tm = 19.586; // the duration of the orbit insertion burn
double dr = 1.0; // the small distance displacement for calculating derivatives of the final eccentricity vector with resepct to r;
// OUTPUTS
double e[4]; // the radial and transverse components of the eccentricity vector at the end of the orbital insertion burn
// e[0] holds the radial component of the eccentricity vector
// e[1] holds the transverse component of the eccentricity vector
// e[2] holds the derivative of the radial component with respect to time
// e[3] holds the derivative of the transverse component with respect to time
double g[3]; // a vector to hold other useful information to be used once the Newton-Raphson iteration has been completed
// g[0] holds the the angle variable after orbit insertion
// g[1] holds the spacecraft acceleration due to thrust immediatel prior to the end of the insertion burn
// g[2] holds the spacecraft's specific angular momentum
double f0[2]; // a 2-vector declaration to hold the values of the orbital eccentricity vector at the current N-R iteration
double jac[2][2]; // a 2 x 2 array declaration to hold the Jacobian
double det; // the determinant of the Jacobian
int flag = 1;
// perform the Newton-Raphson (N-R) root-finding algorithm
int status;
while ((sqrt(f0[0] * f0[0] + f0[1] * f0[1]) > 1.e-8) || flag == 1) {
flag = 0;
status = f(r + dr, tm, e, g);
jac[0][0] = e[0];
jac[1][0] = e[1];
status = f(r - dr, tm, e, g);
jac[0][0]-= e[0]; jac[0][0] = jac[0][0] / dr / 2;
jac[1][0]-= e[1]; jac[1][0] = jac[1][0] / dr / 2;
status = f(r, tm, e, g);
jac[0][1] = e[2];
jac[1][1] = e[3];
f0[0] = e[0];
f0[1] = e[1];
det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
r = r - (+ jac[1][1] * f0[0] - jac[0][1] * f0[1]) / det;
if (r <= SMAJ * (1 - ECC)) {
r = SMAJ * (1 - ECC) + 100;
}
if (if (r >= SMAJ * (1 + ECC) && a > 0) {) {
r = SMAJ * (1 + ECC) - 100;
}
tm = tm - (- jac[1][0] * f0[0] + jac[0][0] * f0[1]) / det;
}
status = f(r, tm, e, g);
f0[0] = e[0];
f0[1] = e[1];
// print the solution
// 'r' is the orbital radius of the spacecraft at the start of the orbit insertion burn (this can be converted to a time measure.)
// 'tm' is the duration (seconds) of the orbit insertion burn
// 'ecc' is the final orbital eccentricity at the end of the orbit insertion burn. This should be zero for a circular orbit
// 'r_c' is the orbital radius at the end of the burn (metres). For entry into a circular orbit, it also the radius of the final circular orbit.
// 'tbp' is the time before (current) periapsis when the orbit insertion burn should start.
printf("r %2.10f\n", r );
printf("tm %2.10f\n", tm );
printf("\n");
printf("ecc %2.10f\n", sqrt(f0[0] * f0[0] + f0[1] * f0[1]) );
printf("r_c %2.10f\n", g[2] * g[2] / MU / (1 + f0[0]));
printf("tbp %2.10f\n", sqrt(SMAJ/MU)*(sqrt(SMAJ*SMAJ*(-1 + ECC*ECC) + 2*SMAJ*r - r*r) - 2*SMAJ*atan(sqrt((SMAJ*(-1 + ECC) + r)/(SMAJ + SMAJ*ECC - r)))));
printf("\n");
return 0;
}
Last edited:
Enjo
Mostly harmless
- Joined
- Nov 25, 2007
- Messages
- 1,665
- Reaction score
- 13
- Points
- 38
- Location
- Germany
- Website
- www.enderspace.de
- Preferred Pronouns
- Can't you smell my T levels?
Excellent. I will refactor it then, and include this in BTC.
Of course there will be small differences for short burns, but since it will be possible to simulate low thrust burns, it naturally extends the usage of BTC and TransX.
Of course there will be small differences for short burns, but since it will be possible to simulate low thrust burns, it naturally extends the usage of BTC and TransX.
Enjo
Please note that I have amended the code of the post time-stamped 11-02-16 04:39 AM to fix two minor things:
1. a lingering reference to 'a' which should have been 'SMAJ'; and
2. the calculation of the 'tbp' - time before periapsis. Previously, it was correct only for elliptical orbits. I have now added another line that calculates 'tbp' for hyperbolic orbits.
Also, just to emphasise: the code is currently only applicable to retrograde circularisation burns 'close to' orbital periapsis. Although the algorithm is essentially the same for prograde circular burns close to apoapsis, I will need to make some tweaks to make it work in this case. I would suggest that you just focus on periapsis circularisation first. And if that works well, I will provide modified code suitable for prograde apoapsis circularisation burns.
Please note that I have amended the code of the post time-stamped 11-02-16 04:39 AM to fix two minor things:
1. a lingering reference to 'a' which should have been 'SMAJ'; and
2. the calculation of the 'tbp' - time before periapsis. Previously, it was correct only for elliptical orbits. I have now added another line that calculates 'tbp' for hyperbolic orbits.
Also, just to emphasise: the code is currently only applicable to retrograde circularisation burns 'close to' orbital periapsis. Although the algorithm is essentially the same for prograde circular burns close to apoapsis, I will need to make some tweaks to make it work in this case. I would suggest that you just focus on periapsis circularisation first. And if that works well, I will provide modified code suitable for prograde apoapsis circularisation burns.
Last edited:
jedidia
shoemaker without legs
Guys, if you finally bring feasible navigation for low-thrust trajectories to orbiter, I'll love you so much!
OK, Jedidia
I assume you mean going somewhat beyond low-thrust orbit insertion - say Earth to Mars transit?
And when you say low thrust, can you give an example of thrust characteristics - i.e.,
m_0 - mass of spacecraft;
gam - mass flow-rate of craft
ve - effective exhaust velocity of propellant
Cheers
I assume you mean going somewhat beyond low-thrust orbit insertion - say Earth to Mars transit?
And when you say low thrust, can you give an example of thrust characteristics - i.e.,
m_0 - mass of spacecraft;
gam - mass flow-rate of craft
ve - effective exhaust velocity of propellant
Cheers
Last edited:
Enjo
Mostly harmless
- Joined
- Nov 25, 2007
- Messages
- 1,665
- Reaction score
- 13
- Points
- 38
- Location
- Germany
- Website
- www.enderspace.de
- Preferred Pronouns
- Can't you smell my T levels?
I knew YOU would come up
I might need some time to work on this, as my wife returns from her business trip tomorrow and will stay for 2 weeks. I plan to code this from time to time in this time span.
jedidia
shoemaker without legs
Well, any transit, really. It doesn't have to be a complete trajectory solver too, just a tool that can work with low acceleration for ejection and insertion, allowing to plan a low-thrust trajectory with some fiddling (i.e. a tool that can resolve the trajectory based on initial input of a delta-v vector, but not necessarily calculate an ideal trajectory from nothing but a source and a target). We almost got there some years ago, but mindblast lost interest in the project, and it wasn't open-source, so it's pretty much out of date by now.
As flexible as possible, but generally when I say "low thrust" I'm thinking in a range between a mili-G and a centi-G.
Oh did you? Is my desire for practical low-thrust navigation in Orbiter that well-known already? :lol:
Well, I just always thought that it was a pitty that the most realistic spaceflight simulator on the market didn't have any propper tools to simulate the most realistic future way to get from planet to planet...
I would be long busy doing this myself, if only my math wasn't stuck at lower high-school level (I tried to improve it once, but it didn't work out...). If there's any non-math related assistance I can lend, just give me a call.
Enjo
Mostly harmless
- Joined
- Nov 25, 2007
- Messages
- 1,665
- Reaction score
- 13
- Points
- 38
- Location
- Germany
- Website
- www.enderspace.de
- Preferred Pronouns
- Can't you smell my T levels?
Yes, in the TransX dev thread, but I couldn't help with that.
I didn't think that it would be as hard as the proposed solution, but it seems that Keith eats math books for breakfast.
---------- Post added at 04:27 PM ---------- Previous post was at 06:33 AM ----------
Keith:
Something's not right. The function is supposed to return the time to burn, based on DeltaV needed and orbital / spacecraft parameters, while in your code you're dependent on some tm, which you describe as "the duration of the orbit insertion burn". This is exactly the value we want to learn, not pass into the algo
Also, in your example, the following routine is not executed at all:
PHP:
// carry out the integration step nine times - with each integration step spanning 1000 seconds
while (t >= DELTAT) {
status = taylorInt.Run(er, et, th, h, f);
status = taylorInt.updateInitialValues(er, et, th, h, f, DELTAT);
t -= DELTAT;
}[EDIT]
or is it so, that we're still dependent on the current rocket equation and use its result to refine it?
Last edited:
Yes, it is a process of refinement. The underlying process for finding the burn start time and duration is an iterative one. This means that given a first 'best guess', the algorithm provides a better guess; and then, based on than new, 'best guess' it calculates an even better guess. And so on until according to some criterion, the algorithm says that the guess is 'good enough' and stops. In the code, the stopping criterion is that the final insertion orbit has an eccentricity of less than 10^-8.
So, what does this mean for BTC? Let's suppose that BTC starts is powered up and starts its calculations with a refresh rate of, say, 0.1 seconds. On the very first calculation, BTC's calculation produces a reasonable guess of burn start time and duration. This initial 'best guess' is passed to the algorithm that I have written for refinement. It runs and produces an improved guess that is then reported, presumably, on the MFD display. But on subsequent cycles, the 'best guess' is the previously calculated value, not the value that BTC would have calculated. Essentially, BTC's current algorithm is used to bootstrap the new process. Or, at least, that's how I imagine that it will work.
With respect to the code block:
PHP:
// carry out the integration step nine times - with each integration step spanning 1000 seconds
while (t >= DELTAT) {
status = taylorInt.Run(er, et, th, h, f);
status = taylorInt.updateInitialValues(er, et, th, h, f, DELTAT);
t -= DELTAT;
}Yes, it does run. To test I put a dummy printf statement in the code block and, sure enough, the printf statement repeatedly printed the informative message "Got to here!". It runs, of course, because the algorithm is passed some initial guesses for 'tm' and 'R'. These are:
PHP:
int main(int argc, const char * argv[]) {
// INPUTS
double r = 3380938.0; // the 'best guess' of the orbital radius at the start of the orbit insertion burn
double tm = 19.586; // the 'best guess' duration of the orbit insertion burnDoes that help?
Last edited:
Enjo
Mostly harmless
- Joined
- Nov 25, 2007
- Messages
- 1,665
- Reaction score
- 13
- Points
- 38
- Location
- Germany
- Website
- www.enderspace.de
- Preferred Pronouns
- Can't you smell my T levels?
Alright. I get the idea of initial guess, but I don't like this bit:
Could we make a deal, that you prepare another version in such a way, that you don't care how long it takes to compute - meaning that you get the initial guess at every iteration, and return the refined value at every iteration? I will test it in BTC and it this doesn't work, we'll get back to the version which spares CPU.
The more general idea is that it's always good to treat complex stuff like you created as "black boxes", which are fed only a series of inputs, and return one output, ideally not storing any state information (refined initial guess) between calls.
This is also loosely related with the PEG discussion - yes, it was written a long time ago, but PEG as a whole means separating calculation loops into major (less frequent) cycles and minor (more frequent) cycles, and that's why it has been given an antique badge. There's no such limitation now. In Launch MFD I find roots of major + minor loops of PEG on every MFD refresh.
This idea comes from underestimation of today's processing power. I could almost guarantee you, that after you integrated via Taylor's series, a modern PC will handle the rest smoothly. After all, the only major iterative process will be root finding, but it's not a naive search either. To get an idea of how fast you can calculate, please have a look at TransX's Auto-Min feature. It calculates Nelder-Mead optimization 10 times per MFD refresh (because this is how much is needed for convergence), and that can be calculated ~100 times per second, so making a total of 1000 iterations per second.
Could we make a deal, that you prepare another version in such a way, that you don't care how long it takes to compute - meaning that you get the initial guess at every iteration, and return the refined value at every iteration? I will test it in BTC and it this doesn't work, we'll get back to the version which spares CPU.
The more general idea is that it's always good to treat complex stuff like you created as "black boxes", which are fed only a series of inputs, and return one output, ideally not storing any state information (refined initial guess) between calls.
This is also loosely related with the PEG discussion - yes, it was written a long time ago, but PEG as a whole means separating calculation loops into major (less frequent) cycles and minor (more frequent) cycles, and that's why it has been given an antique badge. There's no such limitation now. In Launch MFD I find roots of major + minor loops of PEG on every MFD refresh.
Last edited:
Actually, I've given no thought to how fast modern CPUs execute. Generally, I presume them to be fast enough for current purposes.
The issue is: "What is the best available first guess?". At the outset, the best available first guess is the value calculated by the current BTC algorithm; subsequent best first guesses are the previous outputs of the routine since orbital parameters will have changed negligibly in 0.1 seconds or so (or whatever the MFD refresh rate is). But, if you wish, feel free to repeatedly use BTC's current algorithm as a starting point.
At the moment, I can see no need for another version of the code. Whatever policy re first guesses you want to use, the algorithm works 'as is' and will already produce a refined value each time the MFD is refreshed (i.e., every time the algorithm is run.). You just need to choose a 'first guess' policy.
Sorry I can't be more helpful, but I'm a little baffled by your comment. We may well be talking at cross-purposes here.
---------- Post added at 10:11 AM ---------- Previous post was at 09:21 AM ----------
Jedida
Below is a scenario that tests the algorithms ability to calculate orbit insertion for a low thrust 'centi-G' Earth orbit insertion. The scenario has a Shuttle PB, a stock Orbiter vessel, approaching Earth with a hyperbolic excess velocity of 5 km/s. The Shuttle PB uses its RCS to provide the retrograde thrust.
According to the algorithm described above, the RCS orbit insertion burn should start 17005.7 seconds before current orbital periapsis when the Shuttle PH has an orbital radius of 104,424 km. The total burn time is 45,863 seconds (i.e., 12 hours 44 minutes and 23 seconds). And the radius of the final circular orbit is 377 km. During the burn, the Shuttle A will complete two orbits around the Earth.
If you want to run this scenario, engage the prograde autopilot; switch engines to "RCS RE" in BTC; enter a Delta-V amount of 7.795 m/s; and when the Time to Periapsis reaches 17.00k seconds, engage the RCS by hitting the "BRN" button in BTC. Although executing this scenario is akin to watching paint dry, the Shuttle PB will insert itself into the circular orbit with some fidelity.
Overall, I think the current algorithm can calculate orbit insertion burns for 'centi-G' level thrusts. And there is no reason to believe it can't do the same for orbit escape burns as well. However, for mill-G thrusts, the gravity model in the code may well need to include other gravitating bodies since the insertion burn may start outside or at the edge of Earth's SOI at 0.5 million km or more. This far out, the Sun's gravitational field becomes a factor - as does, indeed, the Moon's.
Code:
BEGIN_DESC
END_DESC
BEGIN_ENVIRONMENT
System Sol
Date MJD 52007.2633725410
Help CurrentState_img
END_ENVIRONMENT
BEGIN_FOCUS
Ship SHU
END_FOCUS
BEGIN_CAMERA
TARGET SHU
MODE Cockpit
FOV 40.00
END_CAMERA
BEGIN_HUD
TYPE Orbit
REF AUTO
END_HUD
BEGIN_MFD Left
TYPE Orbit
PROJ Ship
FRAME Ecliptic
REF Earth
END_MFD
BEGIN_MFD Right
TYPE User
MODE BurnTimeMFD
END_MFD
BEGIN_SHIPS
SHU:ShuttlePB
STATUS Orbiting Earth
RPOS 13103960.028 -0.001 -105731779.682
RVEL 1625.1266 0.0000 5462.7410
AROT 0.000 -16.568 -0.000
AFCMODE 7
PRPLEVEL 0:1.000000
NAVFREQ 0 0
END
END_SHIPS
Last edited:
Longjap
Active member
- Joined
- Jun 8, 2011
- Messages
- 191
- Reaction score
- 41
- Points
- 28
Sorry to chime in guys but jedidia, does this wish also include low-thrust trajectory planning by means of solar sails? If this is a possibility to implement I'll double on that love to the Keithth's and Enjo's of this world.
