Rapid Algorithm for Generating Entry Landing Footprints Satisfying the No-Fly Zone Constraint

An online estimation algorithm of landing footprints based on the drag acceleration-energy profile is proposed for an entry hypersonic vehicle. Firstly, based on the Evolved Acceleration Guidance Logic for Entry (EAGLE), drag acceleration-energy profiles are designed. To track the drag acceleration-energy profile obtained by the interpolation, a drag acceleration tracking law is designed. Secondly, based on the constraint model of the no-fly zone, flying around strategies are proposed for different conditions, and a reachable area algorithm is designed for no-fly zones. Additionally, by interpolating the minimum and maximum drag acceleration profiles, the terminal heading angle constraint is designed to realize the accurate calculation of the minimum and maximum downrange ranges by adjusting the sign of the bank angle. In this way, the distribution of landing footprints is more reasonable, and the boundary of a reachable area is more accurate. The simulation results under typical conditions indicate that the proposed method can calculate landing footprints for different situations rapidly and with the good adaptability.


Introduction
Landing footprints provide key information for the mission planning, such as the terminal area energy management (TAEM) guidance for a shuttle entry [1]. Due to a vehicle's high lift-drag ratio and the strong maneuverability, it forms a large range of landing footprints. The reachable area [2] formed by landing footprints is the maximum maneuvering range under given initial conditions, path constraints, and terminal constraints. During the entry process, the environment is complex, the flight state varies widely, and the shapes and sizes of reachable areas differ. Different from the case in a conventional mission planning, the initial condition has a great influence on a reachable area for an entry vehicle, and it has a great uncertainty. It is necessary to quickly obtain the reachable area based on landing footprints to select an appropriate landing site. Therefore, the reachable area needs to be calculated online, which requires the rapidity.
Hypersonic entry trajectory planning technology is the foundation of solving the reachable area. The main purpose of the entry trajectory planning system is to steer the vehicle flying from an entry point to a desired terminal area energy management interface. Since a well-planned trajectory is the key to the entry mission, a large amount of research on developing trajectory planning methods can be found. The sequential convex programming (SCP) method can be used to generate a reference trajectory under multiple constraints. Wang [3][4][5] used a convex quadratically constrained quadratic programming (QCQP) method to generate the optimal tracking guidance command. With the consideration of probabilistic constraints, a differentiable chance-constraint approximation method [6,7] is proposed. Therefore, the transformed optimal control model becomes solvable for trajectory optimization methods. Based on the optimal trajectories generated off-line by an indirect method, a deep neural network (DNN) model [8][9][10][11] is designed to enable the capability of optimal action predictions. Lu et al. [12] and Cheng et al. [13] introduced the Full Numerical Predictorcorrector Entry Guidance (FNPEG) into entry missions without the need for a reference trajectory or the vehicledependent planning. Chai et al. [14] proposed an optimal trajectory generation method by constructing an integrated framework with an inner gradient solver.
As the missions become more and more diversified, the vehicle needs to avoid threat areas and politically forbidden areas [15]. Thus, the vehicle needs to consider the constraint of the no-fly zone. In the trajectory optimization, the constraint of the no-fly zone has been studied [16][17][18][19]. Yang et al. [20] used an additional bank reversal scheme to shape the entry trajectory at the right time to avoid no-fly zones. Gao et al. [21] improved the tentacle-based bank angle transient lateral strategy to avoid static and dynamic no-fly zones. Yu et al. [22,23] proposed an analytical iteration scheme to plan the bank reversal sequence, for constraints of multiple no-fly zones and the flight-time. They also employed a new crossrange formula to the schedule proper bank reversal under the no-fly zone and terminal position constraints.
For a hypersonic vehicle, all possible trajectories can be obtained by trajectory optimization methods, and the reachable area can be formed. A well-designed solution to the reachable area should possess characteristics of the constraint satisfaction, the reliability, and the robustness. The existing reachable area algorithms can be divided into three categories: (1) Trajectory optimization techniques are employed directly in landing footprint generation algorithms, generally by a pseudospectral method (PSM) [24,25]. Fahroo et al. [26,27] employed the Legendre pseudospectral method for landing footprints. Other trajectory optimization methods such as the Gauss pseudospectral method [28][29][30] can also be applied to the landing footprint problem. The pseudospectral method has high accuracy, but the calculation takes a long time. (2) A simplified mathematical model based on the quasi-equilibrium glide condition (QEGC) can also generate landing footprints. Lu et al. [31,32] proposed a closestapproach method to obtain the maximum crossrange trajectory at a particular downrange by setting up a nonreachable virtual target. The centerpiece of this method is a closedform near-optimal bank angle control law that can solve landing footprints rapidly and reliably. However, this method is limited by a narrow model based on the assumption of the QEGC. Based on the same model, Liu et al. [33] obtained the maximum crossrange trajectory by the convex optimization. Li et al. [34] also proposed the selection scheme of a virtual target to make the method more applicable. Ngo and Doman [35] solved landing footprints under the failure of the control components based on the QEGC. The method only can obtain landing footprints corresponding to the maximum boundary due to the limitations of the simplified mathematical model of the quasi-equilibrium glide condition, and it takes a long time to discern the parameters. (3) Based on the EAGLE and the core idea that the trajectory range is inversely proportional to the drag acceleration, the entry corridor described by the drag acceleration profile can be designed. The maximum and minimum range trajectories and corresponding drag acceleration profiles can be obtained, respectively. Different drag acceleration profiles and different range trajectories can be obtained by linear interpolation. By keeping the sign of the bank angle positive or negative during a whole flight, Saraf et al. [36] can obtain terminal landing footprints. The method requires less calculation time, but its accuracy is lower than that of the pseudospectral method.
To solve landing footprints, a rapid algorithm for generating entry landing footprints satisfying the no-fly zone constraint is designed based on [36]. In this paper, we focus on the following three aspects of contribution: (1) we designed the drag acceleration profile tracking law and introduced the terminal heading angle constraint to realize the accurate calculation of initial landing footprints. (2) Considering the constraint of the no-fly zone, we proposed a fly-around strategy. Different relative trajectories were transformed into corresponding terminal heading angle constraints to realize the fast calculation of reachable areas under the constraint of the no-fly zone. (3) Based on the typical nonlinear model of an entry vehicle, the simulation of the proposed method was carried out. The results show that the method can obtain landing footprints based on the direct physical concept under the no-fly zone constraint. Compared with the traditional trajectory optimization, the calculation speed is faster and satisfies the requirement of the online calculation.
Although entry trajectory optimization problems for a hypersonic vehicle have been significantly addressed in the past years, the reachable area considered in this paper focuses on the generation algorithm of multiple trajectories. Different from [37] in which the no-fly zone constraint is transformed into the violation degree for the trajectory optimization, a fly-around strategy for the no-fly zone constraint is designed by adjusting the sign of the bank angle to avoid the threaten area in this paper. Compared with the reachable area algorithm [36], the no-fly zone constraint is added to make the application of the algorithm closer to the actual scene. In addition, the reachable area algorithm without the no-fly zone constraint is improved to distribute landing footprints more reasonable. Therefore, the results of the reachable area are more accurate compared with previous studies.
This paper is organized as follows: The landing footprint problem is briefly described, and the simulation model is introduced in Section 2. The solution to landing footprints is designed in Section 3. The verification of the simulation model is presented in Section 4. Thereafter, this paper closes with a summary of the most important conclusions in Section 5.

Problem Formulation
2.1. Entry Dynamics. Since the time and the energy of an unpowered hypersonic entry vehicle are monotonically changing, they can be chosen as independent variables. Due to constraints of the terminal velocity and altitude, the initial energy and terminal energy are limited in a small range. Therefore, the drag acceleration-energy profile is designed for the trajectory planning in the following sections in this paper. If the time is chosen as the independent variable, the 2 International Journal of Aerospace Engineering velocity can be calculated by integrating the velocity differential equation. However, considering the consistency with the drag acceleration-energy profile, the inverse mechanical energy E is an appropriate independent variable, which is determined by E = μ/r − V 2 /2, where r is the radial distance from the Earth's center to the vehicle, V is the velocity of the vehicle, and μ is the gravitational constant designed as 3.98606e14. Based on a sphere and irrational earth mode, the 3-dimensional equations of motion [38] are shown as follows: ψ′ = cos ψ tan ϕ cos γ rD where θ is the longitude, ϕ is the latitude, ψ is the heading angle, γ is the flight path angle, and σ is the bank angle. The acceleration of gravity g is modeled as follows: The aerodynamic lift acceleration L and the aerodynamic drag acceleration D are given by where q = ρV 2 /2 is the dynamic pressure, S is the reference area of the vehicle, C L is the aerodynamic lift coefficient, C D is the aerodynamic drag coefficient, α is the angle of attack, Ma = V/c is the Mach number, c is the velocity of sound, and m is the mass of the vehicle. For the entry trajectory planning, the angle of attack is usually specified as a function of the Mach number. The atmospheric density ρ based on the exponential atmosphere model is given by where ρ 0 = 1:225 kg/m 3 is the atmospheric density at sea level, h = r − r 0 is the altitude, r 0 = 6378136:3 m is the equatorial radius of the Earth, and h s = 7600 is the atmospheric scale height.

Constraints.
The vehicle has maximum path constraints for the heating rate, dynamic pressure, and aerodynamic acceleration [39], which are given by where _ Q is the heating rate, k Q is the thermal model constant, and n is the normal acceleration. Moreover, _ Q max , q max , and n max are the maximum values.
The terminal constraints for the altitude and the velocity are given by where Δh and ΔV are thresholds of the terminal constraint, r f and V f are the expected terminal values, and E f is the terminal inverse mechanical energy. When an area is not allowed to be passed due to regional factors, a no-fly zone should be established, as shown in Figure 1. The center of the no-fly zone is designed on the surface of the earth at the position Nðr nfz , θ nfz , ϕ nfz Þ, where r nfz , θ nfz , and ϕ nfz are the radial distance from the Earth's center to the center of the no-fly zone, the longitude, and the latitude of the center of the no-fly zone, respectively.
Thus, in the geocentric coordinate system, coordinates of the position N and the vehicle are given by x nfz = r nfz cos ϕ nfz cos θ nfz , y nfz = r nf z cos ϕ nfz sin θ nfz , where x nfz , y nfz , and z nfz are coordinates of the vehicle and x, y, and z are coordinates of the vehicle in the geocentric coordinate system. According to relations of the geocentric coordinate system and the north-up-east (NUE) coordinate system [40], we can determine the position of the vehicle in the NUE coordinate system as follows: where x 1 , y 1 , and z 1 are coordinates of the vehicle in the NUE coordinate system, and α nfz is the angle between the short half axis and the north direction.
The model of the no-fly zone constraint [41] is established as an infinite elliptical cylinder in the NUE coordinate system, which is given as where a and b are the long and short half-axes of the ellipse, respectively.

Algorithm Description
Lateral parameters of the trajectory are determined by the bank angle σ. Generally, the angle of attack is determined by a designed profile. Once the amplitude of the bank angle is determined, we can obtain L and L cos σ by (7). The amplitude of ψ ′ also can be obtained by (4), in which L sin σ is given by Therefore, lateral parameters of the trajectory are totally determined by the sign of the bank angle. Different drag acceleration profiles can be tracked by changing the magnitude of the bank angle. The sign of the bank angle should be kept consistent so that the landing footprint problem is translated into a maximum crossrange problem for different drag acceleration profiles. The terminal altitude and velocity constraints are satisfied by adjusting the sign of the terminal bank angle. Supplementary trajectories for maximum crossrange and minimum crossrange drag acceleration profiles make the reachable area more precise. For the no-fly zone constraint, lateral parameters need to be controlled by an appropriate strategy for the sign of the bank angle.

Drag Acceleration Profile Tracking Law.
In this section, we design a law to track drag acceleration profiles [42]. For D, we take the derivative with respect to E, which gives rise to Taking the derivative of D ′ with respect to E yields By substituting (5) into (24), we obtain L cos σ = ðD ″ − bÞ/a, which gives rise to where International Journal of Aerospace Engineering A linear feedback control law is used to track drag acceleration profiles, as follows: where ξ is the damp coefficient, ω n is the natural frequency of the tracking law, D is the real drag acceleration, and D ref is the reference drag acceleration. We can obtain D″ by (28) and substitute it into (25) to obtain the drag acceleration profile tracking law denoted as In this paper, simulations are carried out on the Common Aero Vehicle (CAV) [43]. Let the input reference drag acceleration D ref = 13 m/s 2 , the initial altitude h 0 = 35 km, the attack angle α = 20 ∘ , and the constant velocity V = 3000 m/s. Drag acceleration command tracking processes are shown in Figures 2 and 3. The tracking system can depress the overshoot but slow down the tracking speed for a higher damp coefficient ξ and a constant natural frequency ω n . If the damp coefficient is constant, the higher the natural frequency, the faster the tracking speed. However, if the natural frequency is too large, the high overshoot may be caused. Thus, tracking law coefficients are designed ξ = 2 and ω n = 5e − 6. If there is a white noise with a mean of 0 and a standard deviation of 0.5 in the drag acceleration command, the tracking process is shown in Figure 4, and the tracking system can filter the noise. The overall tracking speed is slow because the initial state of the tracking system is far away from the reference drag acceleration command, so the vehicle will take some time to fly to the appropriate altitude to track the reference command.

Equilibrium Glide Condition.
We assumed that the influence of the angle of attack on the range is ignored, and the angle of attack profile is determined by the Mach number, as follows: where α max and α max L/D are the maximum angle of attack and the maximum lift drag ratio angle of attack. Ma α1 and Ma α2 are the turning points of the attack profile according to the vehicle parameters, as shown in Figure 5.

International Journal of Aerospace Engineering
To achieve a smooth maximum downrange trajectory, we use the equilibrium glide condition by setting γ ≤ ε, where ε is a small vector. The bank angle tracking law is designed to obtain the relationship between the drag acceleration D and the inverse energy E at the lowest boundary of the drag profiles through numerical integration of the equations of motion. We use the feedback linearization to keep the flight path angle close to zero during the flight. To track the reference flight path angle γ ref , the feedback control law can be written as where K γ is a tracking coefficient.
To reduce bounces of the trajectory, we set the reference path angle γ ref = -0:01 ∘ . The reference path angle can be tracked in the initial phase, but the phugoid motion cannot be avoided during a whole flight with a fixed bank angle. The magnitude of the bank angle has been constrained to be 90 ∘ at an appropriate switching time to satisfy the terminal velocity constraint.
A simulation is carried out with parameters in section 4 to study the effect of different values for the tracking coefficient K γ on the equilibrium glide condition, as shown in Figure 6. The excessive coefficient may cause the divergence of the tracking system. According to Figure 6, we set K γ = 1 in this paper.

Path Constraint
Translation. The heat rate, dynamic pressure, and load constraints are translated into the corresponding drag profiles. By combining these drag acceleration profiles, the minimum downrange drag acceleration profile can be determined. We take the dynamic pressure constraint as an example to explain the transformation process.
(1) According to the initial and terminal states, the initial and terminal energy are obtained by E = μ/r − V 2 /2. We discretize the independent variable E into N intervals by ΔE = ðE f − E 0 Þ/N and obtain N + 1 node points denoted as fE 0 , E 1 , ⋯, E N g The dynamic pressure constraint is transformed into an equality constraint: f = q − q max (3) The values of r and V satisfying E = E i and f = 0 are determined by the secant method [44]. The detailed iterative process of the secant method is as follows: (a) In the kth iteration, r ðk−1Þ is the radial distance from the center of the Earth to the vehicle for the last iteration, V ðk−1Þ is obtained by E = μ/ r ðk−1Þ − ðV ðk−1Þ Þ 2 /2, and ρ ðk−1Þ is obtained by (6).
Especially in the first iteration, we can calculate these parameters according to the initial state of the vehicle. Based on the above definition, the (k − 1)th equality constraint becomesf ðk−1Þ = q ðk−1Þ − q max where ζ is a sufficiently small vector. If the condition in step d is met, we go to step e; otherwise, we go back to step a and start the next iteration using k = k + 1 (e) Once r ðk+1Þ is determined, V ðk+1Þ is also determined by the energy formula. A converged solution fr ðk+1Þ , V ðk+1Þ g satisfying E = E i and f = 0 is found (4) By performing iterative processes for each node fE 0 , E 1 , ⋯, E N g, we can obtain the minimum downrange drag acceleration profile The magnitude of the bank angle is constrained to 0 ∘ at an appropriate switching time to satisfy the terminal velocity constraint. The path constraints are translated into the minimum downrange drag acceleration profile, as shown in Figure 7. Since the minimum downrange drag acceleration profile cannot be differentiable at switching points where D ref ′ and D ref ″ cannot be obtained, a cubic spline is used to realize the second derivative of the drag acceleration profile continuously at the switching points.
For profiles between the maximum downrange and minimum downrange drag acceleration profiles, the interpolation equation is given by where D max and D min are the maximum and minimum drag International Journal of Aerospace Engineering acceleration profiles corresponding to the minimum and maximum range, respectively. i presents the ith drag acceleration profile, and n is the number of drag acceleration profiles.

Supplementary
Trajectories. Since the bank angle is kept positive or negative during flight, the interval of landing footprints between the maximum crossrange trajectories is considerable. For maximum downrange trajectories, we apply an interpolation between the heading angles ψ max 1 = ψ 0 and ψ max 2 = ψ egcf , where ψ 0 is the initial entry heading angle and ψ egcf is the terminal heading angle of the maximum downrange trajectory; thus, the supplementary terminal constraint ψ max is given by where i ψ represents the ith supplementary trajectory and n ψ is the number of supplementary trajectories. Moreover, another supplementary terminal constraint is given by ψ min = ψ max − δψ, and δψ is the threshold of the heading angle. Similarly, for minimum downrange trajectories, we apply an interpolation between the heading angles ψ min 1 = ψ 0 and ψ min 2 = ψ consf , where ψ consf is the terminal heading angle of the minimum downrange trajectory. Thus, the supplementary terminal constraint ψ min is given by Moreover, another supplementary terminal constraint is given by ψ max = ψ min + δψ.
For the supplementary terminal constraints, supplementary trajectories are designed by adjusting the sign of the bank angle, as given by where σ i−1 is the bank angle at the previous moment and σ i−1 is the bank angle at this moment. δψ is the threshold of the heading angle. If the value of δψ is too large, the result may be inaccurate. On the contrary, if the value is too small, it will cause the rapid reversal of the bank angle, which may cause the out of control.

Flying-Around Strategy for the No-Fly Zone Constraint.
For the no-fly zone constraint, according to the relationship among the no-fly zone, the position of the vehicle, and the heading angle of the vehicle, we designed a strategy to prevent the vehicle from entering the no-fly zone. In Section 2.2, we established a model of a no-fly zone. In this paper, the elliptical no-fly zone is simplified as a circular zone by the assignment of a = R z and b = R z in (21), in which R z is the radius of a circular no-fly zone. Assuming that the center of the no-fly zone Nðθ nfz , ϕ nfz Þ and the position of the vehicle Mðθ, ϕÞ are in the same horizontal plane, the relationship between N and M is shown in Figure 8. MA is the tangent from M to the no-fly zone, and A is the tangent point.
ψ A is the angle between the tangent MA and the north, and ψ Z is the angle between MN and the north. ψ is the redefined heading angle, where ψ = -ψ. ψ Z is solved by the following equation:

International Journal of Aerospace Engineering
If the vehicle passes the no-fly zone from the right side, the boundary of the heading angle is given by δ ψ is the threshold of the redefined heading angle, and j MNj is the distance between the point M and point N. The selection of the value for δ ψ is the same as δψ in Section 3.4.
If the vehicle passes the no-fly zone from the left side, the boundary of the heading angle is given by The sign of the bank angle is determined similarly to the equation in Section 3.4.

Landing Footprint Solution.
In this section, the reachable method is designed based on the condition that the no-fly zone is located on the right side of the initial heading direction of the vehicle. Under the condition of a left no-fly zone, a similar strategy can be designed but is not described in detail in this paper.
According to the relationship between the no-fly zone and landing footprints, the solution is divided into five cases, as follows: Case 1, without no-fly zones; Case 2, the no-fly zone is located entirely outside the area of landing footprints; Case 3, the no-fly zone is located entirely inside the area of landing footprints; Case 4, the no-fly zone intersects the edges of the supplementary minimum downrange landing footprints, while the center of the no-fly zone is located outside the area of landing footprints; and Case 5, the no-fly zone intersects the edges of the supplementary minimum downrange landing footprints, while the center of the no-fly zone is located inside the area of landing footprints. For the cases in which the no-fly zone intersects the edges of the supplementary maximum downrange or interpolated maximum crossrange landing footprints, the no-fly zone has a small effect on landing footprints, which is not described in this paper. The set of landing footprints is denoted as J, and the trajectory is denoted as j. The steps based on EAGLE for solving the landing footprints for each case are described as follows: (a) Initial landing footprints: interpolated drag acceleration profiles between the maximum downrange and minimum downrange drag acceleration profiles determined by (34) are tracked by the drag acceleration profile tracking law presented in Section 3.1. By keeping the bank angle positive or negative during the whole flight, we can obtain sets of interpolated landing footprints J 1 and J 4 . The trajectory is denoted as j 2 or j 4 when it corresponds to the positive or negative bank angle tracking from the maximum downrange drag acceleration profile. Similarly, the trajectory is denoted as j 1 or j 5 when it corresponds to the positive or negative bank angle tracking from the minimum downrange drag acceleration profile. The trajectory with the terminal heading angle ψ 0 corresponding to the maximum downrange drag acceleration profile is denoted as j 3 and that corresponding to the minimum downrange drag acceleration profile is denoted as j 6 (b) Supplementary landing footprints: we can obtain the set of interpolated landing footprints J 2 or J 3 constrained to the supplementary terminal heading angle by (35), in which the value of ψ egcf is ψ 2 (the terminal heading angle of j 2 ) or ψ 4 (the terminal heading angle of j 4 ). The sign of the bank angle is given by (37). Similarly, we can obtain the set of interpolated landing footprints J 5 or J 6 in the same way, corresponding to the minimum downrange drag acceleration profile   International Journal of Aerospace Engineering Case 2.
(a) Initial landing footprints: sets of landing footprints J 2 , J 3 , J 4 , and J 5 are obtained by steps a and b in Case 1. For trajectories not affected by the no-fly zone in J 1 , we denote the trajectory closest to the no-fly zone from the left as j 7 and closest from the right as j 11 . Landing footprints calculated from j 7 to j 2 are denoted as J 7 , and landing footprints from j 11 and j 1 are denoted as J 16 (b) Modified landing footprints for the left side: use the left flying-around strategy in Section 3.5 for trajectories between j 7 and j 1 and denote these landing footprints as J 8 . The terminal heading angles of the flying-around trajectories j 8 and j 9 are denoted as ψ 8 and ψ 9 , corresponding to the original trajectories j 7 and j 1 . Where the trajectory j 9 and the no-fly zone intersect at the point of tangency, we denote the heading angle here as ψ 10 . The interpolation between ψ 9 and ψ 10 is regarded as the terminal heading angle constraint. Thus, we can obtain the set of landing footprints J 9 by the same method used in Section 3.4. Similarly, we can obtain set J 10 by the interpolation of ψ 0 and ψ 10 (c) Modified landing footprints for the right side: landing footprints calculated from j 11 to j 1 are denoted as J 16 . Use the right flying-around strategy in Section 3.5 for the trajectories between j 11 and j 1 and denote these landing footprints as J 13 . The terminal heading angle of flying-around trajectory j 13 or j 14 is denoted as ψ 13 or ψ 14 , corresponding to the original trajectory j 11 or j 1 . The terminal heading angle of the trajectory j 11 or j 1 is denoted as ψ 11 or ψ 1 . Where the trajectory j 13 or j 14   (a) Initial landing footprints: sets of landing footprints J 2 , J 3 , J 4 , J 5 , J 6 , J 7 , and J 16 are obtained through steps a and b in Case 2 (b) Modified landing footprints for the left side: we can obtain set J 8 through a process which is similar to step b in Case 2 but slightly different. If a trajectory will not reach the no-fly zone at the point of tangency, the landing footprint of this trajectory will no longer belong to J 8 . In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in J 8 is denoted as j 9 (c) Modified landing footprints for the right side: the trajectories j 11 , j 12 , j 13 , set J 11 , and set J 12 are given by the same step in Case 2. We can obtain set J 13 by a process that is similar to step c in Case 2 but slightly different. If a trajectory will not reach the no-fly zone at the point of tangency, the landing footprint of this trajectory will no longer belong to J 13 . In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in J 13 is denoted as j 14 (d) Part of the no-fly zone: select points p 1 and p 2 closest to the landing footprints of j 9 and j 14 on the edge of      (d) Modified landing footprints: the trajectory corresponding to the landing footprint closest to the no-fly zone from the left side is denoted as j 17 and that closest from the right side is denoted as j 18 . The landing footprints of the trajectories between j 6 and j 17 are denoted as J 17 , and the landing footprints of the trajectories between j 18 and j 1 are denoted as J 18 (e) Part of the no-fly zone: select points p 1 and p 3 closest to the landing footprints of j 9 and j 17 on the edge of the no-fly zone and denote the points between p 1 and p 3 as P 1 . Select points p 2 and p 4 closest to the landing footprints of j 14 and j 18 on the edge of the no-fly zone and denote the points between p 2 and p 4 as P 2 (f) Final landing footprint: the area of the landing footprints consists of two closed figures constructed by connecting sets The above solution is based on the assumption that the no-fly zone affects sets of landing footprints J 1 , J 4 , J 5 , and J 6 . If the no-fly zone has an impact on J 2 or J 3 , the key algorithm is the same as that in the above study and is not described in this paper. The set of landing footprints corresponds to Figure 13.

Numerical Examples
In this section, a 3-degrees model of the CAV is used for the landing footprint simulation. To verify the algorithm,    International Journal of Aerospace Engineering MATLAB is used to model the problem. The mass of the vehicle is 907.2 kg, and the aerodynamic reference area is 0.484 m 2 . The path constraints are set as _ Q max = 2000 kW/m 2 , q max = 100 kPa, n max = 3 g, and k Q = 1e-4. The initial states of the vehicle are shown in Table 1. The terminal constraints for the altitude and the velocity are 3000 m/s and 30 km, respectively. And thresholds of the terminal constraints are designed as Δh = 100 m and ΔV = 50 m/s. The angles of attack profile parameters are designed as α max = 20 ∘ , α max L/D = 20 ∘ , Ma α1 = 12, and Ma α2 = 12. The tracking law coefficients for the drag acceleration command tracking process are designed ξ = 1 and ω n = 5e − 6. Similarly, the tracking law coefficient for the equilibrium glide condition is designed K γ = 1. The threshold of the heading angle in Section 3.4 and Section 3.5 is designed as δψ = 2 ∘ . The number of drag acceleration profiles is designed as n = 10, and the number of supplementary trajectories is designed as n ψ = 10. Thus, for Case 1, the number of total trajectories is 80. For comparison, the number of total trajectories for the Gauss pseudospectral method is also set as 80.

Without No-Fly Zones.
In this section, the no-fly zone constraint is not added, corresponding to Case 1, and the effectiveness of the method proposed in this paper based on EAGLE is verified. Figures 14 and 15 show the tracking ability of the minimum and maximum acceleration profiles. At the initial point and the turning point of the attack angle profile, there are little lags in the tracking process, but the overall tracking is still good. Table 2 shows the computation times of different methods. The solving time for the method proposed in this paper is 13.42 s, which is suitable for the online rapid engineering, in contrast to the 429.93 s by the Gauss pseudospectral method, which is suitable for accurate offline engineering. Figure 16 compares landing footprints obtained by the EAGLE and Gauss pseudospectral methods, and we can see that landing footprints obtained by EAGLE are significantly Tracking value Reference value Figure 14: Minimum drag acceleration profile tracking ability. Tracking value Reference value Figure 15: Maximum drag acceleration profile tracking ability.   11 International Journal of Aerospace Engineering smaller. This is because the trajectory based on the drag acceleration profile has not been optimized, the calculation process of landing footprints obtained by EAGLE has some approximate simplification steps, and the equilibrium glide condition is added. Figure 17 shows altitude profiles from the equilibrium glide condition to transformed path       constraints. Due to the flight capacity and terminal constraints, the maximum range trajectory maintains the equilibrium glide condition only for a while, and the bouncing of the trajectory cannot be avoided during the whole process. Figures 18 and 19 show drag acceleration profiles and corresponding lift acceleration profiles. Drag acceleration profiles        Figure 24 shows the dispersion of the terminal altitude within 12 m, and Figure 25 shows the dispersion of the terminal velocity within 10 m/s, which satisfy terminal constraints. Figures 26-28 show all path constraint profiles, which also satisfy proposed path constraints.

With
No-Fly Zone. In this section, the no-fly zone constraint is added. Centers of designed no-fly zones are set as ð1:4 ∘ , 10 ∘ Þ, ð3 ∘ , 12:5 ∘ Þ, ð2 ∘ , 10:6 ∘ Þ, and ð2 ∘ , 11:1 ∘ Þ, corresponding to Cases 2-5 in Section 3.6, respectively. The radius of no-fly zones is designed as 40 km, and the angle between the short half axis and the north direction is designed as α nfz = 0 ∘ . The number of the interpolated trajectories is designed as 5, which are supplementary trajectories for the flying-around strategy in Section 3.6. By applying the flying-around strategy described in Section 3.5 and the landing footprint solution presented in Section 3.6, we can obtain final modified landing footprints shown in Figures 29-32. To see the effect more clearly, landing footprints near the no-fly zone are partially magnified. Except for the Case 3, final landing footprints consist of two closed figures. Table 3 shows the computation times of different cases. Trajectories affected by the no-fly zone mean flying-around trajectories and supplementary trajectories for the flying-around strategy. The solving time for these cases is short, which is suitable for the online rapid engineering.

Conclusions
In this paper, a fast method to generate the reachable area for entry vehicles is developed, which is subjected to inequality entry path constraints, terminal constraints, and the no-fly zone constraint. The reachable area algorithm follows the basic reference profile tracking scheme and consists of three modules: drag acceleration-energy profiles and the corresponding tracking law are designed to achieve initial landing footprints; a flying-around strategy is designed to avoid the no-fly zone; the solution of the reachable area is proposed by the reasonable distribution of landing footprints. Compared with the Gauss pseudospectral method, the rapidity and the effectiveness of the proposed solution are comprehensively analyzed. The results show that the solution has the good applicability and dynamic characteristics. Therefore, the online reachable area solution developed in this paper provides the key information for the trajectory planning of hypersonic vehicles. In addition, the method has some shortcomings in terms of accuracy. Future studies will try to find a highly accurate online landing footprint solution for no-fly zones.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.