Introduction

In this lab we use Maple to demonstrate an important application of linear algebra. When an airplane flys (nonstop) from Chicago to Timbuktu it takes the shortest route between the cities, which is a great circle route. At each point along the route, the pilot needs to know what compass heading to use. We use linear algebra to compute this information.

Getting Started

As usual, we have to load some packages and define some functions to help us in our computations. First, we load the familiar linear algebra and plotting packages.

> with(plots):

> with(linalg):

To make the computations easier to follow, we will compute and display numerical approximations at each step. This is less accurate than the symbolic approach, but we can partially compensate for this by using a high degree of precision in our calculation.

> Digits := 20;

The UnitVector function defined below computes a numerical approximation to the unit vector in the same direction as the vector v. i.e. the direction vector corresponding to v .

> UnitVector := v -> evalm((1/norm(v,2))*v);

The Computation

Here are the latitudes and longitudes of two cities -- Chicago and Timbuktu. Latitudes describe the cities' east-west position on the globe, and are between -Pi/2 and Pi/2 . Longitudes describe north-south position and range between -Pi and Pi . West longitudes and south latitudes are negative.

> lat1:= 41.5*Pi/180;

> long1 := -87.45*Pi/180;

> lat2:= 16.49*Pi/180;

> long2 := -2.59*Pi/180;

Here are the position vectors of the two cities (where the radius of the earth is assumed to be 1). It is not hard to see how the z coordinate is related to lattitude, but can you see where the formulae for the x - and y -coordinates came from?

> Vector1 := evalf([cos(long1)*cos(lat1), sin(long1)*cos(lat1), sin(lat1)]);

> Vector2 := evalf([cos(long2)*cos(lat2), sin(long2)*cos(lat2), sin(lat2)]);

We use the dot product to compute the angle between the position vectors.

> TheAngle := arccos(dotprod(Vector1,Vector2));

Then we convert to degrees and multiply by 60 to get the distance between cities in nautical miles. (Sixty nautical miles make up one degree of longitude on the globe. A nautical mile is about 15% longer than a regular mile.)

> TheDistance := evalf(60*TheAngle*180/Pi);

The vectors perpendicular to the two position vectors are also perpendicular to the great circle joining the two cities. We compute one of them. This vector points from the center of the earth toward Peru; it is indicated by a black + in the picture below.

> NormalVector := UnitVector(crossprod(Vector2,Vector1));#

Then we find the unit vector pointing due east from the first location:

> East1 := evalf([-sin(long1), cos(long1), 0]);

Computing the angle between these two vectors allows us to compute the angle between the vector "east from Chicago" and the great circle joining the two cities. This tells us the compass heading from the first city to the second along the great circle route, as an angle from due north; we must decide whether the heading is east or west of north. (Why does Vector2 come first in the cross product?)

> Azimuth:= evalf(arccos(dotprod(East1, NormalVector))*180/Pi);

We can use the sign of the z -coordinate of NormalVector to decide if we are heading east or west. (Why?)

> if NormalVector[3] < 0 then CompassHeading = Azimuth else CompassHeading = 360 - Azimuth fi;

The information we have just computed tells the pilot what compass heading to use just after takeoff when traveling from Chicago to Timbuktu. Over time, the heading of the plane will have to change, and this computation will be repeated with Chicago's lattitude and longitude replaced by the current lattitude and longitude of the airplane.