Navigation Project

Introduction

A Mercator Projection maps a point on the earth's surface with longititude θ (degrees, positive for east of Greenwich) and latitude λ (degrees, -90<= λ <=90, positive for the northern hemisphere) to a rectangle in such a way that 

(i) curves ("lines") of constant latitude (the "parallels") are mapped to horizontal straight lines, 
(ii) curves of constant longitude ("meridians") are mapped to vertical straight lines, 
(iii) the meridians are mapped to equally spaced vertical lines;
(iv) the map is conformal---ie,  a short (ie, infinitesimal)  line on the earth's surface making an  angle α with the meridian maps to a line making the same angle α with the mapped meridian.

Recall that the usual spherical coords (r, θ, ϕ) use r=radius of earth,  θ = longitude, ϕ = colatitude, so that ϕ = 90-λ. The advantage of cordinates (r, θ, λ) is that the orientation of the pair (θ, λ) is correct (see next diagram) , rather than inverted in the vertical direction.

Show[
        Graphics[Point[{0, 0}], Axes -> True, 
          AxesLabel -> {"θ",  "λ"},  
          PlotRange -> {{-180, 180}, {-90, 90}}]];

Let's agree to use {u,  v} as coordinates on the Mercator map, with 

u = θ
v = h[λ]

where the function h is to be determined by the requirement of conformality.

The following picture shows the essence of conformality and gives the basic relationships linking distance on the sphere and. The graphic is supposed to represent a short directed line segment starting at an arbitrary point on the earth's surface (other than one of the poles), corresponding to a change Δθ (degrees) in longitude and Δλ (degrees) in latitude. The base of the triangle is actually an arc of a circle (the latitude at λ) having radius r Cos[λ Degree] Δθ Degree, and the altitude is an arc of a circle of radius r, hence having length r Δλ Degree. Thus the angle α (degrees) the line segment makes with the positive (easterly) directions satisfies Tan[α Degree]= (Δλ/Δθ ) / Cos[λ Degree].

arc[A_, B_, ra_] := 
      Module[{c, d, 
          u}, 
d = (A + B)/
            2;    
u = {-(B - A)[[2]], (B - 
                    A)[[1]]}/
            Sqrt[(B - A) . (B - 
                    A)] ;   
c = 
          d + Sqrt[
                ra^2 - (B - A) . (B - A)/
                    4] u;    
Circle[c, 
          ra, {ArcTan[Sequence @@ (A - c)], 
            ArcTan[Sequence @@ (B - c)]}]]

A1 = {0, 0}; A2 = {1, 0}; A3 = {1,  .5};
\(\(Show[
        Graphics[{arc[A1, A2, 4], arc[A2, A3, 5], 
            arc[A1, A3, 
              5], 
 (*\(Line[{\ {0, 0}, {1, 
                    0}, {1,  .5}, {0, 
                    0}}]\)\(,\)*) 
			Text[α, {0.2, 
                0.05}], 
			Text[Δs, {0.5, 
                0.29}], 
			Text["r Δλ Degree", {1.1, 0.25}], 
			Text["r Δθ Degree Cos[λ Degree]", {0.5, \(- .05\)}], 
Text["Small triangle on earth surface", { .5,  .6}]		
		}, 
Axes -> False, \ 
          PlotRange -> {{0, 1.3}, {\(- .1\),  .8}}, 
          AspectRatio -> 1
	]];\)\)

\(\(Show[
        Graphics[{Line[{\ {0, 0}, {1, 0}, {1,  .5}, {0, 
                  0}}], 
			Text[α, {0.2, 
                0.05}], 
			Text["Δv", {1.1, 
                0.25}], 
			Text["Δu", {\ 
                0.5, \(- .05\)}], 
Text["Small triangle on Mercator map", { .5,  .6}]			
		}, Axes -> False, \ 
          PlotRange -> {{0, 1.3}, {\(- .1\),  .8}}, 
          AspectRatio -> 1
	]];\)\)

\(
\)

As Δv=h'[λ] Δλ (up to infinitesimals of smaller order) and Δu=Δθ, we conclude from conformality that

(1.1) h'[λ] =1/Cos[λ Degree].

We also conclude fromn the previous picture that 

(1.2) ds=dλ/Sin[α].

Conformality made the Mercator projection useful to early navigators, as they could set a bearing based on drawing a straight line on the map rather than trying to continually change bearing, as is required when travelling a great circle.

Questions

In all that follows, let A, B have coordinates {(Subscript[θ, 0,  ]) (Subscript[λ, 0])}, 
       {(Subscript[θ, 1]),  (Subscript[λ, 1])} respectively (all in degrees.) Take the earth to be a sphere of radius r = 3963 (miles).

(a) Find an explicit expression for v in terms of λ.

(b) Using dot product, set down a Mathematica function greatc[A_,B_] giving the great circle distance between A and B.

(c) Use (1.2) to derive a Mathematica function rhumb[A_,B_] for distance along the shortest "rhumb line" from an arbitrary point A to another arbitrary point B on the earth's surface. ( A rhumb line joining A and B is a curve on the earth corresponding to a straight line segment  from A to B on the Mercator map. There are in general many rhumb lines joining A and B. You may take it that the shortest of the possible rhumb lines is the one that stays in the same hemisphere of longitude. You will also need to distinguish two cases, according to whether A and B lie on the same parallel.)

(d) What is the maximum possible difference between the shortest rhumb line distance and the great circle distance? Under what circumstances is the maximum attained?[This will involve maximizing over several variables. Doing a reasonable eyeball from graphs (withy aptly chosen viewpoints) is acceptable.]

(e) Write Mathematica functions that will, on a map of the world (using WorldPlot), mark points every 500 miles along the trajectories from given points A, B. Print the distances along both trajectories on the map itself. You should use the WorldRotation option to center the map longitudinally half way between A and B, choosing the correct hemisphere for "half way". Be aware that WorldPlot uses coordinates in form lat, long, both measured in minutes.

(f) Use the same sort of construction to devise an animation of the trajectories of two airplanes, one flying the great circle from Sydney to  London, the other the rhumb line. Assume the planes both fly at 500 miles/hour, and update the animation every hour.

Useful facts

A rhumb line on the map from {u0, v0} to {u1,v1} can be written in  the parametric form 

 u= (u1-u0) t + u0,
 v = (v1-v0) t + v0,
 as the parameter t varies from 0 to 1.

Relation between {x,y,z} recangular coordinates and {r,θ,λ}: (not quite standard spherical polar coords.)

x = r Cos[λ] Cos[θ]
y = r Cos[λ] Sin[θ]
z = r Sin[λ]

Great Circle

r  =  3963

Define first a function to transform a point on earth surface with coords 
A={θ, λ} to 3 dimensional rectangular coordinates in the 
form of a list {x,y,z} so we can handle dot product to find angle 
 between two points.  Then recall that Dot[v,w] = |v| |w| Cos[α] where 
 α is the angle between v and w, and use it to calculate the length of a great circle path.
 Check that it works on antipodal points: result should be half of earth circumference, π r, approx 12500.

Rhumb Line

Using the arc length element given above, we can find the rhumb line distance.  We first need to find h[λ]. We can do this by integrating, since we already know h'[λ].  Caution: be careful with converting between degrees and radians.

If we want to go from A = {(Subscript[θ, 0]), (Subscript[λ, 0])} to B = {(Subscript[θ, 1]), (Subscript[λ, 1])} along a rhumb line, we could take the line on the Mercator mat having parametric form:

u = (Subscript[u, 0])+ ((Subscript[u, 1])- (Subscript[u, 0])) t
	v = (Subscript[v, 0])+ ((Subscript[v, 1])-(Subscript[v, 0]))t
	
with 0 ≤ t ≤ 1.

This particular rhumb line may not however give the shortest of the possible rhumb lines from A to B. (Consider for example A={179,0}, B={-179,0}.) To get the shortest rhumb line, we should allow other possible representations of B, adding multiples of 360 (degrees) to its longitude. We will take it as evident that the shortest rhumb line will correspond to the choice of representation of (Subscript[θ, 1])  that is within 180 degrees of (Subscript[θ, 0]) (In other words,  the shortest rhumb line is the one that remains in the same hemisphere of longitude as the starting point.) For example, this would mean that the B above should be replaced by {181,0}. Thus, we have a different parametrization of the paths, expressed as follows:

u = (Subscript[u, 0])+ c t
where c = (Subscript[u, 1]) - (Subscript[u, 0]) = (Subscript[θ, 1]) - (Subscript[θ, 0]) if   | (Subscript[θ, 1]) - (Subscript[θ, 0])|  ≤ 180, else c = (Subscript[θ, 1]) - (Subscript[θ, 0]) + 360n where n is chosen so that -180 < c ≤ 180.  There is a compact expresion for c, namely

	c = Mod [((Subscript[θ, 1]) - (Subscript[θ, 0]) + 180), 360] - 180.  To check, consider:

Mod[-352,360]

In computing the rhumb line distance, we can use the fact that if the rhumb line makes an angle α with East, then 

s = r Δλ Degree / Sin[α] 

gives us the arclength travelled as long as α is non-zero. On the other hand, the case α = 0 is easy, since in that case we move along a parallel and so

s = \(TraditionalForm
      \`r\  . change\ in\ longitude\ in\ 
        radians . \ Cos[λ\_0\ Degree]\).

Now find the maximum difference between the rhumb line distance and the great circle distance.  You may want to use Plot[] and Plot3D[] to estimate, choosing an advantageous viewpoint.

Trajectories

We now need to work out the great circle trajectory on the Mercator map.Given points A, B, the unit vectors a, b in their respective  directions may not be perpendicular, but if they are, the circle from A to B is given in terms of angle w as parameter by

c=Cos[w] a + Sin[w] b

If a and b are not perpendicular, we form the unit vector bb perp to a and in the same plane as a and b, and pointing to the same side as b, by

bb = (a ×b) × a

Then the great circle starting at A and going through B is given by 

c[w_] := Cos[w] a + Sin[w] bb

with w the angle at the center.

We need to convert from polar form in the lists A, B to rectangular form in a, b in order to carry this out.  We have to watch for cases where A and B are parallel or antiparallel.

To convert back to {longitude, latitude}form, we need functions that take unit vector {x,y,z} form back to long, lat.

We can now draw the parametric curve on the Mercator map corresponding to a great circle from given points A, B.

Write a pattern that makes a list of points on the great circle trajectory and a pattern that makes a list of points from the rhumb line trajectory.

We'll use this to show position on a Mercator map of the world. For this, use WorldPlot, where the units are minutes rather than degrees, requiring us to rescale everything by 60. We also need to rotate the worldplot using the WorldRotation option to make the gc and the rhumb line paths from A to B centered with respect to longitude. Handle this by WorldRotation-> {0,0,- (A[[1]]+deltalong[A,B]/2)}. The the center of this map is now the new origin. Letting cc=A[[1]]+deltalong[A,B]/2, the new coords of a point with latitude, longitude (θ,λ) (degrees) are {60 (θ + cc), 60 h[λ]}

Now show the two different trajectories between pairs of world cities.