Equations used in flight code

Fitting a radius

In the theory section, an equation was produced for the change in ground track angle or track angle,

dχ/dt = (1 / R_E) { ω - (A + u) cos(θ) sin(φ) }.

However, knowing the track angle at time t and time t+1, still requires these points to be connected. This can be done by fitting curve of radius R, where you can show that the arc length, d(Arc), is given by R dΧ (R is dimensionless, because our Mercator space is dimensionless). And because (R/S) d(Arc) = G dt, where G is the ground speed then

S/R = (1/G) (dχ/dt)

and so the equation for (dχ/dt) can be replaced by an equation for (1/R), which is largely what is done in the code. The case where 1/R is very small, when there's no change in heading, is dealt with later.

The advantage of knowing R is that is makes connecting two points easier. Because if two points are connected by a radius R, and the first point is on bearing χ, while the second point is on bearing χ + Δχ then it's possible to draw a diagram (on p61 of my personal notes) and show that

Δx = R (sin (χ + Δχ) - sin(χ))

Δy = R (cos(χ) - cos(χ + Δχ))

where Δx and Δy are the change in Mercator space. The distance travelled, or the length of the arc, is R Δχ in Mercator space - or Δχ R/S on the Earth's surface.

Forward in time equation

It should be noted that the conditions at time t are used to determine the path taken from time t until time t+1, which is a forward in time equation.

If we allowed a step to be very large, it would be possible for a path approaching, say, a jet stream to take a measurement t just before the jet stream and a measurement t+1 just after the jet stream - and so would completely miss the jet stream all together. So it is important that the step shouldn't be too large.

Size of step

It would seem sensible to take a step similar to the grid spacing of the input wind data. Much smaller than this and the benefits are likely to be fairly negligible, because neightbouring timesteps will just be using slightly different proportions of the same grid points.

The input wind data has grid points spaced by Δφ and Δλ, so we take an angle Δ which is representative of smallest of these two angles, given by

1/Δ = 1/Δφ + 1/Δλ.

And so the distance between the two points, as the crow flies, is R_E Δ. The distance travelled by the plane is slightly longer because it travels in an arc rather than a straight line. In Mercator space, if the plane travels an arc R Δχ, it travels 2R sin(Δχ/2) as the crow flies. Hence

R_E Δ = (2R/S) sin(Δχ/2)

where

Δχ = 2 sin^-1 [Δ / (2R cos(φ))].

And the time taken to complete this is the distance travelled, (R ΔΧ / S) by the ground speed G, so

Δt = R ΔΧ / (S G).

No change in heading (1/R = 0)

When there is no change in heading, Δχ will be zero and R is ∞. But this will actually cover any case where 1/R is sufficiently small for the code to risk a floating point exception, when dividing by it. In this case the distance travelled in Mercator space is S R_E Δ and so

Δx = Δ cos(χ) / cos(φ),

Δy = Δ sin(χ) / cos(φ).

Each step the progression in time is the distance travelled divided by the ground speed

Δt = R_E Δ / G.

Ground speed and nose heading

Cross-winds mean that the air velocity, A, is not in the same direction as the velocity with respect to the ground - the ground speed vector, G. Let u_track be wind vector in the direction of G, and v_track is the cross-wind (the component at right angles and left of u_track). Note that this is different to u and v used before because they are reference by the nose heading.

The direction of the plane, given by G, must be the direction of the sum of its air velocity, A, and the cross-wind velocity, v_track. By definition v_track is at right angles to G and so the length of the component A + v_track in the direction of G is given by Pythagoras' theorem as |A + v_track| is √ (|A|² - (v_track)²|). Adding the tailwind to this gives the full length of the ground speed vector

|G| = |u_track| + √ (|A|² - (v_track)²|).

Be removing the wind vector, w = u_track + v_track, from the ground speed vector, we can determine the direction of A and hence determine the nose heading,

θ = tan^-1(A_y/A_x) = tan^-1[(|G|sin(χ) - V) / (|G|cos(χ) - U) ]

where A = A_xi + A_yj, and w = Ui + Vj.

Determining U, V and ω

For a given point, it is unlikely that a grid point will exactly match the location, and so the nearest four points are averaged to calculate U and V. Nearest points are given a higher weighting in the average. This is done by calculating the east and north weighting factors,

w_E = (x - x_W) / (x_E - x_W)   &   w_N = (y - y_S) / (y_N - y_S)

where the grid points, starting with the north-west point and moving clockwire, are located at (x_W, y_N), (x_E, y_N), (x_E, y_S) and (x_W, y_S), and the given point is at (x, y). And the subscript N indicates the nearest latitude north of given point, subscript E the nearest longitude East of grid point and so on for S and W. U at (x,y) is then calculated as

U = (1 - w_E) w_N U(x_W, y_N) + w_E w_N U(x_E, y_N) + w_E (1 - w_N) U(x_E, y_S) + (1 - w_E) (1 - w_N) U(x_W, y_S)

and V is calculated in a similar way. The velocity gradient dU/dy is determined from calculated the difference in U between neighbouring latitude levels and so is stored on a levels half between the actual latitude y levels, known as a half level. Similarly, dV/dx are stored as level half way between the longitude x levels. However, the method for determining their value at (x,y) is the same as used for U above, but obviously the weights will be different as the variables are located at different places.