5.1 Introduction

Overview

To recap what we have done so far:

• we have fixed a horizontal route; in our illustration this was ['WETSI', 'DAYVU', 'CRACK', 'TANIE', 'ICADI', 'BATTY', 'AUGEY', 'KECKI', '60N50', 'NOWEL', 'OLLEO']

• we have created aircraft performance lookup functions for climb, cruise and descent 

The next remaining things we have to figure out are:

• fuel used (per leg and total)

• time required (per leg and total)

• TOC (top of climb) and TOD (top of descent) positions 

Spoiler alert: we are going to solve the remaining values by computing the trajectory backwards from the departure waypoint all the way to the arrival waypoint. In this chapter we are going to set up the mechanism to compute this trajectory. We will compute the trajectory in 10 second increments. We compute in reverse because we know the desired fuel level at the arrival point, but not the departure point. 

Note: we haven't forgotten about the weather ingredient, but we are saving that for a later chapter. When we set up the trajectory computation we are going to leave room for a wind vector input. 

Computation Strategy

Here is a brief outline of how the computation will work. 

Step 1

We start by creating an initial state vector. This is a collection of all the relevant aircraft attributes at the initial point of the computation, i.e., the arrival waypoint. By attributes we mean things like: 

• latitude, longitude, and altitude

• horizontal and vertical speed

• fuel level and aircraft weight 

• distance and direction to next and final waypoints

• phase (climb, cruise, or descent)

Note: since we are computing the trajectory in reverse, the initial value for the phase attribute is 'descent'. 

Step 2: 

We estimate the position of TOC using our estimate_climb_dist_m function (created in last chapter, see 5.5)

Step 3: 

We create a function called advance_state_vector, which has:

• inputs: a state vector, and a step time (10 seconds)

• output: an updated state vector

As the function name suggests, its job calculate the updated state vector attributes of the aircraft after one timestep has elapsed. 

Note 1: We are computing in reverse, so even though we say "advance", we are really going backwards in time. 

Note 2: This function has logic checks to update the phase field:

• it will update to 'cruise' when 30,000 ft MSL is reached

• it will update to 'climb' when the trajectory has reached the estimated climb distance 

Step 4: 

We start a loop to advance the state vector in 10 second increments until the end condition is met. The end condition is: the state vector has reached within 4 km of the last waypoint, i.e., the departure waypoint. Each new state vector will get added to a list of vectors, forming a trajectory. 

Step 5: 

The trajectory (which is just a list of state vectors) will be sent to a function called assess_trajectory(). As the name suggests, its job is to judge whether or not the computed trajectory is acceptable or not. The judgement is based on the altitude of the final state vector: if it is within the tolerance (1000 ft), the trajectory is accepted. If not, the function will issue a revised climb distance estimate and ask the loop to recompute the last portion of the trajectory. The revised estimate is based on the difference between the trajectory's final altitude and the target altitude. 

Visually, the process we describe in Step 5 looks like this:

Tolerances

Where does the 4 km horizontal tolerance come from? 

A plane traveling at 240 knots is traveling about 4 nautical miles per minute, or 2 nautical miles per 30 seconds; this is about 4 km in 3 iterations, since each iteration is 10 seconds; we feel is a reasonable tolerance.

Where does the 1000 ft vertical tolerance come from? 

At the end of the trajectory, we are expecting vertical speeds of about 4000 ft / min; this works out to about 600 feet per iteration since each iteration is 10 seconds; we feel this is a reasonable tolerance.

Details on the state vector computation

In this next section I will show the code for the trajectory computation, which will include the code that updates the state vector in 10 second increments. Before showing the code I will give a few explanatory notes on one tricky aspect of this state vector calculation. 

We are throwing around words like "computation", "trajectory", and "state vector". To clarify: this is not a physics-based trajectory simulation in which we try to model all the forces, moments, mass, acceleration, etc. of the aircraft. We don't have enough information to construct such a model. For the most part, we are just looking up values with the aircraft performance functions we created earlier and applying them in 10 step increments. 

However, there is one tricky part which deserves some explanation prior to presenting the code: the relationship between airspeed, groundspeed and wind speed. To untangle this relationship let us create three velocity vectors:

• a: aircraft velocity 

• b: wind velocity

• c: groundspeed vector 

The relationship between them is: 

a + b = c 

If wind (b) is zero, then the aircraft velocity and groundspeed vector are the same (a = c) 

If wind is nonzero, then we must find two unknowns:

• the direction of vector a; i.e., we need to find the heading the aircraft should fly in order to correct for the wind. In aviation this is called the wind correction angle. 

• the magnitude of vector c, i.e. the resultant groundspeed 

Note that we know the direction of the groundspeed because we assume it is always directed towards the next waypoint.

Graphically, the relationship between the three vectors looks like this:

We can find our unknowns (the direction of vector a and the magnitude of vector c) by taking the steps that follow. 

Step 1:

First, we find θ. Considering the dot product of two vectors, we can say that:

b ⋅ c = |b||c| cos θ

We don't have the magnitude of c, but since we know the direction, we can just use the unit vector ĉ instead:

b ⋅ ĉ = |b||ĉ| cos θ 

The magnitude of a unit vector is one, so we can rewrite this as:

b ⋅ ĉ = |b|cos θ 

Step 2:

From the Law of Cosines, for our triangle we can write:

a2 = b2 + c2 - 2bc cosθ

We need to rearrange this to solve for c. We can't really isolate c but we can rearrange it into the form of the quadratic equation: 

c2 − 2b cosθ c +  b2 − a2 = 0

Here the coefficients are:

• coefficient_1 = 1

• coefficient_2 = 2b cosθ

• coefficient_3 =  b2 − a2

Using the quadratic formula we can find the roots of this equation, i.e. find the value of c. When working with the quadratic equation there is a fancy thing called the discriminant, which can tell us whether the quadratic equation will have one root, two roots, or zero roots. Based on our triangle, we expect to always have one positive root for our values of wind and airspeed. 

Once we have the magnitude of c, we can create vector c by multiplying the unit vector by c: 

c = c ĉ 

Then solve for the direction of a by vector subtraction:

a = c −  b

In the next section we will show how we coded this process into Python.