Modelling of aircraft trajectories for emergency landing using kinematoid chains

Most methods that compute trajectories for un- or low-powered flight operate under simplifying assumptions such as constant curve radii and wind conditions. Likewise, changes of air density with altitude that lead to significant differences between equivalent airspeed (EAS) und true airspeed (TAS) are often not considered. Some approaches are based on Dubins paths, which are introduced in Dubins LE (Am J Math 79: 497, 1957). They combine three sections to form a trajectory, which is the shortest from a given start to an end position (In the original work, the position extended by the heading, is referred to as configuration. Since configuration has a different connotation in our context, we use the term state, which can contain other parameters in addition to the position and the heading, e.g. orientation, configuration of landing flaps, landing gear, etc.). A maximum distant landing spot can be reached this way. Often, the targeted landing spot is closer to the aircraft. If it is approached using a Dubins Path, the excess height must be dealt with. Here, we present a method addressing the problem directly, namely finding a trajectory which reduces the excess height over its entire length. Furthermore, it takes spatial and temporal changes of wind and air density into account. Several conditions influence the final shape of the trajectory. For example, avoidance of obstacles and predefined areas is easily achieved. Our method is motivated by kinematic chains, which are used in robotics and computer animation. We extend and modify this principle by incrementally transferring start and end states of the trajectory, modelled as state vectors, into each other. The resulting intermediate states form ends of chain links. To connect initial and final states through the resulting chain, we solve the inverse kinematic problem known from robotics. We extend it by several conditions, which are derived from the flight mechanical characteristics of the modeled aircraft on the one hand and from the desired properties of the trajectory on the other. Using practical examples, we will show the performance of this method, which we have efficiently implemented on off-the-shelf hardware. The method is suitable for systems that assist in the event of engine failures as well as for modeling planned un- or low-powered flights like continuous descent approaches or return flights of space gliders.


Introduction
Un-or low-powered, flights are carried out in aviation for various reasons. On one hand, there are emergencies in which the engine(s) no longer deliver power. On the other hand, there are deliberate low-powered and gliding flights, for example as part of a continuous descent operations (CDO) or the return flight of an unpowered spacecraft.
In all these cases, the gliding flight trajectory should be modeled as accurately as possible. Given conditions such as the changing air density during the descent and complex wind conditions should also be considered. Furthermore, the method should allow the user to influence the shape of the resulting trajectory. It should be able to run efficiently on affordable hardware because it shall operate in flight on an embedded system in real time and has to converge in a deterministic time frame. In our consideration, we assume that the available potential and kinetic energy of the aircraft is sufficient to reach the landing field. This paper is structured as follows: In Sect. 2, we give an overview of existing methods and discuss their advantages and disadvantages. Section 3 describes in detail the method we have developed. The underlying data structures and further developments of the original approach are presented.
In this work, we deliberately limited ourselves to dissipating the existing energy of the aircraft by means of a suitable trajectory without using aerodynamic aids. On one hand, the performance of the method as such can be studied better. On the other hand, these aids are not always available in emergencies.
In Sect. 4, we show how interesting use cases can be calculated with the presented method. Section 5 summarizes what has been achieved while Sect. 6 details how the properties of kinematoid chains could be further explored and their scope broadened, for example to better model configuration changes as well as powered flight.

Related work
Emergency landing assistance is systematically researched for about 20 years. We can distinguish two branches of methods, that provide either advice about the reachability range or effective gliding paths toward possible landing sites. The first category is used as simplified emergency support in commercial navigation solutions like Garmin's SmartGlide [2], ForeFlight [3] and SkyDemon [4]. Those systems provide a glide range ring that is drawn just around the current aircraft's position while adjusting displacements according to the wind and surrounding surface. Unfortunately, the difference between the initial heading of the aircraft at the emergency point and the emergency runway heading is not considered by those glide advisors. Even if an emergency landing site lies within the glide ring range, additional turns are required to adjust the difference between the aircraft heading at the time when the emergency occurs and the heading of the emergency landing strip. Unfortunately, the altitude loss is not considered during those turns. As described in [5], the energy footprint 1 of the airplane will become approximately elliptical. Thus, the glide ring will shrink during the glide along a path that has to be selected solely by the pilot. In contrast, current research focuses on the critical task of computing an effective gliding path that will safely bring the airplane over the threshold of an emergency landing field. By such a 3D landing approach division, the excess altitude must not only be wisely used to reach the threshold in an appropriate height but also to compensate the wind offset.
In [1] Lester E. Dubins described a method for finding a shortest path between two two-dimensional vectors. He found out that there are in total six paths to connect the start and destination vector by a combination of circle segments and straight tangents to it. In the context of landing assistance, we have to extend Dubins original approach to three dimensions. While the aircraft glides laterally along a 2D Dubins path, it consumes potential energy of the excess altitude. Fortunately, the loss of altitude is proportional to the 2D path length and can be easily described by the glide ratio of the aircraft when flying with the speed for optimal glide ratio in no wind conditions. The Dubins path is always of minimum length and its path length can be much shorter than the available glide length, which is the aircraft's current altitude over ground times the glide ratio of the aircraft. To adapt the approach path to the available glide length one could change the radii of the circles, extend the tangent (straight) segments between circles and/or extend the final approach appropriately. Examples for adaptions of this kind can be found in [6,7]. In [8], it is proposed to reduce slight excess altitude by lengthening a straight final approach. In case of larger excesses of potential energy, further waypoints are added to the Dubins path to fly detours.
Reference [9] varies the strategy to dissipate different magnitudes of potential energy. In case of small surpluses, the approach trajectory is lengthened. In case of a large amount of excess height, the aircraft is first flown close to the landing spot and then reduces the altitude by flying full circles. The advantage over the detour strategy is an earlier proximity to the landing spot. The number of full circles can be reduced in the event of incorrectly estimated gliding performance or wind conditions. Since full circles can only be used to reduce height in discrete steps, a further procedure is required to dissipate the remaining height. In the method of Ref. [9], a combination with an extension of the final approach is proposed.
For finding a valid glide trajectory, it is important to compensate the influence of the wind by choosing the route appropriately. Related work either just considered a combination of primitive manoeuvres in windless situations [10] or their authors integrated the wind directly into Dubins-based route calculation. Examples for the last mentioned methods are modifications of Dubins curves in the earth coordinate system. Due to the wind's influence so called trochoids develop, see [6] for example. Unfortunately, trochoids must be calculated starting from the emergency's location as well as from the threshold of the destination runway. Moreover, a connecting tangent must be fitted at both ends which can be numerically demanding and thus time consuming.
References [11,12] present an alternative solution that is much faster and that can be used with any trajectory calculation method. By estimating the time for the glide trajectory and moving the target of the approach opposite to the current wind vector we get a substitute target. The route to that target can easily be calculated as in the windless case and due to the blowing wind the aircraft eventually arrives at the requested runway threshold. Instead of defining the flight path relative to the earth, it is defined relative to the air mass. In this way, the moving target approach can easily model wind vector fields.
There are also various other path planning techniques from the field of robotics, e.g. the Rapid-exploring Random Tree (RRT) methods and its variants, introduced 1998 by LaValle [13]. Unfortunately, there is a big computational demand and the resulting trajectories are not smooth enough for gliding approaches. Also, Genetic Algorithms [14] have been proposed but due to the huge computational demand of those techniques, it is impossible to apply them under real-time conditions. Reinforcement Learning (RL) has been proposed for a wide variety of applications in aviation [15].
The method presented here has the following advantages over aforementioned ones: • Methods based on Dubins paths first compute the shortest path to the target which they then extend as required. In contrast, our method addresses the problem directly by computing a trajectory of appropriate length. • Parameters that vary over the trajectory are taken into account. An example is the consideration of changing air density so that the relationship between EAS and TAS can be modeled realistically. • Wind conditions of any complexity, in particular changes in direction and magnitude over altitude, can be modeled. The complexity of wind conditions has no influence on the calculation complexity. • The method is able to respect complex constraints influencing the shape of the resulting trajectory. Among other things, this makes it possible to avoid obstacles.

Our method
A method for calculating glide trajectories must be able to generate a trajectory to a landing spot within a reachable area, given height and (unpowered) flight characteristics of an aircraft. If the landing spot lies at the edge/rim of the area, the trajectory corresponds to the shortest path to that spot. Methods based on Dubins paths achieve this and only this. In most cases, however, it is unlikely that a chosen landing spot can be approached via the shortest route. The probability that the landing spot is inside the area that can be reached from the current position and height is far higher than a target position exactly on the edge. The shortest path to an interior landing spot would arrive there with excess height but the related excess energy must be eliminated before landing.
Here, we want to model trajectories that reduce existing height as precisely as possible over their entire length.
These are end-of-trajectory conditions and conditions that must be met throughout the flight.
The search is therefore not for the shortest trajectory but for the one that reaches the landing spot at the requested height and fulfills a given set of conditions. These are conditions regarding the end of the trajectory as well as the course of the flight. They pertain to aircraft configuration as well as flight state and include operational limits among others.
The problem to be solved is the following: An initial and an end state of a trajectory, as well as a sufficiently accurate aircraft performance model (APM) are given. Initial and final states can be represented as vectors in a state space Σ which is a Cartesian product. The components of this state space are at least: position (X, Y and Z coordinates), heading and a time stamp. For precise models, these can be expanded to include variables such as horizontal and vertical speed, bank angle, aircraft configuration (e.g. whether flaps are set and state of the landing gear), and environmental information (air pressure, wind direction and speed).
The initial state is to be transformed into the final state in such a way that each intermediate state, which also represents a state vector in the state space Σ , satisfies a set of conditions that are still to be defined. These can be conditions on the individual state itself, but also on the relation between subsequent states. Examples of conditions on individual states are compliance with maximum and minimum speeds. However, compliance with a constant fixed speed can also be made a condition, for example the speed for best gliding. The speed for the best gliding in turns depends on the aircraft configuration and bank angle. This example shows that the condition itself can depend on the parameters of the state vector.
An example of a condition on relations of subsequent states that must always be met is the reachability of the X, Y and Z coordinates with given speeds and temporal distance of subsequent states.
One can think of the sequence of consecutive states as a chain. A chain link (hereafter chain segment) is formed by a state and the transition to the subsequent state. Since-as mentioned above-the accessibility of the subsequent state is a basic condition, the individual states form a coherent structure. As will be shown, this is very similar to kinematic chains that are known and well studied in kinematics, robotics and computer graphics. However, since we give up elementary properties of classical kinematic chains, such as a fixed length of the chain segments, we call this structure a kinematoid chain. (See Fig. 1

.)
Glide trajectories can be modeled by kinematoid chains as follows: We assume that no wind is present. In Sect. 3.3, we will show how we can also take complex wind conditions into account with a method based on the moving target principle.
We start by defining n chain segments. Since we are modeling a continuous trajectory using a discrete data structure, a discretization error is to be expected. To minimize it, the number of chain segments should be chosen sufficiently large. We will later develop the structure of the kinematoid chain segments in such a way that small discretization errors can be realized even with relatively few chain segments. The total altitude difference Δalt between the initial and final state is now divided into n equidistant sections. Each of these sections is modeled by a chain segment, as visualized in Fig. 1. We determine the length of each chain segment, i.e. the distance between two consecutive states, from the assumed horizontal speed (for example the speed for the optimal glide ratio) and the gliding behavior defined in the APM. This gliding performance is dependent, among other things, on the bank angle. Since we do not yet know the final curvature of the trajectory, we will initially assume gliding performance in straight flight without banking. As the trajectory gradually changes towards its final shape in an iterative process, the chain segment length will be adapted. If the bank angle for the respective chain segment increases, the sink rate on this section increases and the height section, which remains constant, is reached earlier which reduces the length of the chain segment.
At the beginning of an iteration, we line up all chain segments starting from the initial state in the direction of the heading of the initial state. The position of the end of the chain will usually differ from the desired end state. (Fig. 2) To connect the end of the chain with the final state while keeping the beginning of the chain at the initial state and maintaining the connectedness of the chain, the so-called inverse kinematic problem has to be solved. There are many well-studied methods from kinematics and robotics. The FABRIK. 2 algorithm [16] proves to be well suited for the present case. Here, it has been expanded to include (2) Δalt seg = alt begin − alt end n compliance with a freely definable set of conditions. It is an iterative process. After initialization, the FABRIK algorithm runs in two phases, the forward and the backward phase. After the kinematoid chain has been initialized in the direction of flight starting from the aircraft as shown in Fig. 2, the first backward phase is started. The head of the nth chain segment is connected to the landing spot, its tail points towards the head of the remaining chain (Fig. 3a). Now, all chain segments are re-hooked one after the other from one part of the chain to the other. The tail of a chain segment that has just been re-hung always points in the direction of the head of the next chain segment to be re-hung (Fig. 3b). When all chain segments have been re-hung, the head of the chain is connected to the landing spot, but there is a gap between the tail of the chain and the aircraft position (Fig. 3c). Now, the forward phase begins. The first chain segment is moved back to the aircraft position with its tail. Its head points to the rest of the chain. After the forward phase for the entire chain has been traversed, the chain will be reconnected to the aircraft position at the tail, but a gap will have reappeared between its head and the landing position. In contrast to the original form of the algorithm, the following calculations are performed for each chain segment seg ∈ {0, ..., n − 1} : The mean altitude of a chain segment is calculated by: We use a simple APM that, for each aircraft configuration that defines the EAS for the best glide performance and the associated glide ratio E(0) and its reciprocal E(0) = 1 E(0) for straight flight without bank. The glide ratio given a bank angle is then calculated as follows: However, more complex APMs could also be used, for example an APM, that additionally considers the impact of the bank angle on the optimal EAS.
The true airspeed (TAS) is calculated from the EAS and the height alt(seg) via the air density. The air density can either-similar to the wind data-be retrieved from data services such as [17] or determined approximately according to the ICAO Standard Atmosphere, see [18].
with: TAS(seg) average true airspeed of segment 0 air density at sea level (seg) average air density of segment EAS(seg) equivalent airspeed For a very simple model, the relationship between EAS and TAS for a given pressure altitude can be approximated described in [19].
Using the flight path angle , we derive a description of the mean descent rate along a segment as a function of glide ratio E and TAS: The time required for the aircraft to fly the distance of the respective chain segment is calculated from the rate of descent and the height difference to be overcome per chain segment (which by definition is the same and constant for all chain segments).
The length of the chain segments under no-wind conditions is the product of TAS and time.
If, in the current iteration step of the FABRIK algorithm, a smaller curve radius is modeled by a chain segment than in the previous iteration step, this results in an increased bank angle. Due to the higher sink rate, this means that the height difference to be overcome by the chain segment is flown along a shorter horizontal distance. The chain segment is thus shorter than in the previous iteration step.
The minimum radius for a coordinated turn is calculated from the TAS and the specified maximum bank angle with: alt(seg) pressure altitude in feet The maximum course change that can be modeled with this chain segment is calculated from this minimum flyable radius and the length of the chain segment. We are still assuming here that we are modeling the trajectory with straight chain segments. The length of the chain segment corresponds to l straight in Fig. 4. The maximum possible course change ΔΨ max (seg) is then calculated by: We will show in section 3.2 that a chain segment can also be modeled using segments of a circular arc. In this case, the length of the chain segment corresponds to l circ in Fig. 4. ΔΨ max (seg) is calculated according to This ΔΨ max (seg) is the maximum possible course change (in radians) under the given conditions (specified: EAS, max , calculated: alt(seg), TAS(seg)).
If the FABRIK algorithm computes a larger course change, it is limited to ΔΨ max (seg) . The curvature of the chain generated by the FABRIK algorithm is thus limited according to the individual parameters of the respective chain segment.
The described method calculates the individual values for the coordinates, the horizontal and vertical speed, bank angle, aircraft configuration for each chain segment and stores the values in the data structure of the respective chain segment. This ensures that each individual chain segment as well as the relations between subsequent links correspond to the defined set of conditions and the underlying APM.
After each iteration, a residual error remains, either between the initial state and the beginning of the chain, or between the final state and the end of the chain. In [16], it was shown for the original FABRIK algorithm that this error decreases strictly monotonically with the number of iterations. Therefore, the method converges (Fig. 3(d)). This proof is still to be provided for the extension to include the conditions described above. However, all the experiments we conducted showed a convergence. Of course, this only applies for cases in which at least one trajectory exists under the defined conditions. We stop the iterations when the residual error is sufficiently small and output the trajectory found. Otherwise, the method will be terminated unsuccessfully after a maximum number of iterations. We calculate the residual error from the maximum of two possible gaps. One between the aircraft position and the chain tail (tail error, Fig. 3c) and the other between the landing spot and the chain head (head error, analog to tail error).

Data structure
Classic kinematic chains are defined by straight chain segments and the joints located between them (Fig. 5).
This data structure has a disadvantage: To be able to calculate the heading change from one chain segment to the next, not only the heading of the current chain segment must be accessed, but also the heading of the next chain segment. However, this may not have been calculated yet.
To solve this problem, we changed the data structure of a chain segment so that an entity now consists of two halves of the original chain segments and the joint in between. (See Fig. 6.) However, this data structure, which now consists of two straight lines and an articulated joint, still suffers from the discretization error mentioned above. The magnitude of this error depends on the number n of chain segments and the shape of the trajectory. In principle, the discretization error can be made arbitrarily small by increasing the number of chain segments. However, in the case of long and complex trajectories, this can lead to a very high number of chain segments. This in turn is associated with a considerable computational effort. A more elegant way of approximating the trajectory, which always runs continuously in the real  world, with less error is to model the chain segments not with straight but with circular arc sections. They approximate real trajectories much better than straight lines and are similarly easy to handle (Fig. 7).
These circular arc chain segments are now lined up in such a way that, at the connection point of two chain segments, their first derivative is identical. This C 1 -continuity is shown in in Fig. 8.
As can be seen in Fig. 9 using the example of a trajectory in the form of an Archimedean spiral, 3 the discretization error is drastically reduced compared to using straight chain segments. The blue curve shows the deviation of the length of the chain from the actual length of the trajectory when using straight chain segments, the magenta colored one when using circular arc segments as a function of the number n of chain segments.

Trajectory shaping
The calculation of the parameters for each chain segment described above can now be expanded to include additional conditions that offer the opportunity to shape the trajectory. From the mostly infinite number of solutions to the inverse kinematic problem, a solution adapted to the desired use case is determined. For a trajectory for an emergency landing approach, it is desirable not to fly turns unnecessarily tight so that the approach can be corrected with tighter turns if necessary. This can be achieved by making the condition for the maximum permissible bank angle dynamic. For this purpose, in the first iteration steps, a significantly smaller bank angle than the maximum possible is allowed. If no trajectory exists between the initial and final state the permissible bank angle may be slightly increased in each iteration step of the FABRIK algorithm. The kinematoid chain is, therefore, stiffer at the beginning and adapted to its full flexibility during the iteration process.   In addition, it is possible to produce a straight final approach of a certain length. This is achieved by making the permissible maximum bank angle dependent on the distance between the chain segment and the final state (landing spot). For all chain segments below a certain distance, the permissible heading change is set to 0. This can still be modified by gradually reducing the number of the straight final segments if otherwise no trajectory can be formed. A straight final is only formed if possible and to the extent that one can afford it.
The form of a trajectory can also be influenced in a targeted manner towards certain way points or away from defined areas like obstacles. To avoid certain areas, a trajectory is first calculated using the above-mentioned method. If one was successfully formed, all chain segments that are inside an area to be avoided are moved by a certain amount away from it. The avoidance area is thus modeled as a force repelling the chain segments. (See Fig. 10.) Then, the FABRIK algorithm is executed and restores all integrity constraints of the chain. It must be ensured that the changes caused by the rejection are not completely reversed. After several iteration steps, in which the repulsive forces and the FABRIK algorithm alternately act on the chain, a trajectory is created which avoids the obstacle.
In [11], the so-called ≫Moving Target Principle≪ (MTP) is proposed. In its original form, the offset vector, resulting by the wind drift along the glide path, is calculated from the average wind vector and the exposure time. Then, the landing point is shifted (in the earth reference system) in the negative direction of this displacement vector and a virtual landing point is obtained. The virtual landing point is approached as in no-wind conditions since the wind during gliding flight shifts the aircraft towards the desired landing spot. This saves complex path calculations under the influence of wind, which lead to cycloid curves (trochoids) in the earth reference system. While path calculations based on Dubins paths usually only consider a constant wind vector, the MTP can easily be applied piecewise to wind vector fields of any complexity in connection with kinematoid chains. For this purpose, an individual wind offset vector is calculated for each chain segment according to its position, height, length and local wind vector and added to form a total vector. A virtual approach point is then calculated from this total vector as described above. At the time the trajectory starts, this virtual approach point is offset from the actual approach point against the wind direction. The air mass, which is moving during the approach, will have displaced this point exactly to the landing spot at time of touch down. We demonstrate this in examples shown in Sect. 4. For this purpose, the trajectory calculated in the air reference system was transformed segment by segment back into the earth reference system. Our method is thus able to take any complex wind distribution into account. If these are not known from measurements, forecast models or standard gradients can also be used.

Results
The described method is now applied to calculate trajectories for reported flights with emergency landing situations. We assume a simple APM with a fixed EAS for best glide. Configuration changes (deploying flaps or landing gear) were not modeled. n = 100 was chosen as chain length each time. The algorithm was implemented in Java and executed on a tablet with the Android operating system. Running times ranged from 10 to 100 ms. This makes it possible to run the algorithm periodically with an update rate of 10-100 Hz. The trajectory is thus always calculated from the current position of the aircraft.

US Airways Flight 1549
US Airways Flight 1549 [20] is a well-known and intensively studied forced gliding flight. In numerous works, such as [8] and [21], methods based on Dubins paths were used to show that a return flight to La Guardia might have been possible. Our experiments confirmed this. We used historical data from [22] which includes wind data for the position, altitude and time of the incident, see Fig. 11. We assumed a glide ratio E(0) of 17.25 and an EAS of 210 kt.
We calculated trajectories starting from the position and heading given in [8] for the flight 16 s after both engines failed to two accessible runways of the departure airport La Guardia and to a runway of Teterboro Airport to the west. We allowed a maximum bank angle of 45 • . The trajectories found with our methods are shown in Fig. 12. The prevailing Fig. 10 Principle of avoiding obstacles (not to scale). a Trajectory without obstacle. b Trajectory disturbed by obstacle. c Trajectory was "repaired" by further FABRIK runs wind from northwest has little effect on the trajectories to runways 13 and 22 of La Guardia airport. The tailwind component predominates on one part of the trajectory, and the headwind component on another. The effects almost cancel each other out. Ignoring the head wind towards Teterboro Airport's runway 24, it is within gliding distance. If you consider the actual, adverse wind conditions, this runway is only barely reachable. A 1/4 NM straight final was modeled for the approach to La Guardia's runway 13 With La Guardia runway 22, a relatively long straight final approach arises by itself. Forcing a straight final for the approach to runway 24 at Teterboro would have rendered that RWY unreachable.

Longest glide: Air Transat Flight 236
A long gliding flight from cruise to touch down was Air Transat Flight 236, also known as ≫The Azore Glider≪ , [23]. On its way from Canada to Portugal both engines of the Air Transat Airbus A330 failed about 120 km northeast of the Azores Island of Terceira at an altitude of 34,500 ft.
This flight is used to illustrate further advantages of our method. While aerodynamic properties of an aircraft essentially correlate with the EAS, other parameters regarding to the air reference system, such as the curve radius at a given bank angle, depend solely on the TAS. While the difference between EAS and TAS can be neglected for flights at low altitudes (e.g. US Airways Flight 1549), it is significant for flights at usual airliner cruising altitude. Table 1 lists the relevant parameters for selected segments of the trajectory modeled as a kinematoid chain. The altitude differences per segment remain constant by definition, as well as the EAS of approximately 108.0 m/s ( ≈ 210kt ) specified by the APM. In the early phase of the gliding flight at an altitude of more than 10,000 m, the TAS is far higher than the EAS and thus determines the minimum turning radius at a given bank angle. During the flight, the TAS gets ever closer to the EAS. With our method, it is possible to model air density related differences between EAS and TAS without any problems.
Both wind speed and direction vary greatly in this example (Fig. 13). Since the flight time of each chain segment can be calculated in addition to the altitude and speed, the influence of the wind can be modeled precisely-regardless of the complexity of the lateral, vertical and temporal distribution of wind speed and direction.
As described above, a wind offset vector and thus a landing spot in the air reference system (virtual RWY in Fig. 14) is calculated from the individual chain segments and the wind distribution. For our simulation, the maximum bank angle was limited to 5 • . From a distance of 8 NM from the landing spot, the bank angle was limited to 0 • and thus a straight final approach was forced.
According to the Accident Investigation Report [23], the aircraft arrived much too high at the beginning of an 8 NM final at an altitude of approximately 13,000 ft. After the engine failure, the pilot flew on a direct course to the runway, probably to be sure to reach it. This tendency seems to be typical in such situations. (See [24].) A full circle was flown to reduce excess height. Since full circles only allow height to be reduced in discrete (and relatively large) steps, they are only conditionally suitable for this task. According to the investigation reports, the remaining altitude was reduced by S-curves. The presented method has the advantage that it can calculate a trajectory that has exactly the curvature that ensures arrival at the threshold with the correct height and heading.

Emergency landing scenario in complex terrain
We would like to demonstrate the ability of our method to keep a trajectory away from a specific area using the example of an emergency landing situation for a single-engine, general aviation airplane after an engine failure. It shows such a situation in the East Tyrolean Alps. An engine failure occurs near an airfield. The APM used assumes a speed of 73 kt and a glide ratio E(0) of 11.
Without obstacle avoidance, a mountain-intersecting trajectory would be generated. The mountain is modeled as a cylindrical obstacle and the trajectory is deformed using the methods described in Sect. 3.3. In the example, a single cylindrical obstacle is sufficient. However, more complex  Figs. 15 and 16.) However, the detour around the obstacle does not lead to an extension of the trajectory in a glide flight. On the contrary, the more curved curve leads to a higher maximum and average bank angle. This reduces the average glide ratio and leads to a slight temporal and spatial shortening of the trajectory, see Table 2. This corresponds to the behavior of a real flight and need not be explicitly modeled. It rather is the result of a realistic modeling of flight behavior, which is inherent to our method.

Conclusions
The technique we present allows an efficient calculation of glide trajectories. It directly addresses the actual problem, namely the calculation of a trajectory between a starting state and an end state. Conventional methods based on Dubins paths, on the other hand, first calculate the shortest path and then-if necessary-extend it or reduce glide properties of the aircraft by changing its configuration. Our method realistically models the increasing air density during descent and the associated change in TAS if EAS stays constant. Wind vectors that vary in time and space can be precisely considered. The complexity of the wind distribution does not affect the complexity of our method. This makes it suitable for the calculation of forced gliding flights in emergency situations (Azure Glider, Gimly Glider, US 1549). Residual inaccuracies can be arbitrarily reduced for practical purposes by increasing the number n of chain segments. To reduce the residual error even for small n, we introduced circular arc chain segments instead of straight ones. Various approaches to influence the shape of the generated trajectories were presented. It is thus possible to avoid certain areas, given sufficient length of the trajectory. These can be both obstacles, which must be avoided in any case, as well as avoidance areas for noise reduction. In the case of the latter, the deforming influence on the trajectory is only allowed to such an extent that other, higher-priority conditions are still met.
If a trajectory must be calculated to a landing spot that is close to the starting state of the trajectory, it is possible-as already suggested in other publications-to reduce part of the altitude by means of full circles. Here, too, our method proves to be advantageous. Since the air density and thus the relations between EAS and TAS are calculated segment by segment, full circles can be modeled realistically over a large range of altitudes. While keeping bank angle and EAS constant, due to the decreasing TAS these full circles are deformed into a spiral that narrows downwards.
The method is easy to implement as an online algorithm and runs on commercial off-the-shelf hardware. 4

Fig. 16
Trajectory of an emergency landing scenario to LOKL airfield (East Tyrol, Austria) with and without obstacle avoidance, obstacle model shown

Future work
Our current and future work concentrates on proving the convergence of the method and identifying its limits. Furthermore, it should be proved that at least one trajectory can be efficiently calculated within a certain time for each combination of initial and final state for which one or more glide trajectories exists. We will extend the ability to avoid obstacles to more complex situations. In addition, we would like to extend the method to flight trajectories for powered flights. In the work presented, we did not change the configuration (landing gear, flaps, spoilers, etc.) of the aircraft during approach. This corresponds in part to emergency situations, in which the hydraulic system fails in addition to the engines. However, for intentionally initiated gliding flights, for example as part of CDO, it is necessary to be able to model different configurations and switch between them-according to the operating procedures of the aircraft. We will extend our method by this property.
So far, we have assumed that the flight properties of the modeled aircraft, in particular the glide ratio, is precisely known. In practice this will not always be the case. If the actual glide ratio deviates from the model, the trajectory actually flown will also deviate from the modeled one. Since our algorithm periodically recalculates the trajectory from the current position during the flight, the part of the trajectory that is still in front of the aircraft-as well as the time of landing-is constantly shifting. In future work, we want to derive correction parameters for the APM from deviations between the predicted and actual flight path.
Funding Open Access funding enabled and organized by Projekt DEAL. No funding was received to assist with the preparation of this manuscript.
Data availability statement Not applicable.

Conflict of Interest
The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.