The Two-Point Boundary-Value Problem for Rocket Trajectories

Abstract

The two dimensional gravity turn problem is addressed allowing for the effects of variable rocket mass due to propellant consumption, thrust and thrust vector angle, lift and drag forces at an angle-of-attack and atmospheric mass density varying with altitude; Coriolis and centrifugal forces are neglected. Three distinct analytical solutions are obtained for constant: propellant flow rate, thrust, thrust vector angle, angle-of-attack and acceleration of gravity; the lift and drag are assumed to be proportional to the square of velocity, and the mass density is assumed to decrease exponentially with altitude. The method III uses power series of time for the horizontal (downrange) and vertical (altitude) coordinates; the method II replaces the altitude as variable by the atmospheric mass density and method I by its inverse. Thus the three solutions have distinct properties, e.g., I and III converge best close to lift-off and II close to burn-out. The three solutions: I, II, III, can be applied in isolation (or matched in combination) to the single-point boundary-value problem (SPBVP) of finding the trajectory with given initial conditions at launch. They can also be used as pairs in six distinct ways (I + II, I + III, II + III or reverse orders) to solve the two-point boundary-value problem (TPBVP), viz.: from given conditions at launch achieve one (not more) specified condition at burn-out, e.g., ã desired horizontal velocity for payload release. Each of the six distinct combinations of methods of addressing the TPBVP shares three features: (i) it can determine if there is a solution, viz. if the rocket has enough performance to reach the desired burn-out condition; (ii) if the desired burn-out condition is achievable it can calculate the complete trajectory from launch to burn-out; (iii) it can determine the range of achievable burn-out conditions, e.g., the minimum and maximum possible horizontal velocity at burn-out for given initial conditions at launch.

1. Introduction

The optimization [1] and reconstruction [2] of rocket trajectories is a research area with renewed interest for small satellite launchers [3,4], including dynamics [5] and control [6] aspects. Flight Dynamics [7,8,9] is the basis for the optimization of trajectories of single-stage [10,11] and multistage [12] rockets based trajectory arcs [13]: the trajectory arcs are calculated by numerical [14] or analytical [15,16] methods and applied to expendable [17,18] and reusable reusable [19] launchers. The trajectory optimization [10,11,12,13,14,15,16,17,18,19] may be coupled to control optimization [20].

In an earlier paper [21] the calculation of rocket trajectories in the atmosphere was considered taking into account the following effects: (i) the variable mass associated with propellant consumption; (ii) the effect of thrust at an angle to the rocket axis; (iii) the lift and drag for flight at an angle-of-attack; (iv) the proportionality of the aerodynamic forces on the square of the velocity and on the atmospheric mass density; (v) the dependence of the atmospheric mass density on the altitude. The latter effect has not been considered in most of the literature on analytic calculation of rocket trajectories [22,23,24,25,26,27,28,29]: in particular, if time is eliminated to express the trajectory in terms of velocity and flight path angle, the inclusion of the dependence of aerodynamic forces on altitude through the atmospheric mass density, introduces altitude explicitly as a third variable, related to the velocity and flight path angle by an integral relation.

A second difficulty with using velocity and flight path angle as variables, is that the latter is undetermined if the initial velocity is zero. This difficulty can be overcome by taking the initial time a little later than lift-off, e.g., at the start of the gravity turn; however, if the velocity is still small and the flight path angle imprecise, this will introduce an initial error, which can propagate and accumulate through the whole trajectory calculation. Thus, in the earlier paper [21] the equations of motion were written in an earth fixed reference frame, which has two advantages: (i) there is no singularity even for zero initial velocity, which is actually simpler than starting with initial motion, although that case is obviously also included; (ii) the dependence of the aerodynamic forces on altitude through the atmospheric mass density does not add a new variable to the problem. It introduces another nonlinearity, viz: (i) the first nonlinearity is the dependence of aerodynamic forces on the velocity, viz. on the square in the simplest case; (ii) the second nonlinearity is the linear dependence of the aerodynamic forces on the atmospheric mass density, which is a nonlinear function of altitude, e.g., an exponential function in the case of an isothermal atmospheric layer.

For the equations of motion in Cartesian coordinates, the components of the acceleration are linear, and the equations of motion allow: (a) the thrust to depend on altitude, due to the variation of atmospheric pressure; (b) the thrust angle to the rocket axis to vary with time, e.g., vectored thrust; (c) the angle-of-attack to vary with time, i.e., pitch control; (d) the fuel flow rate to vary with time, i.e., engine throttling. For definiteness and simplicity these four quantities are taken as constant in most calculations, although they could be varied step-by-step in a numerical code. The effects neglected altogether are: ($\alpha$) the variation of the acceleration due to gravity over an altitude range small compared with the Earth’s radius; ($\beta$) the Coriolis and centrifugal forces, which become significant at larger velocities. In a preceding paper [21] four methods of trajectory calculation were discussed in detail, with a minimum of algebraic clutter, by considering the simplest case of a vertical climb, without lift and with thrust aligned to the rocket axis; since account was taken of mass dependence on time, and drag dependence on velocity and altitude, similar methods apply in the more general case considered in the present paper: a gravity turn, with lift and drag allowing for an angle-of-attack and thrust vector at an angle to the rocket axis.

The Single-Point Boundary-Value Problem (SPBVP) is extended from the one-dimensional vertical ascent trajectory in the preceding paper [21] to the two dimensional gravity turn in the present paper, taking the equations of motion in a matrix form (Section 2.1), to apply three similar methods of solution: (Section 2.1) method III expressing altitude and downrange as power series of time; (Section 2.3) method II using altitude as a variable to the atmospheric mass density, the latter as a function of burned propellant mass fraction; (Section 2.2) method I using the inverse ratio of atmospheric mass densities as a function of remaining propellant mass fraction. The general approach to the two-point boundary-value problem (TPBVP) addresses (Section 3) three issues: (Section 3.1) existence of solutions; (Section 3.2) calculation of the trajectory; (Section 3.3) performance envelope at burn-out. This general approach can be applied to any combination of the three preceding methods, viz. a total of three combinations I + II (Section 4), III + II (Section 5) and III + I (Section 6); the number of combinations is doubled to six by choosing one method for the “upward” solution from lift-off and another for the “downward” solution from the burn-out, with matching of both coordinates (altitude and downrange) and velocities (horizontal and vertical) at an intermediate time. By choosing the matching time it is possible to specify one (not, more) burn-out condition (Section 7), e.g., the desired horizontal velocity for payload launch.

2. Three Methods of Calculation of Gravity Turn

The equations of motion in an earth-fixed Cartesian frame are solved for a gravity turn by three methods: (III) direct expansion of the coordinates as Taylor series of time (Section 2.1); (I) use of the atmospheric mass density as as an altitude variable (Section 2.2) as a function of the remaining propellant mass fraction; (II) use of the inverse of the atmospheric mass density as altitude variable (Section 2.3) as a function of the burned propellant mass fraction.

2.1. Method III: Iterative Use of Initial Conditions in Taylor Series

The equations of motion of a rocket (Figure 1) in an earth-fixed reference frame take the form ([21], Equation (25)):

$$\begin{array}{cc}[m_0-c(t-t_0)]{\left({\dot{x}}^2+{\dot{z}}^2\right)}^{1/2}\left[\begin{array}{c}\ddot{x}\\ \ddot{z}+g\end{array}\right]=\hfill & \left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]\left[\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right]\hfill \\ & \hfill -\frac{1}{2}{\rho }_{0}Sexp[-(z-z_0)/\ell ]\left({\dot{x}}^2+{\dot{z}}^2\right)\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]\left[\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right],\end{array}$$

(1)

where:

(i)

the variation of mass with time, Equation (2b), corresponds to a constant, Equation (2a), propellant rate:

$$\begin{array}{cc}\hfill c& \equiv -\frac{\mathrm{d}m}{\mathrm{d}t}:\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill m\left(t\right)& =m_0-c(t-t_0);\hfill \end{array}$$

(2b)

(ii)

the acceleration of gravity is assumed to be uniform, and the altitude coordinate z is taken opposite to it; (iii) the thrust matrix, Equation (3a) is constant for thrust T, angle-of-attack $\alpha$ and angle $\epsilon$ of the thrust with the rocket axis:

$$\begin{array}{cc}\hfill T_{ij}& =T\left[\begin{array}{cc}cos(\alpha +\epsilon )& -sin(\alpha +\epsilon )\\ sin(\alpha +\epsilon )& \phantom{-}cos(\alpha +\epsilon )\end{array}\right]\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill F_{ij}& =\left[\begin{array}{cc}cos\alpha & \phantom{-}sin\alpha \\ sin\alpha & -cos\alpha \end{array}\right]\left[\begin{array}{cc}{C}_{\mathrm{D}}& \phantom{-}{C}_{\mathrm{L}}\\ C_{\mathrm{L}}& -C_{\mathrm{D}}\end{array}\right],\hfill \end{array}$$

(3b)

(iii)

the aerodynamic matrix (3b) is constant for constant angle-of-attack $\alpha$ and lift $C_{\mathrm{L}}$ and drag $C_{\mathrm{D}}$ coefficients; (iv) the aerodynamic forces are proportional to the square of the velocity and to the atmospheric mass density; (v) the latter decays exponentially (4b) on the scale height ℓ for an isothermal atmosphere with scale height (4a), viz.:

$$\begin{array}{cc}\hfill \ell & \equiv R\theta /g,\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill \rho \left(z\right)& ={\rho }_{0}exp\left[-(z-z_0)/\ell \right],\hfill \end{array}$$

(4b)

where $\theta$ is the temperature and R the gas constant.

The equations of motion (1) with Equation (3a) and (3b) are to be solved subject to initial conditions specifying the coordinates (5a) and velocity components (5b) at initial time:

$$\begin{array}{cc}\hfill t=t_0:\left\{x\left(t_0\right),z\left(t_0\right)\right\}& =\left\{x_0,z_0\right\},\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill \left\{\dot{x}\left(t_0\right),\dot{z}\left(t_0\right)\right\}& =\left\{u_0,w_0\right\}.\hfill \end{array}$$

(5b)

Instead of (5b) the total velocity $v_0$ and flight path angle ${\gamma }_0$ could be specified initially as:

$$\begin{array}{cc}\hfill v_0& =\sqrt{u_0^2+w_0^2},\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill tan{\gamma }_0& =w_0/u_0.\hfill \end{array}$$

(6b)

The flight path angle ${\gamma }_0$ is indeterminate in Equation (6b), for zero initial velocity $v_0=0 \mathrm{m}/\mathrm{s}$ in Equation (6a), whereas the initial conditions in Equation (5b) are always determinate. Substituting the initial conditions, Equations (5a) and (5b) in the equations of motion (1) specifies:

$$m_0{v}_0\left[\begin{array}{c}{\ddot{x}}_0\\ {\ddot{z}}_0+g\end{array}\right]=\left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]\left[\begin{array}{c}{u}_0\\ w_0\end{array}\right]-\frac{1}{2}{\rho }_{0}S{v}_0^2\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]\left[\begin{array}{c}{u}_0\\ w_0\end{array}\right],$$

(7)

the initial accelerations $({\ddot{x}}_0,{\ddot{z}}_0)$.

Differentiating the equations of motion (1) with regard to time:

$$\begin{array}{cc}[m_0-c(t-t_0)]\hfill & \left[\begin{array}{c}\stackrel{⃛}{x}\\ \stackrel{⃛}{z}\end{array}\right]-c\left[\begin{array}{c}\ddot{x}\\ \ddot{z}+g\end{array}\right]={\left({\dot{x}}^2+{\dot{z}}^2\right)}^{-3/2}(\dot{z}\ddot{x}-\dot{x}\ddot{z})\left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]\left[\begin{array}{c}\phantom{-}\dot{z}\\ -\dot{x}\end{array}\right]\hfill \\ & -\frac{1}{2}{\rho }_{0}Sexp[-(z-z_0)/\ell ]{\left({\dot{x}}^2+{\dot{z}}^2\right)}^{-1/2}\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]\left[\begin{array}{c}2{\dot{x}}^2\ddot{x}+{\dot{z}}^2\ddot{x}+\dot{x}\dot{z}\ddot{z}\\ 2{\dot{z}}^2\ddot{z}+{\dot{x}}^2\ddot{z}+\dot{x}\dot{z}\ddot{x}\end{array}\right]\\ & +\frac{1}{2}({\rho }_{0}S/\ell )exp[-(z-z_0)/\ell ]{\left({\dot{x}}^2+{\dot{z}}^2\right)}^{1/2}\dot{z}\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]\left[\begin{array}{c}\dot{x}\\ \dot{z}\end{array}\right],\hfill \end{array}$$

(8)

and setting $t=t_0$ specifies the $O\left(t^3\right)$ coefficients:

$$\begin{array}{cc}{m}_0\left[\begin{array}{c}{\stackrel{⃛}{x}}_0\\ {\stackrel{⃛}{z}}_0\end{array}\right]=c\left[\begin{array}{c}{\ddot{x}}_0\\ {\ddot{z}}_0+g\end{array}\right]+v_0^{-3}(w_0\hfill & {\ddot{x}}_0-u_0{\ddot{z}}_0)\left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]\left[\begin{array}{c}\phantom{-}{w}_0\\ -u_0\end{array}\right]\hfill \\ & -\frac{1}{2}{\rho }_{0}S{v}_0^{-1}\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]\left[\begin{array}{c}2u_0^2{\ddot{x}}_0+w_0^2{\ddot{x}}_0+u_0{w}_0{\ddot{z}}_0-v_0^2{u}_0{w}_0/\ell \\ 2w_0^2{\ddot{z}}_0+u_0^2{\ddot{z}}_0+u_0{w}_0{\ddot{x}}_0-v_0^2{w}_0^2/\ell \end{array}\right],\hfill \end{array}$$

(9)

in the Taylor series expansion:

$$\left[\begin{array}{c}x\left(t\right)\\ z\left(t\right)\end{array}\right]=\left[\begin{array}{c}{x}_0\\ z_0\end{array}\right]+(t-t_0)\left[\begin{array}{c}{u}_0\\ w_0\end{array}\right]+\frac{1}{2}{(t-t_0)}^2\left[\begin{array}{c}{\ddot{x}}_0\\ {\ddot{z}}_0\end{array}\right]+\frac{1}{6}{(t-t_0)}^3\left[\begin{array}{c}{\stackrel{⃛}{x}}_0\\ {\stackrel{⃛}{z}}_0\end{array}\right]+O\left({(t-t_0)}^4\right),$$

(10)

of coordinates of the trajectory versus time. This cubic approximation can be applied step-by-step along the trajectory, allowing the thrust $T\left(t\right)$, its angle with the rocket axis $\epsilon \left(t\right)$ and the angle-of-attack $\alpha \left(t\right)$ to vary at each step, as well as the scale height $\ell \left(t\right)$ in the mass density.

2.2. Method I: Mass Fraction of Burned Fuel as the Time Variable

The equations of motion (1) may be put in a dimensionless form using a dimensionless time variable, Equation (11a) (Equation (34b) in [21]) the mass of burned fuel as a fraction of lift-off mass and making both altitude, Equation (11b) ([21], Equation (34a)) and downrange distance $x-x_0$ dimensionless, Equation (11c), dividing by the atmospheric scale height, viz.:

$$\begin{array}{cc}\hfill \tau & \equiv 1-\frac{m\left(t\right)}{m_0}=c(t-t_0)/m_0,\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill \zeta & \equiv (z-z_0)/\ell ,\hfill \end{array}$$

(11b)

$$\begin{array}{cc}\hfill \chi & \equiv (x-x_0)/\ell .\hfill \end{array}$$

(11c)

Introducing the dimensionless variables Equation (11a)–(11c) and denoting by prime the derivatives with regard to $\tau$

$$\begin{array}{cc}\hfill \left\{{\chi }^{\prime },{\zeta }^{\prime }\right\}& \equiv \left\{\frac{\mathrm{d}\chi }{\mathrm{d}\tau },\frac{\mathrm{d}\zeta }{\mathrm{d}\tau }\right\}=\frac{m_0}{c\ell }\left\{\frac{\mathrm{d}x}{\mathrm{d}t},\frac{\mathrm{d}z}{\mathrm{d}t}\right\}=\frac{m_0}{c\ell }\{\dot{x},\dot{z}\},\hfill \end{array}$$

(11d)

$$\begin{array}{cc}\hfill \left\{{\chi }^{″},{\zeta }^{″}\right\}& \equiv \left\{\frac{{\mathrm{d}}^2\chi }{\mathrm{d}{\tau }^2},\frac{{\mathrm{d}}^2\zeta }{\mathrm{d}{\tau }^2}\right\}=\frac{m_0^2}{c^2\ell }\left\{\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2},\frac{{\mathrm{d}}^{2}z}{\mathrm{d}{t}^2}\right\}=\frac{m_0^2}{c^2\ell }\{\ddot{x},\ddot{z}\},\hfill \end{array}$$

(11e)

the equations of motion (1) become Equation (12b):

$$\begin{array}{cc}\hfill a& =m_0^{2}g/\left(c^2\ell \right):\hfill \end{array}$$

(12a)

$$\begin{array}{cc}(1-\tau )\left[\begin{array}{c}{\chi }^{″}\\ {\zeta }^{″}+a\end{array}\right]\hfill & ={\left({\chi }^{\prime 2}+{\zeta }^{\prime 2}\right)}^{-1/2}\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\left[\begin{array}{c}{\chi }^{\prime }\\ {\zeta }^{\prime }\end{array}\right]-e^{-\zeta }{\left({\chi }^{\prime 2}+{\zeta }^{\prime 2}\right)}^{1/2}\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{\chi }^{\prime }\\ {\zeta }^{\prime }\end{array}\right],\hfill \end{array}$$

(12b)

involving: (i) the gravity parameter, Equation (12a) (Equation (36a) in [21]) unchanged; (ii) the thrust parameter from Equation (36b) in [21], is replaced by the matrix:

$$\begin{array}{cc}\hfill b_{ij}& =\frac{m_0}{c^2\ell }{T}_{ij}=b\left[\begin{array}{cc}cos(\alpha +\epsilon )& -sin(\alpha +\epsilon )\\ sin(\alpha +\epsilon )& \phantom{-}cos(\alpha +\epsilon )\end{array}\right],\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill b& \equiv \frac{m_{0}T}{c^2\ell };\hfill \end{array}$$

(13b)

(iii) the aerodynamic parameter in [21], Equation (36c), is replaced, Equation (14a), by an aerodynamic matrix:

$$\begin{array}{cc}\hfill f_{ij}& =\frac{{\rho }_{0}S\ell }{2m_0}{F}_{ij}=f\left[\begin{array}{cc}cos\alpha & \phantom{-}sin\alpha \\ sin\alpha & -cos\alpha \end{array}\right]\left[\begin{array}{cc}{C}_{\mathrm{D}}& \phantom{-}{C}_{\mathrm{L}}\\ C_{\mathrm{L}}& -C_{\mathrm{D}}\end{array}\right],\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill m_{★}& ={\rho }_{0}S\ell ,\hfill \end{array}$$

(14b)

$$\begin{array}{cc}\hfill f& =\frac{m_{★}}{2m_0},\hfill \end{array}$$

(14c)

involving both the drag $C_{\mathrm{D}}$ and lift $C_{\mathrm{L}}$ coefficients, and the reference mass, Equation (14b) ([21], Equation (38a)). The reference mass Equation (14b) is the mass of a cylinder with cross-sectional area S, length equal to the scale height ℓ and uniform mass density ${\rho }_0$ equal to the value at sea level.

Using a dimensionless altitude variable, Equation (15a) ([21], Equation (9)) the atmospheric mass density at the initial altitude divided by that at an arbitrary altitude

$$\eta \equiv {\rho }_0/\rho \left(z\right)=exp\left[(z-z_0)/\ell \right]=e^{\zeta },$$

(15a)

and using Equations (15b) and (15c):

$$\begin{array}{cc}\hfill {\zeta }^{\prime }=& \frac{\mathrm{d}}{\mathrm{d}\tau }(log\eta )=\frac{{\eta }^{\prime }}{\eta },\hfill \end{array}$$

(15b)

$$\begin{array}{cc}\hfill {\zeta }^{″}=& \frac{{\mathrm{d}}^2}{\mathrm{d}{\tau }^2}(log\eta )=\frac{\mathrm{d}}{\mathrm{d}\tau }\left(\frac{{\eta }^{\prime }}{\eta }\right)=\frac{{\eta }^{″}}{\eta }-{\left(\frac{{\eta }^{\prime }}{\eta }\right)}^2,\hfill \end{array}$$

(15c)

the exponential nonlinearity in the equations of motion (12b) is replaced (see [21]), Equation (40a) and (40b)) by algebraic nonlinearities of degree up to three:

$$\begin{array}{c}(1-\tau )\eta \left[\begin{array}{c}{\chi }^{″}\eta \\ {\eta }^{″}\eta -{\eta }^{\prime 2}+a{\eta }^2\end{array}\right]={\left({\eta }^{\prime 2}+{\eta }^2{\chi }^{\prime 2}\right)}^{-1/2}{\eta }^3\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\left[\begin{array}{c}{\chi }^{\prime }\\ {\eta }^{\prime }\end{array}\right]\hfill \\ \hfill -{\left({\chi }^{\prime 2}{\eta }^2+{\eta }^{\prime 2}\right)}^{1/2}\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{\chi }^{\prime }\\ {\eta }^{\prime }\end{array}\right].\end{array}$$

(16)

Note that if $x-x_0=\mathrm{const}.$, or ${\chi }^{\prime }=0$, then Equation (16) reduces to [21], Equation (41) of vertical climb. The solution of Equation (16) is taken in the form of Equation (17b) ([21], Equation (49) for altitude and Equation (17a) for the downrange:

$$\begin{array}{cc}\hfill \chi \left(\tau \right)& =(u_0/v_{★})\tau +\sum _{n=2}^{\infty }{X}_n{\tau }^n,\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill \eta \left(\tau \right)& =1+(w_0/v_{★})\tau +\sum _{n=2}^{\infty }{Z}_n{\tau }^n,\hfill \end{array}$$

(17b)

since: (i) the initial position Equation (5a) corresponds to $\chi =0$ in Equation (11c) and $\eta =1$ in Equation (15a); (ii) the first order derivatives are given by Equation (17c),

$$\left\{\frac{\mathrm{d}\chi }{\mathrm{d}\tau },\frac{\mathrm{d}\eta }{\mathrm{d}\tau }\right\}=\left\{\frac{1}{\ell }\frac{\mathrm{d}x}{\mathrm{d}t},exp\left(\zeta \right)\frac{1}{\ell }\frac{\mathrm{d}z}{\mathrm{d}t}\right\}\frac{\mathrm{d}t}{\mathrm{d}\tau }=\frac{m_0}{c\ell }\left\{\dot{x},exp\left[(z-z_0)/\ell \right]\dot{z}\right\},$$

(17c)

and in the limit of initial time, Equation (17d),

$$\underset{\tau \to 0}{lim}\left\{\frac{\mathrm{d}\chi }{\mathrm{d}\tau },\frac{\mathrm{d}\eta }{\mathrm{d}\tau }\right\}=\frac{m_0}{c\ell }\{u_0,w_0\}=\frac{\{u_0,w_0\}}{v_{★}},$$

(17d)

specifies the first-order coefficients in Equations (17a) and (17b), and involves the reference velocity, Equation (18a) ([21], Equation (47b)), that appears explicitly in Equations (18b) and (18c):

$$\begin{array}{cc}\hfill v_{★}& \equiv c\ell /m_0,\hfill \end{array}$$

(18a)

$$\begin{array}{cc}\hfill (x-x_0)/\ell & =(t-t_0)u_0/\ell +\sum _{n=2}^{\infty }{X}_n{\left[c(t-t_0)/m_0\right]}^n,\hfill \end{array}$$

(18b)

$$\begin{array}{c}\hfill exp\left[(z-z_0)/\ell \right]=(t-t_0)w_0/\ell +\sum _{n=2}^{\infty }{Z}_n{\left[c(t-t_0)/m_0\right]}^n,\end{array}$$

(18c)

where in Equations (18b) and (18c) Equations (11a)–(11c) were used. The reference velocity $v_{★}$ in Equation (18a) corresponds to the linear momentum for the initial mass $m_0{v}_{★}=c\ell$ being equal to the propellant mass consumption per unit time multiplied by the atmospheric scale height. The remaining coefficients $(X_n,Z_n)$ in Equations (17a) and (17b) for $n=2,3,\dots$ are obtained by substitution in the equation of motion (16) and equating to zero the coefficients of successive powers ${\tau }^2,{\tau }^3,\dots$.

This process is illustrated by substituting the trajectory, Equations (17a) and (17b) in the equation of motion (16) up to order $\tau$, viz. the first equation becomes:

$$\begin{array}{cc}\hfill (1-\tau ){(1+w_0\tau /v_{★})}^2& (2X_2+6X_3\tau )\hfill \\ & =\left\{G^{-1/2}{(1+w_0\tau /v_{★})}^3\left[b_{11}{b}_{12}\right]-G^{1/2}\left[f_{11}{f}_{12}\right]\right\}\left[\begin{array}{c}{u}_0/v_{★}+2X_2\tau \\ w_0/v_{★}+2Z_2\tau \end{array}\right],\hfill \end{array}$$

(19a)

and the second:

$$\begin{array}{c}(1-\tau )\left\{{(1+w_0\tau /v_{★})}^2(2Z_2+6Z_3\tau )-(1+w_0\tau /v_{★}){(w_0/v_{★}+2Z_2\tau )}^2+a{(1+w_0\tau /v_{★})}^3\right\}\hfill \\ \hfill =\left\{G^{-1/2}{(1+w_0\tau /v_{★})}^3\left[b_{21}{b}_{22}\right]-G^{1/2}\left[f_{21}{f}_{22}\right]\right\}\left[\begin{array}{c}{u}_0/v_{★}+2X_2\tau \\ w_0/v_{★}+2Z_2\tau \end{array}\right],\end{array}$$

(19b)

where

$$\begin{array}{cc}\hfill G\equiv {\eta }^{\prime 2}+{\eta }^2{\chi }^{\prime 2}={(w_0/v_{★}+2Z_2\tau )}^2& +{(1+w_0\tau /v_{★})}^2{(u_0/v_{★}+2X_2\tau )}^2\hfill \\ & ={(v_0/v_{★})}^2+2\left[(w_0/v_{★}){(u_0/v_{★})}^2+2(X_2{u}_0+Z_2{w}_0)/v_{★}\right]\tau ,\hfill \end{array}$$

(20)

using Equation (6a). Equating powers of ${\tau }^0$ in Equations (19a) and (19b) yields:

$$\begin{array}{cc}\hfill 2X_2& =b_{11}{u}_0/v_0+b_{12}{w}_0/v_0-(v_0/v_{★}^2)(f_{11}{u}_0+f_{12}{w}_0),\hfill \end{array}$$

(21a)

$$\begin{array}{cc}\hfill 2Z_2& ={(w_0/v_{★})}^2-a+b_{21}{u}_0/v_0+b_{22}{w}_0/v_0-(v_0/v_{★}^2)\left(f_{21}{u}_0+f_{22}{w}_0\right),\hfill \end{array}$$

(21b)

and equating powers of $\tau$ yields:

$$\begin{array}{cc}\hfill 6X_3& =2X_2(1-2w_0/v_{★})+2(v_{★}/v_0)(b_{11}{X}_2+b_{12}{Z}_2)\hfill \\ \hfill & -2(v_0/v_{★})(f_{11}{X}_2+f_{12}{Z}_2)+3(w_0/v_0)(b_{11}{u}_0+b_{12}{w}_0)/v_{★}\hfill \\ \hfill & -\left[(w_0/v_0){(u_0/v_{★})}^2+2(X_2{u}_0+Z_2{w}_0)/v_0\right]\left[{\left(v_{★}/v_0\right)}^2(b_{11}{u}_0+b_{12}{w}_0)+(f_{11}{u}_0+f_{12}{w}_0)\right]/v_{★},\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill 6Z_3& ={(w_0/v_{★})}^2(w_0/v_{★}-1)+a(1-3w_0/v_{★})+2Z_2\hfill \\ \hfill & +2(b_{21}{X}_2+b_{22}{Z}_2)v_{★}/v_0-2(f_{21}{X}_2+f_{22}{Z}_2)v_0/v_{★}+3(w_0/v_0)(b_{21}{u}_0+b_{22}{w}_0)/v_{★}\hfill \\ \hfill & -\left[(w_0/v_0){(u_0/v_{★})}^2+2(X_2{u}_0+Z_2{w}_0)/v_0\right]\left[{\left(v_{★}/v_0\right)}^2(b_{11}{u}_0+b_{12}{w}_0)+(f_{11}{u}_0+f_{12}{w}_0)\right]/v_{★}.\hfill \end{array}$$

(22b)

From Equation (20) follow the first-order approximation, Equations (22c) and (22d), that were used in Equations (19a) and (19b) to derive Equations (21a), (21b), (22a) and (22b):

$$\begin{array}{cc}\hfill G^{1/2}& =v_0/v_{★}+\tau (v_{★}/v_0)\left[(w_0/v_{★}){(u_0/v_{★})}^2+2(X_2{u}_0+Z_2{w}_0)/v_{★}\right],\hfill \end{array}$$

(22c)

$$\begin{array}{cc}\hfill G^{-1/2}& =v_{★}/v_0-\tau {(v_{★}/v_0)}^3\left[(w_0/v_{★}){(u_0/v_{★})}^2+2(X_2{u}_0+Z_2{w}_0)/v_{★}\right].\hfill \end{array}$$

(22d)

Thus, the first four coefficients in the trajectory Equations (17a) and (17b) are specified by Equations (21a), (21b), (22a) and (22b). The solution Equations (17a) and (17b) converge best for small $\tau$, i.e., short time or mass of fuel burned is a small fraction of the initial mass.

2.3. Method II: Ratio of Atmospheric Mass Densities as Altitude Variable

A faster convergence close to burn-out conditions is obtained by: (i) using the residual mass fraction, Equation (23) (Equation (55) in [21]) as dimensionless independent variable instead of time

$$1⩾\mu =m\left(t\right)/m_0=1-c(t-t_0)/m_0=1-\tau ⩾m_1/m_0,$$

(23)

where $m_1=m_0-c(t_1-t_0)$ is the residual mass; (ii) replacing altitude as a dependent variable by Equation (24a) (Equation (57) in [21]) the ratio of atmospheric densities at arbitrary and initial altitudes

$$\xi \equiv \frac{1}{\eta }=\rho \left(z\right)/{\rho }_0=exp\left[-(z-z_0)/\ell \right]=e^{-\zeta },$$

(24a)

(iii) using as the other dependent variable the same Equation (11c), downrange distance normalized to the scale height, using as variable instead of the burned mass fraction in $\chi \left(\tau \right)$ the residual mass fraction $\mu =1-\tau$ in $\beta \left(\mu \right)=\chi \left(\tau \right)$ in Equation (24b), and evaluating the derivatives, Equations (24c)–(24f) using $\mathrm{d}\mu =-\mathrm{d}\tau$ from Equation (23):

$$\begin{array}{cc}\hfill \chi \left(\tau \right)& =\chi (1-\mu )=\beta \left(\mu \right)\hfill \end{array}$$

(24b)

$$\begin{array}{cc}\hfill {\chi }^{\prime }& \equiv \frac{\mathrm{d}\chi }{\mathrm{d}\tau }=-\frac{\mathrm{d}\chi }{\mathrm{d}\mu }=-\frac{\mathrm{d}\beta }{\mathrm{d}\mu }=-{\beta }^{\prime },\hfill \end{array}$$

(24c)

$$\begin{array}{cc}\hfill {\chi }^{″}& \equiv \frac{{\mathrm{d}}^2\chi }{\mathrm{d}{\tau }^2}=\frac{{\mathrm{d}}^2\chi }{\mathrm{d}{\mu }^2}=\frac{{\mathrm{d}}^2\beta }{\mathrm{d}{\mu }^2}={\beta }^{″},\hfill \end{array}$$

(24d)

$$\begin{array}{cc}\hfill {\zeta }^{\prime }& \equiv -\frac{\mathrm{d}}{\mathrm{d}\tau }(log\xi )=\frac{\mathrm{d}}{\mathrm{d}\mu }(log\xi )=\frac{{\xi }^{\prime }}{\xi },\hfill \end{array}$$

(24e)

$$\begin{array}{cc}\hfill {\zeta }^{″}& \equiv -\frac{{\mathrm{d}}^2}{\mathrm{d}{\tau }^2}(log\xi )=-\frac{{\mathrm{d}}^2}{\mathrm{d}{\mu }^2}(log\xi )=-\frac{{\xi }^{″}}{\xi }+{\left(\frac{{\xi }^{\prime }}{\xi }\right)}^2.\hfill \end{array}$$

(24f)

Substitution of Equations (23) and (24a)–(24f) in the equations of motion (12b) leads to Equation (25):

$$\begin{array}{cc}\mu \left[\begin{array}{c}-{\beta }^{″}\xi \\ -{\xi }^{″}\xi +{\xi }^{\prime 2}+a{\xi }^2\end{array}\right]={\xi }^2{\left({\xi }^{\prime 2}+{\xi }^2{\beta }^{\prime 2}\right)}^{-1/2}\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\hfill & \left[\begin{array}{c}{\beta }^{\prime }\\ {\xi }^{\prime }\end{array}\right]\hfill \\ & -\xi {\left({\xi }^{\prime 2}+{\xi }^2{\beta }^{\prime 2}\right)}^{1/2}\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{\beta }^{\prime }\\ {\xi }^{\prime }\end{array}\right].\hfill \end{array}$$

(25)

The solution is sought as power series similar to (Equation (60) in [21]) viz.:

$$\begin{array}{cc}\hfill \beta \left(\mu \right)& =\sum _{n=0}^{{}^{\infty }}{P}_n{\mu }^{n+p},\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill \xi \left(\mu \right)& =\sum _{n=0}^{{}^{\infty }}{Q}_n{\mu }^{n+q},\hfill \end{array}$$

(26b)

where the initial conditions $t=t_0$ imply: (i) for the coordinates Equation (5a), then $\mu =1$ in Equation (23), $\chi =0$ in Equation (11c) and $\beta =0$ in Equation (24b), and $\xi =1$ in Equation (24a), viz.:

$$\begin{array}{cc}\hfill 0& =\sum _{n=0}^{{}^{\infty }}{P}_n,\hfill \end{array}$$

(27a)

$$\begin{array}{cc}\hfill 1& =\sum _{n=0}^{{}^{\infty }}{Q}_n,\hfill \end{array}$$

(27b)

(ii) for the initial velocities Equation (5b), then $\mu =1$ in Equation (23), ${\beta }^{\prime }\equiv \mathrm{d}\beta /\mathrm{d}\mu =-\mathrm{d}x/\mathrm{d}\tau =-u_0/v_{★}$ in Equation (28a) from Equations (23), (24b) and (17a) and also ${\xi }^{\prime }\equiv \mathrm{d}\xi /\mathrm{d}\mu =-\mathrm{d}(1\left(\eta \right)/\mathrm{d}\tau ={\eta }^{\prime }/{\eta }^2=w_0/v_{★}$ using Equations (24a), (23) and (17b):

$$\begin{array}{cc}\hfill -u_0& =v_{★}\sum _{n=0}^{{}^{\infty }}{P}_n(n+p),\hfill \end{array}$$

(28a)

$$\begin{array}{cc}\hfill w_0& =v_{★}\sum _{n=0}^{{}^{\infty }}{Q}_n(n+q).\hfill \end{array}$$

(28b)

If the coefficients $(P_n,Q_n)$ are known for $n=0,1,2,\dots ,N$ then Equations (27a) and (27b), (28a) and (28b) can be used to determine them for $n=N+1,N+2$.

The trajectory is specified by Equations (26a) and (26b), where the coefficients $(P_n,Q_n)$ and the exponents $(p,q)$ are to be determined. The exponents $(p,q)$ can be determined from the lowest powers in the trajectory, Equations (26a) and (26b), viz. Equations (29a)–(29c):

$$\begin{array}{cc}\hfill \beta & \sim P_0{\mu }^p,\hfill \end{array}$$

(29a)

$$\begin{array}{cc}\hfill \xi & \sim Q_0{\mu }^q,\hfill \end{array}$$

(29b)

$$\begin{array}{cc}\hfill H& \equiv {\xi }^{\prime 2}+{\xi }^2{\beta }^{\prime 2},\hfill \end{array}$$

(29c)

which imply that Equation (29c) is given to lowest order by Equation (29d)

$$p⩾1:H={\left(Q_0{\mu }^{q-1}\right)}^2\left[q^2+{\left(P_{0}p{\mu }^p\right)}^2\right]\sim {\left(Q_{0}q{\mu }^{q-1}\right)}^2,$$

(29d)

if $p>0$. Taking also the lowest powers in the equations of motion Equation (25) leads to

$$\mu \left[\begin{array}{c}-P_0{Q}_{0}p(p-1){\mu }^{q+p-2}\\ Q_0^2\left(q{\mu }^{2q-2}+a{\mu }^{2q}\right)\end{array}\right]=Q_0\left\{\frac{{\mu }^{q+1}}{q}\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]-Q_{0}q{\mu }^{2q-1}\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\right\}\left[\begin{array}{c}{P}_{0}p{\mu }^{p-1}\\ Q_{0}q{\mu }^{q-1}\end{array}\right];$$

(30a)

on the left hand side (l.h.s.) ${\mu }^{2q}$ is of higher order than ${\mu }^{2q-2}$ and can be omitted, and on the right hand side (r.h.s.) ${\mu }^{q+1}$ is of higher order than ${\mu }^{2q-1}$ for $q+1>2q-1$ or $q<2$ and can also be omitted, simplifying Equation (30a) to Equation (30b),

$$\mu \left[\begin{array}{c}-P_{0}p(p-1){\mu }^{q+p-2}\\ Q_{0}q{\mu }^{2q-2}\end{array}\right]=-Q_{0}q{\mu }^{2q-1}\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{P}_{0}p{\mu }^{p-1}\\ Q_{0}q{\mu }^{q-1}\end{array}\right].$$

(30b)

This can be satisfied by choosing the exponents of Equation (31a):

$$\begin{array}{cc}\hfill p& =1=q,\hfill \end{array}$$

(31a)

$$\begin{array}{cc}\hfill \left[\begin{array}{c}0\\ Q_0\end{array}\right]& =-Q_0\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{P}_0\\ Q_0\end{array}\right],\hfill \end{array}$$

(31b)

leading to the system of Equation (31b), which can be rewritten as:

$$\begin{array}{cc}\hfill f_{11}{P}_0+f_{12}{Q}_0& =0,\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill f_{21}{P}_0+f_{22}{Q}_0+1& =0,\hfill \end{array}$$

(32b)

and thus solved as:

$$\begin{array}{cc}\hfill P_0& =\frac{f_{12}}{f_{11}{f}_{22}-f_{12}{f}_{21}},\hfill \end{array}$$

(33a)

$$\begin{array}{cc}\hfill Q_0& =-\frac{f_{11}}{f_{11}{f}_{22}-f_{12}{f}_{21}},\hfill \end{array}$$

(33b)

to determine the lowest order coefficients in Equations (26a) and (26b).

The leading terms of the trajectory, Equations (26a) and (26b), are determined by Equations (31a) and (31b), viz.:

$$\begin{array}{cc}\hfill \chi \left(\mu \right)& =\mu (P_0+P_1\mu +\dots ),\hfill \end{array}$$

(34a)

$$\begin{array}{cc}\hfill \xi \left(\mu \right)& =\mu (Q_0+Q_1\mu +\dots ),\hfill \end{array}$$

(34b)

and substituting in the equation of motion (25) up to order ${\mu }^2$ specifies the next coefficients, viz.:

$$\begin{array}{cc}\mu \left[\begin{array}{c}-2P_1{Q}_0\mu \\ 2Q_0{Q}_1\mu +Q_0^2\end{array}\right]=H^{-1/2}{\left(Q_0\mu \right)}^2\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\hfill & \left[\begin{array}{c}{P}_0\\ Q_0\end{array}\right]\hfill \\ & -H^{1/2}\mu (Q_0+Q_1\mu )\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{P}_0+2P_1\mu \\ Q_0+2Q_1\mu \end{array}\right],\hfill \end{array}$$

(35)

which involves Equation (29c), calculated to order $\mu$

$$H\equiv {\left(Q_0+2Q_1\mu \right)}^2+{\left(Q_0{P}_1\mu \right)}^2=Q_0\left(Q_0+4Q_1\mu \right).$$

(36)

Equating to zero the coefficients of $\mu$ in Equation (35), the system of Equations (32a) and (32b) is regained, and equating to zero the coefficients of ${\mu }^2$ leads to the system of equations:

$$\left[\begin{array}{c}-2P_1\\ \phantom{-}2Q_1\end{array}\right]=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\left[\begin{array}{c}{P}_0\\ Q_0\end{array}\right]-2Q_0\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{P}_1\\ Q_1\end{array}\right]-3Q_1\left[\begin{array}{cc}{f}_{11}& f_{12}\\ f_{21}& f_{22}\end{array}\right]\left[\begin{array}{c}{P}_0\\ Q_0\end{array}\right],$$

(37)

which determines the coefficients $(P_1,Q_1)$. The next two pairs of coefficients follow from Equations (27a) and (27b), viz.:

$$\begin{array}{cc}\hfill P_0+P_1+P_2+P_3& =0,\hfill \end{array}$$

(38a)

$$\begin{array}{cc}\hfill Q_0+Q_1+Q_2+Q_3& =1,\hfill \end{array}$$

(38b)

and Equations (28a) and (28b), viz.

$$\begin{array}{cc}\hfill -u_0/v_{★}& =P_0+2P_1+3P_2+4P_3,\hfill \end{array}$$

(39a)

$$\begin{array}{cc}\hfill w_0/v_{★}& =Q_0+2Q_1+3Q_2+4Q_3;\hfill \end{array}$$

(39b)

this specifies the trajectory, Equations (26a) and (26b) up to order four:

$$\begin{array}{cc}\hfill (x-x_0)/\ell & =\sum _{n=0}^{\infty }{P}_n{\left[1-c(t-t_0)/m_0\right]}^{n+1},\hfill \end{array}$$

(40a)

$$\begin{array}{c}exp\left[-(z-z_0)/\ell \right]=\sum _{n=0}^{\infty }{Q}_n{\left[1-c(t-t_0)/m_0\right]}^{n+1}.\end{array}$$

(40b)

In the derivation of Equations (40a) and (40b) from Equations (26a) and (26b), Equations (11c), (23) and (24a) were used and also the exponents in Equation (31a). The first four coefficients in Equations (40a) and (40b) are specified by Equations (33a) and (33b), (37), (38a) and (38b) and (39a) and (39b).

3. General Approach to the Two-Point Boundary-Value Problem (TPBVP)

Any of the preceding three methods can be used to show the effect on a gravity turn of a nonzero angle-of-attack and nonzero angle of the thrust with the rocket axis, including (a) the burn-out altitude and downrange, and (b) vertical and horizontal velocities for: (i) a powered first stage; (ii) any subsequent powered upper stage. The methods I, II and III can also be combined in pairs to address the TPBVP. This method is of considerable interest in general (Section 3.1). A first example of the TPBVP is given (Section 3.2 and Section 3.3) combining methods I and III.

3.1. Smooth Matching of Ascending and Descending Trajectories

The TPBVP is of most interest because it specifies the final burn-out conditions of delivery of the rocket payload. The TPBVP may have no solution if the required burn-out conditions exceed the performance capabilities of the rocket. The method to be presented shows (i) if the TPBVP has a solution in which case it (ii) specifies the whole trajectory and (iii) also the performance envelope for payload delivery at burn-out. In order to discuss these various aspects a general approach to the TPBVP is presented first, before proceeding to its detailed implementation in calculations. The approach to the TPBVP in the present paper is based on two distinct solutions of the rocket trajectory, each specifying downrange x, altitude z, horizontal velocity u and vertical velocity w as a function of time t, viz.:

(i)

an ascending solution for increasing time starting at lift-off up to burn-out

$$t_0⩽t<t_1:X\left(t\right)=\left\{x\left(t\right),z\left(t\right),u\left(t\right),w\left(t\right),x_0,z_0,u_0,w_0\right\};$$

(41a)

(ii)

a descending solution for decreasing time starting at burn-out

$$t_1⩾t>t_0:\overline{X}\left(t\right)=\left\{\overline{x}\left(t\right),\overline{z}\left(t\right),\overline{u}\left(t\right),\overline{w}\left(t\right),x_1,z_1,u_1,w_1\right\}.$$

(41b)

The two solutions must match at cross-over time:

$$n=1,2,3,4:X\left(t_2\right)-\overline{X}\left(t_2\right)=M_n(x_0,z_0,u_0,w_0;x_1,z_1,u_1,w_1;t_2),$$

(42)

leading to a system of four equations $M_n=0$ with $n=1,2,3,4$, which depend on nine parameters, viz. the four lift-off and four burn-out parameters plus the matching time.

The four lift-off conditions are fixed, e.g., launch position at the origin of the coordinate system, Equation (43a) and zero initial velocity, Equation (43b):

$$\begin{array}{cc}\hfill x_0& =0=z_0,\hfill \end{array}$$

(43a)

$$\begin{array}{cc}\hfill u_0& =0=w_0,\hfill \end{array}$$

(43b)

for a first stage; for an upper stage the values would not be zero, but would be fixed anyway, so the four compatibility conditions now depend only on five parameters:

$$0=M_n(0,0,0,0;x_1,z_1,u_1,w_1;t_2)\equiv L_n(x_1,z_1,u_1,w_1;t_2).$$

(44)

Leaving free the matching time, it is possible to choose one burn-out condition, viz. horizontal velocity $u_1$, or vertical velocity $w_1$, or altitude $z_1$ or downrange $x_1$. It is not, in general, possible to choose more than one burn-out condition, e.g., if two were chosen then the three remaining of the five variables would have to satisfy four equations, which is impossible (unless they are redundant equations to start with). The specification of more than one end condition at burn-out might be possible as a control or optimization problem, which lies beyond the scope of the present work. Within the present dynamical method of matching an ascending and a descending trajectory only one end-condition can be imposed, viz. it can be:

$$N(x_1,z_1,u_1,w_1)=0,$$

(45)

any combination of burn-out coordinates or velocities.

The TPBVP can be illustrated (Figure 2) as follows: (i) there is a unique ascending trajectory, Equation (41a), specified by the initial conditions, e.g., Equations (43a) and (43b); (ii) there is a family of descending trajectories, Equation (41b), complying with the burn-out condition Equation (45), e.g., ã given horizontal velocity $u_1$ at burn-out. Thus, two cases can arise: (i) if the ascending trajectory Equation (1) never crosses (Figure 2a) any of the descending trajectories during the burn-time $t_0<t<t_1$, the TPBVP has no solution, because the performance of the rocket is insufficient to achieve the desired burn-out condition; (ii) if the ascending trajectory is tangent (Figure 2b) some of the descending trajectories, then the intersection specifies the matching time $t_0<t_2<t_1$, and the range of possible matching times $t_0\equiv t_{2\min }<t_2<t_{2\max }$ specifies the performance envelope of the rocket at burn-out, e.g., the feasible range horizontal velocity $u_{1\min }⩽u_1⩽u_{1\max }$.

3.2. Feasibility of Desired Burn-Out Condition for Payload Launch

The preceding method (Section 3.1) of addressing the TPBVP covers three issues: (i) existence of solution: does the rocket have enough performance to attain the desired burn-out condition?; (ii) calculation of the trajectory: if the solution exists (i) the method specifies the matching time, so the trajectory can be a calculated from the ascending trajectory, Equation (41a)≡(46a), before and the descending trajectory, Equation (41b)≡(46a) after

$$\{x,z,u,w\}=\left\{\begin{array}{}\mathrm{(46a)}& X\left(t\right),\hfill & \mathrm{if}{t}_0⩽t⩽t_2\hfill \mathrm{(46b)}& \overline{X}\left(t\right),\hfill & \mathrm{if}{t}_2⩽t⩽t_1\hfill \end{array}\right.$$

with assured matching at $t=t_2$; (iii) performance envelope: varying the matching time $t_0⩽t_2⩽t_1$ scans the performance envelope of the rocket in terms of burn-out conditions for payload delivery. The more desirable mathematical procedure to solve the TPBVP is: (i) use the desired burn-out condition Equation (45) with the four matching conditions Equation (42), to eliminate all burn-out variables $(x_1,z_1,u_1,w_1)$ leaving only matching time as a root of Equation (47a):

$$\begin{array}{cc}\hfill 0& =L_0\left(t_2\right),\hfill \end{array}$$

(47a)

$$\begin{array}{cc}\hfill t_0& ⩽t_2⩽t_1;\hfill \end{array}$$

(47b)

(ii) if there is at least one real positive root within burn-out time, Equation (47b) the TPBVP has a solution; (ii) for each root of Equation (47a) specifying a matching time, Equation (47b), the corresponding trajectory, Equation (46a,b) can be calculated; (iii) the conditions that Equation (47a) has at least one positive real root $t_2$ in the range of Equation (47b), specify the performance envelope of the rocket at burn-out.

The present method addresses the: (i) existence of the solution; (ii) the calculation of the trajectory and (iii) the range of final conditions, provided that one critical step can be performed: solving the set of four matching conditions, Equation (42) or Equations (43a), (43b) and (44), which are generally nonlinear, together with the burn-out condition, Equation (45), to obtain an equation, Equation (47a), involving only the matching time. The Equation (47a) is usually a combination of power series of the time arising from the trajectories in two of the methods I, II or III used as ascending or descending solutions; it may involve transcendental functions as well but can usually be truncated as a polynomial in t, viz. of degree N higher for better accuracy. The approximation of Equation (47a) by a polynomial of degree N:

$$0=L_0\left(t_2\right)=\sum _{n=0}^N{L}_n{\left(t_2\right)}^n+O\left(t^{N+1}\right)=L_0\prod _{n=0}^N\left(t_2-t^{\left(n\right)}\right),$$

(48)

will have N roots $t^{\left(n\right)}$, viz.: (i) the complex roots will come in conjugate pairs, because the coefficients $A_n$ are real, and are of no interest; (ii) the real negative roots $t^{\left(n\right)}<0$ are also of no interest; (iii) the real positive roots larger than burn time $t^{\left(n\right)}>t_1$ mean that the propellant is exhausted before reaching the desired final condition; (iv) only the real positive roots within the burn-out range $t_0<t^{\left(n\right)}<t_1$ correspond to feasible solutions of the TPBVP.

The method of approach to the TPBVP described in general before, is implemented next using methods I and II for the calculations. The ascending trajectory is calculated by method I in Section 2.2 which is most accurate for short time, viz. the trajectory is specified as a function of time by: (i) the downrange Equation (49a) from Equations (17a) and (11c):

$$\begin{array}{cc}\hfill x\left(t\right)& =\ell \sum _{n=2}^{\infty }{X}_n{\tau }^n,\hfill \end{array}$$

(49a)

$$\begin{array}{cc}\hfill exp\left[z\left(t\right)/\ell \right]& =1+\sum _{n=2}^{\infty }{Z}_n{\tau }^n,\hfill \end{array}$$

(49b)

and the altitude, Equation (49b), from Equations (17b) and (15a); in both Equation (49a) and (49b) the dimensionless time, Equation (11a) is used, as well as the initial conditions, Equation (43a) and (43b), viz. Equation (50a) is the ascending time:

$$\begin{array}{cc}\hfill \tau & =ct/m_0,\hfill \end{array}$$

(50a)

$$\begin{array}{cc}\hfill \mu & =1-ct/m_0=1-\tau ,\hfill \end{array}$$

(50b)

and Equation (23) ≡ Equation (50b) acts as dimensionless descending time. It appears in the method II in Section 2.3 of trajectory calculation, which is most accurate closer to burn-out, viz. the trajectory is specified for downrange, Equation (51a), by Equations (11c), (24b), (26a) and (31a):

$$\begin{array}{cc}\hfill \overline{x}\left(t\right)& =\ell \sum _{n=0}^{\infty }{P}_n{(1-\tau )}^{n+1},\hfill \end{array}$$

(51a)

$$\begin{array}{cc}\hfill exp\left[-\overline{z}\left(t\right)/\ell \right]& =\sum _{n=0}^{\infty }{Q}_n{(1-\tau )}^{n+1},\hfill \end{array}$$

(51b)

and for the altitude, Equation (51b), by Equations (26b), (31a) and (24a); in the first case the variable is Equation (50a) and in the second the variable is Equation (50b), but both solutions are expressible as power series of $\tau$. The ascending, Equations (49a) and (49b), and descending, Equations (51a) and (51b), solutions are used next to solve the TPBVP in the three aspects of: (i) existence of solution; (ii) calculation of the trajectory; (iii) performance envelope for payload launch.

3.3. Performance Envelope at Burn-Out Condition for Payload Launch

The solution of the TPBVP requires several steps beyond the ascending, Equations (49a) and (49b), and descending, Equations (51a) and (51b), solutions viz.:

(i)

the ascending solutions are needed not only for the downrange, Equation (49a), and altitude, Equation (49b), but also for horizontal, Equation (52a), and vertical, Equation (52b), velocities:

$$\begin{array}{cc}\hfill u\left(t\right)& \equiv \frac{\mathrm{d}x}{\mathrm{d}t}=v_{★}\sum _{n=2}^{\infty }{X}_{n}n{\tau }^{n-1},\hfill \end{array}$$

(52a)

$$\begin{array}{cc}\hfill w\left(t\right)& \equiv \frac{\mathrm{d}z}{\mathrm{d}t}=v_{★}{e}^{-z/\ell }\sum _{n=2}^{\infty }{Z}_{n}n{\tau }^{n-1},\hfill \end{array}$$

(52b)

where $\mathrm{d}\tau /\mathrm{d}t=c/m_0$ was used from Equation (50a), and the reference velocity Equation (18a) was introduced;

(ii)

it also appears in the descending trajectory, Equations (51a) and (51b), when the horizontal, Equation (53a), and vertical, Equation (53b), velocities are calculated:

$$\begin{array}{cc}\hfill \overline{u}\left(t\right)& \equiv \dot{\overline{x}}=-v_{★}\sum _{n=0}^{\infty }{P}_n(n+1){(1-\tau )}^n,\hfill \end{array}$$

(53a)

$$\begin{array}{cc}\hfill \overline{w}\left(t\right)& \equiv \dot{\overline{z}}=-v_{★}{e}^{\overline{z}/\ell }\sum _{n=0}^{\infty }{Q}_n(n+1){(1-\tau )}^n;\hfill \end{array}$$

(53b)

(iii)

the descending solutions, Equations (51a), (51b), (53a) and (53b), can be used to determine the burn-out conditions, when the mass, Equation (2b), equals the residual mass corresponding to

$$m_1\equiv m\left(t_1\right)=m_0-c{t}_1=m_0(1-{\tau }_1);$$

(54)

(iv)

substitution of the dimensionless burn-out time ${\tau }_1$ in the descending solution, Equations (51a) and (51b), (53a) and (53b), specifies four dimensionless constants:

$$\begin{array}{cc}\hfill x_1/\ell & =\sum _{n=0}^{\infty }{P}_n{(m_1/m_0)}^{n+1}\equiv I_1,\hfill \end{array}$$

(55a)

$$\begin{array}{cc}\hfill e^{-z_1/\ell }& =\sum _{n=0}^{\infty }{Q}_n{(m_1/m_0)}^{n+1}\equiv I_2,\hfill \end{array}$$

(55b)

$$\begin{array}{cc}\hfill u_1/v_{★}& =-\sum _{n=0}^{\infty }{P}_n(n+1){(m_1/m_0)}^n\equiv I_3,\hfill \end{array}$$

(55c)

$$\begin{array}{cc}\hfill (w_1/v_{★})e^{-{\overline{z}}_1/\ell }& =-\sum _{n=0}^{\infty }{Q}_n(n+1){(m_1/m_0)}^n\equiv I_4,\hfill \end{array}$$

(55d)

which are calculated from the coefficients of the powers series of time specifying, Equations (51a) and (51b), the descending trajectory.

The next steps are (v) to introduce the burn-out conditions, Equations (55a)–(55d), in the descending trajectory, Equations (51a), (51b), (53a) and (53b), viz.:

$$\begin{array}{cc}\hfill \left[\overline{x}\left(t\right)-x_1\right]/\ell & =-I_1+\sum _{n=0}^{\infty }{P}_n{(1-\tau )}^{n+1},\hfill \end{array}$$

(56a)

$$\begin{array}{cc}\hfill I_{2}exp\left\{-\left[\overline{z}\left(t\right)-z_1\right]/\ell \right\}& =\sum _{n=0}^{\infty }{Q}_n{(1-\tau )}^{n+1},\hfill \end{array}$$

(56b)

$$\begin{array}{cc}\hfill \left[\overline{u}\left(t\right)-u_1\right]/v_{★}& =-I_3-\sum _{n=0}^{\infty }{P}_n(n+1){(1-\tau )}^n,\hfill \end{array}$$

(56c)

$$\begin{array}{cc}\hfill I_4\left[\overline{w}\left(t\right)/v_{★}\right]exp\left\{-\left[\overline{z}\left(t\right)-z_1\right]/\ell \right\}& =-(w_1/v_{★})\sum _{n=0}^{\infty }{Q}_n(n+1){(1-\tau )}^n,\hfill \end{array}$$

(56d)

where subtraction was used to obtain Equations (56a) and (56c) and ratios were used to obtain Equations (56b) and (56d); (vi) step (v) has put the descending solution, Equations (56a)–(56d) in the form of Equation (41b), and the ascending solution, Equations (49a), (49b), (52a) and (52b) is already in the form of Equation (41a), so they can be matched as in Equation (42), viz.:

$$\begin{array}{cc}\hfill \sum _{n=2}^{\infty }{X}_n{\tau }^n+I_1& =x_1/\ell +\sum _{n=0}^{\infty }{P}_n{(1-\tau )}^{n+1},\hfill \end{array}$$

(57a)

$$\begin{array}{cc}\hfill I_{2}exp\left(z_1/\ell \right)& =\left\{\sum _{n=0}^{\infty }{Q}_n{(1-\tau )}^{n+1}\right\}\left\{1+\sum _{n=2}^{\infty }{Z}_n{\tau }^n\right\},\hfill \end{array}$$

(57b)

$$\begin{array}{cc}\hfill \sum _{n=2}^{\infty }{X}_{n}n{\tau }^{n-1}+I_3& =u_1/v_{★}-\sum _{n=0}^{\infty }{P}_n(n+1){(1-\tau )}^n,\hfill \end{array}$$

(57c)

$$\begin{array}{cc}\hfill exp\left(z_1/\ell \right)\times I_4\sum _{n=2}^{\infty }{Z}_n{\tau }^{n-1}& =-(w_1/v_{★})exp[2z\left(t\right)/\ell ]\sum _{n=0}^{\infty }{Q}_n(n+1){(1-\tau )}^n,\hfill \end{array}$$

(57d)

where subtraction was used to derive Equations (57a) and (57c), a product in Equation (57b) and a division in Equation (57d); (vii) all expressions Equations (57a)–(57d) are taken at matching dimensionless time, Equation (58a):

$$\begin{array}{cc}\hfill 0& <{\tau }_2=c{t}_2/m_0<1,\hfill \end{array}$$

(58a)

$$\begin{array}{cc}\hfill 0& <t_2<t_1,\hfill \end{array}$$

(58b)

which must lie within the range of burn time, Equation (58b), otherwise the TPBVP has no solution.

Among the Equations (57a)–(57d) only Equation (57d) involves $z\left(t\right)$ that may (viii) be eliminated substituting Equation (49b) leading to Equation (59a)

$$exp\left(z_1/\ell \right)\times I_4\sum _{n=2}^{\infty }{Z}_n{\tau }^{n-1}=-(w_1/v_{★}){\left\{1+\sum _{n=2}^{\infty }{Z}_n{\tau }^n\right\}}^2\left\{\sum _{n=0}^{\infty }{Q}_n(n+1){(1-\tau )}^n\right\};$$

(59a)

further substitution of Equation (57b) in Equation (59a) eliminates $exp(z_1/\ell )$ leaving only $w_1$ in Equation (59b)

$$I_4\left\{\sum _{n=2}^{\infty }{Z}_n{\tau }^{n-1}\right\}\left\{\sum _{n=0}^{\infty }{Q}_n{(1-\tau )}^{n+1}\right\}=-I_2(w_1/v_{★})\left\{1+\sum _{n=2}^{\infty }{Z}_n{\tau }^n\right\}\left\{\sum _{n=0}^{\infty }{Q}_n(n+1){(1-\tau )}^n\right\}.$$

(59b)

The matching conditions are interpreted as follows at the next steps: (ix) if the desired downrange $x_1$ at the burn-out time is specified, then the matching time ${\tau }_2$ is a root of Equation (57a), if the burn-out altitude $z_1$ is specified then ${\tau }_2$ is a root of Equation (57b), if the horizontal velocity $u_1$ is specified it is a root of Equation (57c) and if the vertical velocity $w_1$ is specified at the burn-out then the dimensionless matching time ${\tau }_2$ is the root of Equation (59b). The derivation of Equations (56a)–(56d), (57a)–(57c) and (59b) is an example of the procedure to solve a nonlinear system of four equations obtained by matching Equations (56a)–(56d) to Equations (49a), (49b), (51a) and (51b). Equation (59b) is the most complex equation among the matching equations, when the vertical velocity at the burn-out is specified, and is used next to illustrate the remaining steps of solution of the TPBVP.

The remaining steps in the TPBVP solution are: (x) the equation of a single series on the l.h.s. to the product of three series on the r.h.s. of Equation (59b), can be truncated at orders $m,n,r=1,\dots ,N$ leading to a polynomial of degree N, e.g., for $N=2$ only powers up to two are retained:

$$\begin{array}{cc}\hfill I_4\tau \left(Z_2+Z_3\tau \right)& \left(Q_0(1-\tau )+Q_1{(1-\tau )}^2\right)\hfill \\ \hfill & =-(w_1/v_{★})I_2\left(1+Z_2{\tau }^2\right)\left[Q_0+2Q_1(1-\tau )+3Q_2{(1-\tau )}^2\right]\hfill \end{array}$$

(60)

(xi) the preceding result is a polynomial of the second-degree, Equation (61a),

$$\begin{array}{cc}\hfill A_0& +A_1\tau +A_2{\tau }^2=0,\hfill \end{array}$$

(61a)

$$\begin{array}{cc}\hfill A_0& \equiv (w_1/v_{★})I_2(Q_0+2Q_1+3Q_2),\hfill \end{array}$$

(61b)

$$\begin{array}{cc}\hfill A_1& \equiv I_4(Q_0+Q_1)Z_2-(w_1/v_{★})I_2(2Q_1+6Q_2),\hfill \end{array}$$

(61c)

$$\begin{array}{cc}\hfill A_2& \equiv I_4\left[Z_3(Q_0+Q_1)-Z_2(Q_0+2Q_1)\right]+(w_1/v_{★})I_2\left[Z_2(Q_0+2Q_1+3Q_2)+3Q_2\right],\hfill \end{array}$$

(61d)

with coefficients in Equations (61b)–(61d); (xii) the matching time is specified by the roots

$$0<t_2=\frac{-A_1\pm \sqrt{A_1^2-4A_0{A}_2}}{2A_2}<1-m_1/m_0,$$

(62)

if they are real and lie in the range of burn time; (xiii) the positive real root(s) of Equation (62) specifies ${\tau }_2$ and $t_2=m_0{\tau }_2/c$, so the initial trajectory can be calculated from Equations (49a), (49b), (52a) and (52b) from lift-off to matching time $0⩽t⩽t_2$ and by Equations (51a), (51b), (53a) and (53b) or (56a)–(56d) after matching time till burn-out $t_2⩽t⩽t_1$; (xiv) the performance envelope of the rocket at burn-out in terms of the vertical velocity $w_1$ it can achieve is specified by substituting Equations (61b)–(61d) in the inequalities present in Equation (62). The truncation could be made at a higher degree N for better accuracy.

4. Trajectory for a Given Horizontal Velocity at Burn-Out

One criterion of major interest at burn-out is the horizontal velocity, which will be calculated applying the general method (Section 4) to find the matching time (Section 4.1) and which is applied to specific rocket data (Section 4.2) to assess the application (Section 4.3) to the TPBVP.

4.1. Matching Time as a Root of the Series Solution

If the horizontal velocity $u_1$ at the burn-out is given the matching time ${\tau }_2$ between ascending and descending trajectories is a root of Equation (57c), viz.:

$$N=1:u_1/v_{★}=I_3+\sum _{n=0}^{\infty }{P}_n(n+1){(1-\tau )}^n+\sum _{n=2}^{\infty }{X}_{n}n{\tau }^{n-1}.$$

(63)

The series on the r.h.s. may be truncated at order N leading to a polynomial of degree N, e.g.,: (i) for $N=1$,

$$u_1/v_{★}-I_3=P_0+2P_1(1-\tau )+2X_2\tau =P_0+2P_1+2(X_2-P_1)\tau .$$

(64)

The approximation at the same order of power series with different variables $\tau$ and $1-\tau$, like the passage from from Equation (59b) to Equation (60), or the passage from Equation (63) to Equation (64) is discussed in the conclusion (Section 7). The matching time is

$$0⩽{\tau }_2=\frac{P_0+2P_1+I_3-u_1/v_{★}}{2(P_1-X_2)}⩽1-m_1/m_0,$$

(65)

and the inequalities specify the minimum and maximum horizontal velocities achievable at the burn-out that are the extreme values

$$\begin{array}{cc}\hfill u_{11}& =v_{★}(P_0+2P_1+I_3)\hfill \end{array}$$

(66a)

$$\begin{array}{cc}\hfill u_{12}& =v_{★}\left[P_0+2P_1+I_3-2(P_1-X_2)(1-m_1/m_0)\right].\hfill \end{array}$$

(66b)

The maximum (minimum) among the two values specified in Equations (66a) and (66b) depends on which is larger (smaller) for the particular values of $P_0,P_1,X_2,I_3,m_1/m_0$. The order $N=1$ is the lowest possible, and more accurate results would be obtained for $N=2$ and beyond. However, this case is sufficient for illustration with minimal computation. A feature of the approximation $N=1$ is that there is always a real solution, but it is unknown whether it lies in the burn time range, Equation (65).

If the approximation $N=2$ is made in Equation (63),

$$N=2:u_1/v_{★}-I_3=P_0+2P_1(1-\tau )+3P_2{(1-\tau )}^2+2X_2\tau +3X_3{\tau }^2,$$

(67)

then a polynomial of second degree, Equation (61a), is obtained with coefficients

$$\left\{A_0,A_1,A_2\right\}=\left\{I_3+P_0+2P_1+3P_2-u_1/v_{★},2X_2-6P_2-2P_1,3(P_2+X_3)\right\}.$$

(68)

The roots of Equation (62) will be complex unless the condition of Equation (69a) is met:

$$\begin{array}{cc}\hfill A_1^2& ⩾4A_0{A}_2,\hfill \end{array}$$

(69a)

$$\begin{array}{cc}\hfill 2A_2{\tau }_2& =-A_1+\sqrt{A_1^2-4A_0{A}_2},\hfill \end{array}$$

(69b)

and then only one root, Equation (69b), might be positive, so there is at most one matching time. It is clear that the root of Equation (69b) will be positive if $A_0{A}_2<0$ and then it must lie in the burn-out range Equation (65). The condition $A_0{A}_2<0$ implies that Equation (69a) is satisfied, so a necessary condition for the TPBVP to have a solution is

$$0>A_0{A}_2=3(P_2+X_3)(I_3+P_0+2P_1+3P_2-u_1/v_{★}).$$

(70)

If the first factor is negative

$$P_2+X_3<0:u_1<v_{★}(I_3+P_0+2P_1+3P_2)\equiv u_{1\max },$$

(71)

the horizontal velocity at the burn-out cannot exceed Equation (71).

4.2. Rocket Data Required for Trajectory Calculation

The initial data on environment, propulsion and aerodynamics in

Table 1. Input data for the calculation of rocket trajectories, extended from vertical ascent (Table 2 in [21]) to a gravity turn. The matrices were calculated for $\alpha =2^{\circ }$ and $ϵ=0$ (cf. text)

Symbol	Meaning	Value	Unit
Environment
g	acceleration of gravity	9.81	$\mathrm{m}/\mathrm{s}$
${\rho }_0$	sea level mass density	1.225	$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
ℓ	atmospheric scale height	$2.6\times {10}^4$	$\mathrm{m}$
Propulsion
T	thrust	$1.555\times {10}^7$	$\mathrm{k}\mathrm{g}\mathrm{m}/{\mathrm{s}}^2$
c	propellant flow rate	$2.029\times {10}^3$	$\mathrm{k}\mathrm{g}/\mathrm{s}$
$m_0-m_1$	propellant mass	$2.841\times {10}^5$	$\mathrm{k}\mathrm{g}$
$t_1-t_0$	burn time	$1.4\times {10}^2$	$\mathrm{s}$
Aerodynamics
$m_0$	initial mass	$7.77\times {10}^5$	$\mathrm{k}\mathrm{g}$
S	cross-sectional area	3.76 x 10	${\mathrm{m}}^2$
$C_{\mathrm{D}}$	drag coefficient	0.15	*
$C_{\mathrm{L}}$	lift coefficient	0.10	*
Calculated parameters
$m_{*}$	reference mass	$1.198\times {10}^6$	$\mathrm{k}\mathrm{g}$
a	weight parameter	5.533 × 10	*
b	thrust parameter	$1.129\times {10}^2$	*
f	aerodynamic parameter	$7.706\times {10}^{-1}$	*
$t_{*}-t_0$	reference time	$3.829\times {10}^2$	$\mathrm{s}$
$v_{*}$	reference velocity	6.789 × 10	$\mathrm{m}/\mathrm{s}$
Calculated matrices
trust		aerodynamic
$b_{11}$	1.128 × 102	$f_{11}$	0.118
$b_{12}$	−3.94	$f_{12}$	0.073
$b_{21}$	$3.94$	$f_{21}$	$-0.073$
$b_{22}$	1.128102	$f_{22}$	0.118

* Dimensionless parameter.

is the same as in

Table 2. Application of distinct pairs of methods to TPBVP.

Pair of Methods *
I + II (Section 4)	III + II (Section 5)
$X_2=$ 3.85 × 101	$R_0=0=R_1$
$P_0=3.78$	$S_0=0=S_1$
$Q_0=-6.12$	$R_2=6.82 \mathrm{m}/{\mathrm{s}}^2$
$P_1=$ −1.80 × 102	$S_2=2.41 \mathrm{m}/{\mathrm{s}}^2$
$Q_1=$ 2.17 × 102	$R_3=$ 5.94 × 10−3 $\mathrm{m}/{\mathrm{s}}^3$
$I_3=$ 7.88 × 101	$S_3=$ 6.37 × 10−3 $\mathrm{m}/{\mathrm{s}}^3$
$u_{11}=$ −1.88 × 104 $\mathrm{m}/\mathrm{s}$	$u_{13}=$ −1.88 × 104 $\mathrm{m}/\mathrm{s}$
$u_{12}=$ −7.98 × 103 $\mathrm{m}/\mathrm{s}$	$u_{14}=$ −7.98 × 103 $\mathrm{m}/\mathrm{s}$

* In all cases: $u_1=$ 5.35 × 10³ $\mathrm{m}/\mathrm{s},m_1/m_0=0.634,{\gamma }_0=\pi /4$.

in [21]; the calculated parameters for a one-dimensional trajectory (i.e., vertical climb) are extended in the case of a two-dimensional trajectory (i.e., a gravity turn), in the case of thrust, Equation (13a), and aerodynamic force, Equation (14a), given by matrices, involving: (i) for thrust the angle-attack $\alpha =2^{\circ }$ and thrust vector angle $\epsilon =0^{\circ }$; (ii) for aerodynamic forces only angle-of-attack, drag $C_{\mathrm{D}}$ and lift $C_{\mathrm{L}}$ coefficients. The thrust vector angle was chosen to be $\epsilon =0^{\circ }$ because it is often used for trim only, so it varies over longer time scales, i.e., over short time scales of interest there is a “static” balance; the angle-of-attack takes a small value $\alpha =2^{\circ }$ since at high speeds the structural loads would become too large otherwise.

The data in the Table 1 consists of five sets: (i) acceleration due to gravity g, mass density $\rho$ and scale height ℓ for the International Standard Atmosphere (ISA) at sea level; (ii) thrust T, propellant mass $m_0-m_1$ and burn time $t_1-t_0=(m_0-m_1)/c$ for the Vulcain first stage of the Ariane 5 rocket [30,31]; (iii) lift-off mass $m_0$, cross-sectional area S, with assumed values for the drag $C_{\mathrm{D}}$ and lift $C_{\mathrm{L}}$ coefficients; (iv) the preceding data is used to calculate the reference mass Equation (14b), the weight parameter Equation (12a), the thrust parameter Equation (13b), the aerodynamic parameter in Equation (14c), reference time corresponding to the burn-out time if all mass was propellant $t_{★}-t_0=m_0/c>t_1-t_0$ and reference velocity Equation (18a); (v) the angle of attack $\alpha =2^{\circ }$ is used with thrust parameter b to calculate the components of the dimensionless thrust matrix Equation (13a) with thrust axis $ϵ=0$, and the aerodynamic matrix Equation (14a) involving the aerodynamic parameter f, Equation (14c), and drag $C_{\mathrm{D}}$ and lift $C_{\mathrm{L}}$ coefficients. This is all the data needed for subsequent calculation of trajectories specified in Table 2 and

Table 3. Solution of TPBVP combining methods I and III.

Expression	Value	Unit
$u_1$	$5.35\times {10}^3$	$\mathrm{m}/\mathrm{s}$
$u_1/v_{★}=J_3$	$7.88\times {10}^1$
$1-m_1/m_0$	0.366
$2X_2(1-m_1/m_0)$	$2.81\times {10}^1$
$u_{15}=v_{★}{J}_3$	$5.35\times {10}^3$	$\mathrm{m}/\mathrm{s}$
$2X_2{v}_{★}$	$5.23\times {10}^3$	$\mathrm{m}/\mathrm{s}$
$2R_2{m}_0/c$	$5.23\times {10}^3$	$\mathrm{m}/\mathrm{s}$
$u_{16}$	$5.35\times {10}^3$	$\mathrm{m}/\mathrm{s}$

.

4.3. Combination of Methods I and II for the TPBVP

The lowest order approximation to matching time, Equation (65), uses from the ascending solution only the coefficient (72) in Equation (21a)

$$2X_2=b_{11}cos{\gamma }_0+b_{12}sin{\gamma }_0,$$

(72)

for the initial conditions Equation (43b) with initial flight path angle ${\gamma }_0$, implying that the last term on the r.h.s. of Equation (21a) vanishes but not the first two terms; choosing launch at an oblique angle Equation (73a) leads to Equation (73b)

$$\begin{array}{cc}\hfill {\gamma }_0& =\pi /4:\hfill \end{array}$$

(73a)

$$\begin{array}{cc}\hfill X_2& =\frac{b_{11}+b_{12}}{2\sqrt{2}},\hfill \end{array}$$

(73b)

in Table 2, in the first column concerning matching of methods I and II. Concerning the descending trajectory, the lowest order approximation of Equation (65) to matching time involves only $P_0$ in Equation (33a) and $P_1$ is calculated from Equation (37), viz.:

$$\left[\begin{array}{cc}-2+2Q_0{f}_{11}& 2Q_0{f}_{12}\\ 2Q_0{f}_{21}& 2+2Q_0{f}_{22}\end{array}\right]\left[\begin{array}{c}{P}_1\\ Q_1\end{array}\right]=\left[\begin{array}{cc}{b}_{11}-3Q_1{f}_{11}& b_{12}-3Q_1{f}_{12}\\ b_{21}-3Q_1{f}_{21}& b_{22}-3Q_1{f}_{22}\end{array}\right]\left[\begin{array}{c}{P}_0\\ Q_0\end{array}\right].$$

(74)

Due to the small angle of attack $\alpha =2^{\circ }$, in the thrust matrix Equation (13a) the cross-terms $b_{12}$ and $b_{21}$ are small relative to the diagonal terms $b_{11}$ and $b_{22}$, by about two orders of magnitude as seen in Table 1, which also shows that the terms of the aerodynamic matrix are all comparable in magnitude. Thus, all terms in Equations (33a) and (33b) are retained. Equation (74) is a linear system in $P_1,Q_1$ whose solution is a linear combination of the thrust matriz terms $b_{ij}$ and therefore, we can neglect the cross-terms to obtain the simplified solution as specified in Equations (75a) and (75b):

$$\begin{array}{cc}\hfill P_1& =\frac{f_{12}\left(b_{11}(f_{12}{f}_{21}-3f_{11}{f}_{22})-2b_{22}{f}_{11}^2\right)}{2(f_{11}{f}_{22}-f_{12}{f}_{21})\left(3f_{11}^2-f_{12}{f}_{21}+3f_{11}{f}_{22}\right)},\hfill \end{array}$$

(75a)

$$\begin{array}{cc}\hfill Q_1& =\frac{f_{11}\left(b_{11}{f}_{12}{f}_{21}+b_{22}(f_{11}^2-f_{12}{f}_{21}+f_{11}{f}_{22})\right)}{(f_{11}{f}_{22}-f_{12}{f}_{21})\left(3f_{11}^2-f_{12}{f}_{21}+3f_{11}{f}_{22}\right)}\hfill \end{array}$$

(75b)

Thus, $P_0$ and $P_1$ in Table 2 are calculated from Equations (33a) and (75a), respectively. The burn-out velocity is assumed to be $u_1=5.35 \mathrm{k}\mathrm{m}/\mathrm{s}$, leading Equation (55c) to $I_3=u_1/v_{★}$ in Table 2. This specifies Equation (66a) the value $u_{11}$ in the Table 2 and also from Equation (66b) the value of $u_{12}$. Since both are negative, a positive horizontal velocity between them is not possible, and the combination of methods I and II fails to solve the TPBVP. Thus, a second alternative is considered for TPBVP using another pair of methods.

5. Second Alternative Set of Matching Conditions for the TPBVP

The second solution (Section 5) to the TPBVP (Section 3) differs from the first (Section 4) in using method III (Section 2.1) for the ascending solution instead of method I (Section 2.2), while retaining method II (Section 2.3) for the descending solution.

5.1. Alternative Choices of Ascent Trajectories for Matching up to Burn-Out

Method II in Section 2.3 is retained in Equations (56a)–(56d) for the descending solution; method I in Section 2.2 for the ascending solution from Equations (49a), (49b), (52a) and (52b) is replaced by method III in Section 2.1, viz.:

$$\begin{array}{cc}\hfill \overline{\overline{x}}\left(t\right)& =\sum _{n=0}^{\infty }{R}_n{(t-t_0)}^n,\hfill \end{array}$$

(76a)

$$\begin{array}{cc}\hfill \overline{\overline{z}}\left(t\right)& =\sum _{n=0}^{\infty }{S}_n{(t-t_0)}^n,\hfill \end{array}$$

(76b)

where the coefficients:

$$\begin{array}{cc}\hfill R_n& \equiv x^{\left(n\right)}\left(t_0\right)/n!,\hfill \end{array}$$

(77a)

$$\begin{array}{cc}\hfill S_n& \equiv z^{\left(n\right)}\left(t_0\right)/n!,\hfill \end{array}$$

(77b)

apply as well to the velocities:

$$\begin{array}{cc}\hfill \overline{\overline{u}}\left(t\right)\equiv \dot{\overline{\overline{x}}}& =\sum _{n=1}^{\infty }{R}_{n}n{(t-t_0)}^{n-1},\hfill \end{array}$$

(78a)

$$\begin{array}{cc}\hfill \overline{\overline{w}}\left(t\right)\equiv \dot{\overline{\overline{z}}}& =\sum _{n=1}^{\infty }{S}_{n}n{(t-t_0)}^{n-1}.\hfill \end{array}$$

(78b)

The first four coefficients follow from the initial conditions, Equations (43a) and (43b) applied to Equation (7):

$$\begin{array}{cc}\hfill R_1=R_0& =0=S_0=S_1:\hfill \end{array}$$

(79a)

$$\begin{array}{cc}\hfill m_0\left[\begin{array}{c}2R_2\\ 2S_2+g\end{array}\right]& =\left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]\left[\begin{array}{c}{u}_0/v_0\\ w_0/v_0\end{array}\right],\hfill \end{array}$$

(79b)

and to Equation (9),

$$\begin{array}{cc}3m_0\left[\begin{array}{c}2R_3\\ 2S_3\end{array}\right]=c\left[\begin{array}{c}2R_2\\ 2S_2+g\end{array}\right]+\frac{2}{v_0^3}\hfill & (w_0{R}_2-u_0{S}_2)\left[\begin{array}{cc}{T}_{11}& T_{12}\\ T_{21}& T_{22}\end{array}\right]\left[\begin{array}{c}\phantom{-}{w}_0\\ -u_0\end{array}\right]\hfill \\ & -\frac{1}{2}\frac{S{\rho }_0}{v_0}\left[\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right]\left[\begin{array}{c}2R_2(2u_0^2+w_0^2)+2u_0{w}_0{S}_2-v_0^2{u}_0{w}_0/\ell \\ 2S_2(u_0^2+2w_0^2)+2u_0{w}_0{R}_2-v_0^2{w}_0^2/\ell \end{array}\right],\hfill \end{array}$$

(80)

these Equations (79a), (79b) and (80) specify the first four terms of the trajectory, Equations (76a) and (76b).

Although the initial horizontal and vertical velocities are zero, they are related to the initial total velocity, Equation (6a), and launch angle, Equation (6b), by:

$$\begin{array}{cc}\hfill u_0& =v_{0}cos{\gamma }_0,\hfill \end{array}$$

(81a)

$$\begin{array}{cc}\hfill w_0& =v_{0}sin{\gamma }_0,\hfill \end{array}$$

(81b)

and the launch angle ${\gamma }_0$ can be chosen freely, viz.:

(i)

in Equations (79a) and (79b):

$$\left[\begin{array}{c}2R_2\\ 2S_2+g\end{array}\right]=\frac{T}{m_0}\left[\begin{array}{cc}cos(\alpha +\epsilon )& -sin(\alpha +\epsilon )\\ sin(\alpha +\epsilon )& \phantom{-}cos(\alpha +\epsilon )\end{array}\right]\left[\begin{array}{c}cos{\gamma }_0\\ sin{\gamma }_0\end{array}\right],$$

(82)

where Equation (13a) was used in terms of the angle-of-attack $\alpha$ and thrust angle $\epsilon$;

(ii)

in Equations (80) and (79b),

$$\left[\begin{array}{c}2R_3\\ 2S_3\end{array}\right]=\frac{c}{3m_0}\left[\begin{array}{c}2R_2\\ 2S_2+g\end{array}\right]=\frac{cT}{3m_0^2}\left[\begin{array}{cc}cos(\alpha +\epsilon )& -sin(\alpha +\epsilon )\\ sin(\alpha +\epsilon )& \phantom{-}cos(\alpha +\epsilon )\end{array}\right]\left[\begin{array}{c}cos{\gamma }_0\\ sin{\gamma }_0\end{array}\right].$$

(83)

The last term is $O\left(v_0\right)$ and tends to zero as $v_0\to 0$, and the second term which is singular $O(1/v_0)$ as $v_0\to 0$ is omitted, by excluding that thrust causes a discontinuous third-order time derivative of position $\stackrel{⃛}{x}=\infty =\stackrel{⃛}{z}$ (the result ${\stackrel{⃛}{x}}_0=\infty ={\stackrel{⃛}{z}}_0$ is a feature of “instantaneous” thrust application at time $t=0 \mathrm{m}/\mathrm{s}$, but a real thrust application is gradual).

5.2. Matching Distinct Ascending Solutions to the Same Descending Solution

The second solution to the TPBVP is obtained matching the descending solution, Equations (56a)–(56d), from method II:

$$\begin{array}{cc}\hfill \overline{x}\left(t_2\right)& =\overline{\overline{x}}\left(t_2\right),\hfill \end{array}$$

(84a)

$$\begin{array}{cc}\hfill \overline{z}\left(t_2\right)& =\overline{\overline{z}}\left(t_2\right),\hfill \end{array}$$

(84b)

$$\begin{array}{cc}\hfill \overline{u}\left(t_2\right)& =\overline{\overline{u}}\left(t_2\right),\hfill \end{array}$$

(84c)

$$\begin{array}{cc}\hfill \overline{w}\left(t_2\right)& =\overline{\overline{w}}\left(t_2\right),\hfill \end{array}$$

(84d)

to the ascending solution, Equations (76a), (76b), (78a) and (78b) from method III:

$$x_1/\ell -I_1+\sum _{n=0}^{\infty }{P}_n{(1-\tau )}^{n+1}=(1/\ell )\sum _{n=0}^{\infty }{R}_n{(m_0\tau /c)}^n,$$

(85a)

$$exp\left(-z_1/\ell \right)\sum _{n=0}^{\infty }{Q}_n{(1-\tau )}^{n+1}=I_{2}exp\left[-(1/\ell )\sum _{n=0}^{\infty }{S}_n{(m_0\tau /c)}^n\right],$$

(85b)

$$u_1/v_{★}-I_3-\sum _{n=0}^{\infty }{P}_n(n+1){(1-\tau )}^n=(1/v_{★})\sum _{n=1}^{\infty }{R}_{n}n{(m_0\tau /c)}^{n-1},$$

(85c)

$$w_{1}exp\left[(\overline{z}-z_1)/\ell \right]\sum _{n=0}^{\infty }{Q}_n(n+1){(1-\tau )}^n=-I_4\sum _{n=1}^{\infty }{S}_{n}n{(m_0\tau /c)}^{n-1},$$

(85d)

where the relation with dimensionless time $\tau =ct/m_0$ and the reference velocity Equation (18a) were used. If the burn-out downrange $x_1$ is given the root of Equation (85a) specifies the matching time $\tau ={\tau }_2$ and likewise Equation (85b) is used if the burn-out altitude $z_1$ is specified, and Equation (85c) if the horizontal velocity at burn-out $u_1$ is specified.

If the vertical velocity at burn-out $w_1$ is specified, then the matching time ${\tau }_2$ is the root $\tau ={\tau }_2$ of Equation (85d), where $exp[(\overline{z}-z_1)/\ell ]$ is substituted from Equation (56b)

$$I_2{w}_1\sum _{n=0}^{\infty }{Q}_n(n+1){(1-\tau )}^n=-I_4\left\{\sum _{n=0}^{\infty }{Q}_n{(1-\tau )}^{n+1}\right\}\left\{\sum _{k=1}^{\infty }{S}_{k}k{(m_0\tau /c)}^{k-1}\right\}.$$

(86)

Selecting the horizontal velocity $u_1$ at the burn-out and taking Equation (85c) to the lowest order, viz. the first order, leads to

$$N=1:u_1/v_{★}-I_3-P_0-2P_1(1-\tau )=[R_1+2R_2(m_0\tau /c)]/v_{★}.$$

(87)

From Equation (87), noting Equation (79a), it follows that the matching time is

$$0⩽{\tau }_2=\frac{u_1-v_{★}(I_3+P_0+2P_1)}{2R_2{m}_0/c-2v_{★}{P}_1}⩽1-m_1/m_0.$$

(88)

Thus, the range of horizontal velocities at the burn-out lies between the extreme values

$$\begin{array}{cc}\hfill u_{13}& =v_{★}(I_3+P_0+2P_1),\hfill \end{array}$$

(89a)

$$\begin{array}{cc}\hfill u_{14}& =u_{13}+(1-m_1/m_0)(2R_2{m}_0/c-2v_{★}{P}_1).\hfill \end{array}$$

(89b)

If the term in the second brackets in Equation (89b) is positive, Equation (89c), then Equation (89b) is the highest and Equation (89a) is the lowest limit:

$$\begin{array}{cc}\hfill 2R_2{m}_0& >2v_{★}{P}_{1}c,\hfill \end{array}$$

(89c)

$$\begin{array}{cc}\hfill u_{13}& ⩽u_1⩽u_{14}.\hfill \end{array}$$

(89d)

If the inequality Equation (89c) is reversed then the signs in Equation (89d) are also reversed. The second solution (Section 5.1 and Section 5.2) of the TPBVP will be tested next (Section 5.3) with the same data as the first TPBVP in Section 4.

5.3. Trajectory Matching outside the Burn Range

The extreme values $u_{11}$ in Equation (66a), for the first solution of the TPBVP, and $u_{13}$ in Equation (89a), for the second solution of the TPBVP, coincide $u_{12}=u_{13}$. The other extreme value $u_{14}$ in Equation (89b) specified by the combination of methods III+II is distinct $u_{14}\ne u_{11}$ from the extreme value $u_{11}$ in Equation (66b) specified by the combination of methods I+II and is calculated next. This requires the coefficient $R_2$, that appears together with $S_2,R_3,S_3$ in the two lowest orders for the ascending solution, Equations (76a), (76b), (77a) and (77b):

$$\begin{array}{cc}\hfill \overline{\overline{x}}\left(t\right)& =R_2{t}^2+R_3{t}^3,\hfill \end{array}$$

(90a)

$$\begin{array}{cc}\hfill \overline{\overline{z}}\left(t\right)& =S_2{t}^2+S_3{t}^3,\hfill \end{array}$$

(90b)

$$\begin{array}{cc}\hfill \overline{\overline{u}}\left(t\right)& =2R_{2}t+3R_3{t}^2,\hfill \end{array}$$

(90c)

$$\begin{array}{cc}\hfill \overline{\overline{w}}\left(t\right)& =2S_{2}t+3S_3{t}^2,\hfill \end{array}$$

(90d)

has two lowest order nonzero coefficients specified, Equations (82) and (83), by the vector

$$\left[\begin{array}{cc}cos(\alpha +\epsilon )& -sin(\alpha +\epsilon )\\ sin(\alpha +\epsilon )& \phantom{-}cos(\alpha +\epsilon )\end{array}\right]\left[\begin{array}{c}cos{\gamma }_0\\ sin{\gamma }_0\end{array}\right]=\left[\begin{array}{c}0.682\\ 0.731\end{array}\right],$$

(91)

which is calculated assuming: (i) the thrust is aligned with the rocket axis $\epsilon =0$, for low-frequency trim; (ii) the angle-of-attack is small $\alpha =2^{\circ }$ to limit aerodynamic forces; (iii) the launch angle is Equation (73a). These values can be used in Equations (82) and (83)

$$\frac{m_0}{T}\left[\begin{array}{c}2R_2\\ 2S_2+g\end{array}\right]=\left[\begin{array}{c}0.682\\ 0.731\end{array}\right]=\frac{3m_0^2}{cT}\left[\begin{array}{c}2R_3\\ 2S_3\end{array}\right],$$

(92)

to calculate the coefficients $R_2,S_2,R_3,S_3$ in Table 2. The second extreme velocity Equation (89b) specified by the combination of methods III and II gives the same result in Table 2 as the combination of methods I and II although Equation (66a) uses different parameters. Thus, the method II combined with either the method I in Section 4 or method III in Section 5, leads to the same negative extreme value $u_{11}=u_{13}$ and $u_{12}=u_{14}$ in Table 2 that are both negative and fail to solve the TPBVP. This suggests excluding method II and combining methods I and III for the solution (Section 6) of the TPBVP.

6. Six Trajectories of the TPBVP Using Three Pairs of Solutions

The three methods I, II, III of trajectory calculation lead to three combinations for the TPBVP: (Section 4) matching ascending I to descending II; (Section 5) matching ascending III to descending II; (Section 6) matching III and I, e.g., ascending III and descending I. The number of trajectories solving the TPBVP would double from three to six interchanging the ascending and descending solutions.

6.1. Matching of the Third Pair of Solutions in Two Forms

The TPBVP has been addressed: (i-ii) in Section 4 using I + II method I for the ascending solution and method II for the descending solution (the reverse II + I would also be possible); (iii-iv) in Section 5 using III + II method III for the ascending solution (the reverse II + III would also be possible); (v-vi) the remaining combination is in Section 6 is III+I descending I and ascending III (or the reverse I + III). The method III is already in an ascending form, Equations (76a), (76b), (78a) and (78b). The method I is also in an ascending form Equations (49a), (49b), (52a) and (52b), but it can be put into a descending form by: (i) using the ascending form to calculate burn-out conditions for the downrange and altitude:

$$\begin{array}{cc}\hfill x_1/\ell & =\sum _{n=2}^{\infty }{X}_n{(1-m_1/m_0)}^n\equiv J_1,\hfill \end{array}$$

(93a)

$$\begin{array}{cc}\hfill exp\left(z_1/\ell \right)& =1+\sum _{n=2}^{\infty }{Z}_n{(1-m_1/m_0)}^n\equiv J_2,\hfill \end{array}$$

(93b)

and the horizontal and vertical velocities:

$$\begin{array}{cc}\hfill u_1/v_{★}& =\sum _{n=2}^{\infty }{X}_{n}n{(1-m_1/m_0)}^{n-1}\equiv J_3,\hfill \end{array}$$

(93c)

$$(w_1/v_{★})exp\left(z_1/\ell \right)=\sum _{n=2}^{\infty }{Z}_{n}n{(1-m_1/m_0)}^{n-1}\equiv J_4,$$

(93d)

bearing in mind that the dimensionless time at burn-out ${\tau }_1=1-m_1/m_0$ is unity minus the residual mass fraction.

6.2. Third Set of Matching Condition for TPBVP Trajectories

The next step (ii) is to subtract the burn-out conditions, Equations (93a)–(93d), from the ascending solution, Equations (49a), (49b), (52a) and (52b) to obtain the descending solution of method I:

$$\begin{array}{cc}\hfill \left[x\left(t\right)-x_1\right]/\ell & =-J_1+\sum _{n=2}^{\infty }{X}_n{\tau }^n,\hfill \end{array}$$

(94a)

$$\begin{array}{cc}\hfill J_{2}exp\left\{\left[z\left(t\right)-z_1\right]/\ell \right\}& =1+\sum _{n=2}^{\infty }{Z}_n{\tau }^n,\hfill \end{array}$$

(94b)

$$\begin{array}{cc}\hfill \left[u\left(t\right)-u_1\right]/v_{★}& =-J_3+\sum _{n=2}^{\infty }{X}_{n}n{\tau }^{n-1},\hfill \end{array}$$

(94c)

$$\begin{array}{cc}\hfill J_4\left[w\left(t\right)/w_1\right]exp\left\{\left[z\left(t\right)-z_1\right]/\ell \right\}& =\sum _{n=2}^{\infty }{Z}_{n}n{\tau }^{n-1}.\hfill \end{array}$$

(94d)

Equating the descending solution, Equations (94a)–(94d), by method I to the ascending solution, Equations (76a), (76b), (78a) and (78b), by method III:

$$\begin{array}{cc}\hfill x\left(t_2\right)& =\overline{\overline{x}}\left(t_2\right),\hfill \end{array}$$

(95a)

$$\begin{array}{cc}\hfill z\left(t_2\right)& =\overline{\overline{z}}\left(t_2\right),\hfill \end{array}$$

(95b)

$$\begin{array}{cc}\hfill u\left(t_2\right)& =\overline{\overline{u}}\left(t_2\right),\hfill \end{array}$$

(95c)

$$\begin{array}{cc}\hfill w\left(t_2\right)& =\overline{\overline{w}}\left(t_2\right),\hfill \end{array}$$

(95d)

leads to the third set of matching conditions:

$$\begin{array}{cc}\hfill (1/\ell )\sum _{n=0}^{\infty }{R}_n{(m_0\tau /c)}^n& =x_1/\ell -J_1+\sum _{n=2}^{\infty }{X}_n{\tau }^n,\hfill \end{array}$$

(96a)

$$\begin{array}{cc}\hfill J_{2}exp\left[(1/\ell )\sum _{n=0}^{\infty }{S}_n{(m_0\tau /c)}^n\right]& =exp\left(z_1/\ell \right)\left(1+\sum _{n=2}^{\infty }{Z}_n{\tau }^n\right),\hfill \end{array}$$

(96b)

$$\begin{array}{cc}\hfill \sum _{n=1}^{\infty }{R}_{n}n{(m_0\tau /c)}^{n-1}& =u_1-v_{★}{J}_3+v_{★}\sum _{n=2}^{\infty }{X}_{n}n{\tau }^{n-1},\hfill \end{array}$$

(96c)

$$\begin{array}{cc}\hfill J_4\sum _{n=1}^{\infty }{S}_{n}n{(m_0\tau /c)}^{n-1}& =w_{1}exp\left\{-\left[z\left(t\right)-z_1\right]/\ell \right\}\sum _{n=2}^{\infty }{Z}_{n}n{\tau }^{n-1}.\hfill \end{array}$$

(96d)

The matching time ${\tau }_2$ is a root of Equation (96a) if the downrange $x_1$ is specified at burn-out, of Equation (96b) if altitude $z_1$ is required, of Equation (96c) if horizontal velocity $u_1$ is specified, and of Equation (96d) combined with Equation (94b) viz.

$$J_4\left(\sum _{n=1}^{\infty }{S}_{n}n{(m_0\tau /c)}^{n-1}\right)\left(1+\sum _{k=2}^{\infty }{Z}_k{\tau }^k\right)=w_1{J}_2\left(\sum _{n=2}^{\infty }{Z}_{n}n{\tau }^{n-1}\right),$$

(97)

if the vertical velocity $u_1$ at burn-out is specified.

6.3. Third Trajectory with Specified Horizontal Velocity at Burn-Out

A specified horizontal velocity at burn-out leads to the condition of Equation (96c), which at the lowest order

$$N=1:R_1+2R_2(m_0\tau /c)=u_1-v_{★}{J}_3+2v_{★}{X}_2\tau ,$$

(98)

specifies the matching time

$$0⩽{\tau }_2=(u_1-v_{★}{J}_3-R_1)/(2R_2{m}_0/c-2X_2{v}_{★})⩽1-m_1/m_0,$$

(99)

and also the range of feasible horizontal velocities between the extreme values

$$\begin{array}{cc}\hfill u_{15}& =R_1+v_{★}{J}_3,\hfill \end{array}$$

(100a)

$$\begin{array}{cc}\hfill u_{16}& =u_{15}+(1-m_1/m_0)(2R_2{m}_0/c-2X_2{v}_{★}).\hfill \end{array}$$

(100b)

The data available in Table 2 is sufficient for these calculations, except for $J_3$, which follows from Equation (93c)

$$\begin{array}{cc}\hfill J_3& =2X_2(1-m_1/m_0),\hfill \end{array}$$

(101a)

$$\begin{array}{cc}\hfill J_3& =u_1/v_{★}.\hfill \end{array}$$

(101b)

The exact value of $J_3=78.8$ is given by Equation (101b) and appears in Table 3 using the reference velocity $v_{★}$ from Table 1 and the horizontal velocity at burn-out $u_1$. The lowest order approximation Equation (101a) to the series expansion Equation (96c) is calculated $J_3\sim 28.1$ from the lift-off $m_0$ and propellant mass $m_0-m_1$ in Table 1 and the parameter $X_2$ in Table 2. The first-order approximation is an underestimate 28.1 of the exact value 78.8, corresponding $28.1/78.8=0.36$ to 36% only and showing that better accuracy requires more terms of the series. Bearing in mind Equation (79a), one extreme of the horizontal burn-out velocity Equation (100a) is $u_{15}=v_{★}{J}_3=$5.35103$\mathrm{m}/\mathrm{s}$, which is the desired horizontal burn-out velocity. The other extreme adds an extra term involving the difference of $2R_2{m}_0/c$ and $2X_2{v}_{★}$. These two quantities are calculated in Table 3 using data from Table 1 and Table 2 and coincide within truncation error of 10 $\mathrm{m}$/$\mathrm{s}$. It follows that the second extreme value Equation (100b) coincides with the first Equation (100a) to within the truncation error. In fact, the two values are exactly equal, Equation (102b), due to the identity Equation (102a)

$$\begin{array}{cc}\hfill X_2{v}_{★}& =R_2{m}_0/c:\hfill \end{array}$$

(102a)

$$\begin{array}{cc}\hfill u_{15}& =J_3{v}_{★}=u_{16},\hfill \end{array}$$

(102b)

that can be proved by equating the initial horizontal acceleration calculated:

(i)

by method III using Equation (77a) with $n=2$ leading to Equation (102c)

$$\begin{array}{cc}\hfill {\ddot{x}}_0& =2R_2:\hfill \end{array}$$

(102c)

$$\begin{array}{cc}\hfill {\ddot{x}}_0& =\underset{t\to 0}{lim}\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2}=\ell {\left(\frac{\mathrm{d}\tau }{\mathrm{d}t}\right)}^2\underset{\tau \to 0}{lim}\frac{{\mathrm{d}}^2\chi }{\mathrm{d}{\tau }^2}=2X_2\ell {\left(\frac{c}{m_0}\right)}^2=\frac{2X_2{v}_{★}c}{m_0};\hfill \end{array}$$

(102d)

(ii)

by method I using Equations (11c), (17a), (11a) and (18a) leading to Equation (102d).

The equality of Equations (102c)≡(102d) proves Equation (102a) and hence, from Equation (100a) and (100b) also proves Equation (102b). The three examples of application of the TPBVP method, namely I + II in Section 4, III + II in Section 5 and I + III in Section 6 are compared in the conclusion Section 7.

7. Conclusions

The matching of the methods I and III has led to the unique solution of the TPBVP with the desired horizontal burn-out velocity

$$\mathrm{I}+\mathrm{III}{u}_1=5.35\times {10}^3 \mathrm{m}/\mathrm{s}=u_{15}=u_{16},$$

(103)

equal to the coincident extreme values in the Table 3. From Equations (102a) follows that the matching time Equation (99) is indeterminate ${\tau }_2=0/0$, so that the trajectories I and III can be matched at any time and thus coincide when truncated at the same order. Since the trajectories I and II are matched with the same horizontal and vertical coordinates and velocities, and these boundary conditions determine the whole trajectory, the two trajectories coincide. Although the trajectories I and III coincide, their series solutions converge more rapidly for short Equation (10) and long Equation (18b) and (18c) time respectively. Since both method III and method I involve power series of time, respectively Equations (76a) and (76b) and Equations (17a) and (17b), they lead to the same result if truncated at the same order N; the accuracy of the result will improve as the order N increases and the result is exact as $N\to \infty$ since both series converge.

Since the trajectories I and III coincide, the TPBVP combining either of them with the trajectory II gives the same result Equations (104a) and (104b)

$$\begin{array}{cc}\hfill \mathrm{I}+\mathrm{II}\equiv \mathrm{II}+\mathrm{III}:u_{11}& =u_{13}=-1.88\times {10}^4 \mathrm{m}/\mathrm{s},\hfill \end{array}$$

(104a)

$$\begin{array}{cc}\hfill u_{12}& =u_{14}=-7.98\times {10}^3 \mathrm{m}/\mathrm{s},\hfill \end{array}$$

(104b)

in agreement with Table 2. However, the TPBVP fails in this case because the extreme velocities are negative. In the successful application Equation (103) of the TPBVP the two trajectories I in Equation (10) and III in Equations (18b) and (18c) are ascending power series of time and coincide when both are truncated at the same order. The method II uses a power series not of time but rather of a shifted variable in Equations (40a) and (40b), and thus truncation at the same order does not give the same result; only the limit $N\to \infty$ would lead to the same results for the method II as for methods I and III. Since the TPBVP was applied at the first order, it applies to I + III but not when II is involved in II + I or II + III.

In order to illustrate the latter point, the methods I and II are compared next using the downrange. The horizontal displacement Equation (17a) for the method I leads to the Equations (105a) and (105b):

$$\begin{array}{cc}\hfill X_1& \equiv u_0/v_{★},\hfill \end{array}$$

(105a)

$$\begin{array}{cc}\hfill \chi \left(\tau \right)& =\sum _{n=1}^N{X}_n{\tau }^n,\hfill \end{array}$$

(105b)

that is a power series of time; the method II leads, Equations (24b) and (26a), to a power series Equation (106) of $(1-\tau )$

$$\chi \left(\tau \right)=\sum _{n=0}^N{P}_n{(1-\tau )}^n.$$

(106)

The two trajectories, Equations (106), (105a) and (105b) will coincide as $N\to \infty$ but not if truncated at a finite order N. The trajectory II in Equation (106) can be rewritten as a power series of $\tau$:

$$\begin{array}{cc}\hfill \chi \left(\tau \right)& =\sum _{n=0}^N{P}_n\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right){(-\tau )}^k\hfill \end{array}$$

(107a)

$$\begin{array}{cc}\hfill & =\sum _{k=0}^n{C}_k{\tau }^k\hfill \end{array}$$

(107b)

with coefficients

$$C_k={(-)}^k\sum _{n=0}^N{P}_n\left(\genfrac{}{}{0pt}{}{n}{k}\right),$$

(108)

for example Equations (109a) and (109b)

$$\begin{array}{cc}\hfill n=0& :0=C_0=\sum _{n=1}^{\infty }{P}_n=P_1+P_2+P_3+\cdots ,\hfill \end{array}$$

(109a)

$$\begin{array}{cc}\hfill n=1& :u_0/v_{★}=C_1=-\sum _{n=1}^{\infty }n{P}_n=-P1-2P_2-3P_3-\cdots .\hfill \end{array}$$

(109b)

At zero order $n=0$ from Equations (105b)≡(107b) follow the first equality in Equations (109a), and Equation (108) yields the second and third equalities in Equation (109a); at first order $n=1$ from Equations (105a) and (105b)≡(107b) follows the first equality of Equation (109b), and Equation (108) implies the second and third equalities in Equation (109b) in agreement with Equationss (28a) and (31a).

If the series Equation (106) is truncated at order $N=1$ it leads to Equation (110a)

$$\begin{array}{cc}\hfill P_0+P_1(1-\tau )& =P_0+P_1-P_1\tau \hfill \end{array}$$

(110a)

$$\begin{array}{cc}\hfill & \ne C_0+C_1\tau =(u_0/v_{★})\tau ,\hfill \end{array}$$

(110b)

that is distinct from Equations (107b) and (108) truncated at the same order $N=1$ leading by the first and third equalities in Equation (109a) to the second equality in Equation (109b). An accurate application of the TPBVP to the cases I + II in Section 4 or II + III in Section 5, would require truncation of Equations (105b), (107b) and (108) at the same order and would then lead to the same result as the TPBVP in the case I+II in Section 5. The less accurate approximation Equation (110a) highlights some features of the TPBVP like a range of velocities Equations (104a) and (104b) not found in the accurate approximation leading to a single value Equation (102b). Since both method III and method I involve power series of time, Equations (76a), (76b), (17a) and (17b), respectively, they lead to the same result if truncated at the same order N; the accuracy of the result will improve as the order N increases and the result is exact as $N\to \infty$ since both series converge.

The calculation of rocket trajectories in the atmosphere has been addressed in two cases: (a) for the single-point boundary-value problem (SPBVP) by calculating the trajectory from given initial conditions, using three distinct methods; (b) the combinations of two of the three methods, with one as ascending and the other as descending solution, provide six trajectories for the two-point boundary-value problem (TPBVP). The latter specifies in addition to the initial conditions one (not more) condition at burn-out, e.g., the desired horizontal velocity. The three methods I, II and III of solution of the SPBVP, and six combinations for solution of the TPBVP would lead to the same result if applied with full accuracy $N\to \infty$ but may give different trajectories if truncated at finite order N with shifted time variables (method II) rather than powers of time (methods I and III). The method of solving the TPBVP resolves three issues: (i) existence of the solution: whether the rocket has sufficient performance to be able to reach the desired burn-out condition; (ii) construction of the trajectory: if a solution exists the matching time is determined for the transition from the lower to the upper range allowing the construction of the complete trajectory; (iii) performance envelope: the range of burn-out conditions which is achievable, e.g., for the horizontal velocity at burn-out the minimum and maximum possible values. The solution to the three aspects (i) to (iii) of the TPBVP depend on the initial launch angle, as does the solution of the SPBVP.