Documentation for flight code

Rotating great circle route to equator (deriving α)

The location of any two points on a perfect sphere can be used to define a plane which contains both points and the centre of the sphere. Here the two points are the origin and destination and the Earth is approximated as a perfect sphere. The normal to this plane can be determined by simply taking the cross product of the vector to the origin and the vector to the destination. This plane will be referred to at the great circle plane, because it will also contain every point on the great circle route from orgin to destination. And I'll define the great circle disk as the part of this plane which is contained within the sphere.

Let α be the angle the great circle plane make with the plane through the equator. It will also be the highest latitude, and minus the lowesst latitude, on the great circle disk. The great circle plane can be moved to the Equator by simply rotating it by angle α, about the line where the great circle plane and the plane of the equator intercept.

Angle α can be determined by using the cosine rule and taking the dot product of the normal to the plane, which is cross product of the origin vector with destination vector, and the vector point directly north. Let the origin be defined by latitude, φ₁, and longitude, λ, and destination by (φ₂, λ₂), so that

cos(α) = cos(φ₁) cos(φ₂) cos(λ₂ - λ₁) / sin (γ)

where γ is the angle between the origin and destination and is given by

cos(γ) = cos(φ₁) cos(φ₂) cos(λ₂ - λ₁) + sin(φ₁) sin(φ₂).

Rotating about poles so that zero longitude coincides with highest/lowest part of great circle disk (deriving β)

Having got the great circle plane to the equator it's convient to rotate it about the poles, by angle β, so that zero longitude conicides with either what was the highest or lowest part of the great circle disk. It's fairly unimportant whether the highest or lowest part of the great circle disk is chosen, but I've chosen to take the one which avoids the route crossing 0/360 degrees (it is slightly simpler if there isn't a step change in the longitude of the flight paths).

Let x, y and z be cartesian points in our the normal lat-lon reference frame, which are non-dimensionalised by radius of Earth. Let x', y' and z' be similar to this coordinate sytem but defined for the rotated coordinate system so that

x' = x cos(α) cos(β) + y cos(α) sin(β) + z sin(α)

y' = -x sin(β) + y cos(β)

z' = -x sin(α) cos(β) - y sin(α) sin(β) + z cos(α)

and for reference (not actually needed in this section)

x = x' cos(α) cos(β) - y' sin(β) - z' sin(α) cos(β)

y = x' cos(α) sin(β) + y' cos(β) - z' sin(α) sin(β)

z = x' sin(α) + z' cos(α)

Let φ and λ be the latitude and longitude in the non-rotated frame. And define: x=cos(φ)cos(λ), y=cos(φ)sin(λ) and z=sin(φ). By definition z' must be zero at both the origin and destination, so

-sin(α) cos(β) cos(φ₁) cos(λ₁) - sin(α) sin(β) cos(φ₁) sin(λ₁) + cos(α) sin(φ₁) = 0

-sin(α) cos(β) cos(φ₂) cos(λ₁) - sin(α) sin(β) cos(φ₂) sin(λ₂) + cos(α) sin(φ₂) = 0

which allows a solution for tan(α) to be derived and an equation for β

tan(β) = [cos(φ₁) sin(φ₂) cos(λ₁) - sin(φ₁) cos(φ₂) cos(λ₂)] / [sin(φ₁) cos(φ₂) sin(λ₂) - cos(φ₁) sin(φ₂) sin(λ₁)]

To avoid the route crossing the 0/360 line, we determine if the route is predominantly above or below the equator. This can be done by simply added the height of the origin and the destination together, (z₁ + z₂), which will be positive if route is predominantly above the equator and negative otherwise. If the route is predominantly above the equator, the great circle route will not pass through the minimum point of great circle disk and so λ' is set to zero at the lowest great circle disk height. Otherwise λ' is set to zero at the highest great circle disk height.

Tilt

Even after moving the great circle disk to the Equator for the Mercator projection, the wind speeds which are read in from file will still be in the standard lat-lon coordinates system. So we need to know how the standard lat-lon coorindates system relates to our rotated coordinate system. The unit vectors for the standard lat-lon coordinates system are given by

i_r = cos(φ) cos(λ) i + cos(φ) sin(λ) j + sin(φ) k

i_λ = - sin(λ) i + cos(λ) j

i_φ = -sin(φ) cos(λ) i - sin(φ) sin(λ) j + cos(φ)k

where i, j and k are the standard unit vectors for the x, y and z axies, and i_r * i_λ = i_φ, i_λ * i_φ = i_r and i_φ * i_r = i_λ.

Before defining the unit vectors in the rotated coordinate system is useful to define three unit vectors which will simplify the unit vectors equations, which are

a = cos(α) cos(β) i + cos(α) sin(β) j + sin(α) k

b = -sin(β) i + cos(β) j

c = -sin(α) cos(β) i - sin(α) sin(β) j + cos(α) k

where it can be shown that a * b = c, b * c = a and c * a = b. Deriving the unit vectors for the r, λ' and φ' and using the vectors above produces the equations

i_r = cos(φ') cos(λ') a + cos(φ') sin(λ') b + sin(φ') c

i_λ' = - sin(λ') a + cos(λ') b

i_φ' = -sin(φ') cos(λ') a - sin(φ') sin(λ') b + cos(φ') c

where i_r * i_λ' = i_φ', i_λ' * i_φ' = i_r and i_φ' * i_r = i_λ'. These equations have a similar form the equations for i_r, i_λ and i_φ shown above, but i, j and k have been replaced with the vectors a, b and c - which contain the dependendence on α and β, so that α and β don't explicitly appear in the equations above.

For a given location the plane normal to i_r will contain all the other unit vectors, i_λ, i_φ, i_λ' and i_φ'. The tilt is the angle, T, that the coordinate system (φ',λ') makes with the coordinate sytem (φ,λ), and so from the dot product,

cos(T) = i_λ . i_λ' = i_φ . i_φ'.

And

i_λ . i_λ' = sin(λ) sin(λ') a.i - cos(λ) sin(λ') a.j - sin(λ) cos(λ') b.i + cos(λ) cos(λ') b.j

= cos(α) sin(λ') sin(λ - β) + cos(λ') cos(λ - β)

where between the lowest part and highest part of great circle disk the tilt of the rotated coordinates, (φ',λ'), with respect with standard coordinates, (φ,λ), will be positive, and negative otherwise.

It's possible to derive the equation above from the equation for i_φ . i_φ', but this is far more complicated because of the involement of φ and φ'. It is much easier to show that it is equal to the same angle,

i_φ . i_φ' = i_φ . (i_r * i_λ') = i_λ' . (i_φ * i_r) = i_λ' . i_λ.

For positive α, the tilt T is negative (the great circle route is moving south) between 0° and 180°, and is positive (the great circle route is moving north) between 180° and 360°. And the opposite is true if α is negative.

Having derivied the tilt angle it's simple to switch between angle in the standard coordinate and our rotated coordinates

(angle in φ and λ coordinates) = (angle in φ' and λ' coordinates) + T.