Theory for flight code
The path of least time
The title of Sawyer (1949) is `Theoretical Aspects of Pressure-Patten flying', and is true to its title. It states on p42 that over flat terrain that the equation for the path of least time is given by
dθ'/dt = ∂u/∂n
where θ' is the nose heading of the plane, and is measured in a clockwise direction from North; t is time; u is the tailwind (in same direction as θ'); and n is at right angles and left to the direction given by the nose heading. Note that nose heading is not the direction the plane actual travels - which is the ground track angle χ - but the direction the plane points. It's the direction the plane would travel in the absence of a cross-wind.
The curvature term
Equation in Sawyer (1949)
Sawyer (1949) goes on to say that when the curvature of the Earth is included the equation becomes
dθ'/dt = ∂u/∂n + ((A + u) / S) (∂S/∂n)
where S is the scale factor of the map. A vigorous derivation of the above equation is given in `Washington, Headquarters, Air Weather Service; The minimal flight path. Washington, D.C., 1947', which I've not seen. On pages 42-44 of Sawyer (1949) a pictural derivation of the equation is given.
Mercator projection
I'm going to do even less than that, but say that the equations are firstly derived on a flat surface with an x and y axis at right angles to one another - before adding the curvature correction term. On the Earth this is equivalent to deriving the equations on a Mercator projection where the x and y axis are at right angles to one another and angles between the Mercator projection and the real lat-lon space are preserved, where
x=λ
y=ln [ tan(45 + φ/2)]
where λ is the longitude, x has been set to zero when λ is zero and x and y have been non-dimensionalised by the radius of the Earth, RE (no effort in the theory is made to allow for the Earth not being a perfect sphere).
Hence, for our Mercator projection the scale factor is given by
S = DistMerc / DistReal = 1/ (RE cos(φ))
where DistMerc is any distance on the Mercator projection and DistReal is its real distance on Earth. So at the equator this is just 1/RE, but as we move to the poles then any distance on the Earth becomes much larger and so S must increase. See the Wikipedia page for more details on S in a Mercator projection.
How angles are defined
For work in this area, there seem to be variey of ways of defining angles so that they can either increase in a clockwise or anti-clockwise direction, and they start from diffent points on the compass. When working in Mercator space it's natural to define angles as zero which point along the positive x-axis and to increase in a anti-clockwise direction. I'm going to define all my angles like this. This means that my nose heading θ is effectively zero when pointing East and increases in an anti-clockwise direction (90° is North, 180° is West and 270° is South). And hence my dθ has the opposite sign to the dθ' used by Sawyer (1949).
Deriving curvature term in Mercator space
Having defined our nose heading angle, we can match our x and y axis with the coordinates used by Sawyer above,
| dx = ds cos(θ) - dn sin(θ) | ds = dx cos(θ) + dy sin(θ) | |
| dy = ds sin(θ) + dn cos(θ) | dn = -dx sin(θ) + dy cos(θ) |
where s is in the direction of the nose heading. Hence
(1/S)∂S/∂n = (1/S) (∂y/∂n) (dS/dy) = (1/S) cos(θ) (dS/dφ) (dφ/dy)
where dφ/dy = 1/S, so that
(1/S)∂S/∂n = (1/S2)cos(θ) (dS/dφ) = (1/S2) cos(θ) (sin(φ) / cos2(θ)) = cos(θ) sin(φ)
Hence the equation for the heading becomes
dθ/dt = - (1 / RE) { ∂u/∂n + (A + u) cos(θ) sin(φ) }
where the source terms are now negative because dθ = - dθ'.
Change in ground track angle
The equation above gives the change in nose heading, but what is really required is the change in ground track angle. As shown in Lunnon and Marklow (1992), the ground track angle, Χ is the sum of the nose heading, θ and drift angle, γ, where the drift angle is given by
sin γ = v / A
where v is speed in the direction n and A is the air speed. Given that v << A, γ ≈ v / A, and
dγ/dt = (1/A) (dv/dt) = (1/A) { (ds/dt) ∂v/∂s + (dn/dt) ∂v/∂n } = (1/ (A * RE) { (A + u) ∂v/∂s + v ∂v/∂n }
and given that u,v << A then
dγ/dt ≈ (1 / RE) ∂v/∂s.
And hence the equation for rate of change of ground track angle for the path of least time is given by
dχ/dt = (1 / RE) { -∂u/∂n + ∂v/∂s - (A + u) cos(θ) sin(φ) }
However, as vorticity ω = -∂u/∂n + ∂v/∂s the equation becomes
dχ/dt = (1 / RE) { ω - (A + u) cos(θ) sin(φ) }
and voriticity is a coordinate invariant quantity and so it is unneccessary to rotate it. Easiest to calculate ω as -∂U/∂y + ∂V/∂x where x and y are the Mercator coordinates , and U and V are the velocities in the x and y direction.