Great-circle lines in Equirectangular projection

Question

Just to check that I'm on the right track:

Are all great-circles on the sphere and in equirectangular projection (i.e. latitude, longitude pairs) either:

1. meridians (i.e. going pole to pole)
2. of the form tan latitude = sin360(longitude + rotation) * amplitude + offset

(with additional restrictions on the offset/amplitude combinations -- obviously, all great circle paths with an amplitude of 0 also have offset 0 - the equator).

Or are there great-circle paths that do not fit into this scheme (again, only in a longitude-latitude coordinate system, not on other map projections).

Note: I added the tan above after posting the question, in a reaction to whubers excellent reply. It turns out that offset then always is 0.

Answer 1

Although geodesics do look a little like sine waves in some projections, the formula is incorrect.

Here is one geodesic in an Equirectangular projection. Clearly it is not a sine wave:

(The background image is taken from http://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/Equirectangular-projection.jpg/800px-Equirectangular-projection.jpg.)

Because all Equirectangular projections are affine transformations of this one (where the x-coordinate is the longitude and the y-coordinate is the latitude), and affine transformations of sine waves are still sine waves, we cannot expect any geodesics in any form of the Equirectangular projection to be sine waves (except for the Equator, which plots as a horizontal line). So let's start at the beginning and work out the correct formula.

Let the equation of such a geodesic be in the form

latitude = f(longitude)

for a function f to be found. (This approach has already given up on the meridians, which cannot be written in such a form, but otherwise is fully general.) Converting to 3D Cartesian coordinates (x,y,z) gives

x = cos(l) cos(f(l))
y = sin(l) cos(f(l))
z = sin(f(l))

where l is the longitude and a unit radius is assumed (without any loss of generality). Since geodesics on the sphere are intersections with planes (passing through its center), there must exist a constant vector (a,b,c)--which is directed between the poles of the geodesic--for which

a x + b y + c z = 0

no matter what the value of l might be. Solving for f(l) gives

f(l) = ArcTan(-(a cos(l) + b sin(l)) / c)

provided c is nonzero. Evidently, when c approaches 0, we obtain in the limit a pair of meridians differing by 180 degrees--precisely the geodesics we abandoned at the outset. So all is good. By the way, despite appearances, this uses just two parameters equal to a/c and b/c.

Note that all geodesics can be rotated until they cross the equator at zero degrees longitude. This indicates that f(l) can be written in terms of f0(l-l0) where l0 is the longitude of the equatorial crossing and f0 is the expression of a geodesic crossing at the Prime Meridian. From this we obtain the equivalent formula

f(l) = ArcTan(gamma * sin(l - l0))

where -180 <= l0 < 180 degrees is the longitude of the equatorial crossing (as the geodesic enters the Northern Hemisphere when traveling eastwards) and gamma is a positive real number. This does not include the meridian pairs. When gamma = 0 it designates the Equator with a starting point at longitude l0; we may always take l0 = 0 in that case if we wish a unique parameterization. There are still just two parameters, given by l0 and gamma this time.

---

Mathematica 8.0 was used to create the image. In fact, it created a "dynamic manipulation" in which the vector (a,b,c) can be controlled and the corresponding geodesic is instantly displayed. (That's pretty cool.) First we obtain the background image:

i = Import[
   "http://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/\
    Equirectangular-projection.jpg/800px-Equirectangular-projection.jpg"]

Here is the code in its entirety:

Manipulate[
 {a, b, c} = {Cos[u] Cos[v], Sin[u] Cos[v], Sin[v]};
 Show[Graphics[{Texture[i], 
    Polygon[{{-\[Pi], -\[Pi]/2}, {\[Pi], -\[Pi]/2}, {\[Pi], \[Pi]/2}, {-\[Pi], \[Pi]/2}}, 
     VertexTextureCoordinates -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}}]}], 
  Plot[ArcTan[(a Cos[\[Lambda]] + b Sin[\[Lambda]]) / (-c)], {\[Lambda], -\[Pi], \[Pi]}, 
   PlotRange -> {Automatic, {-\[Pi]/2, \[Pi]/2}}, PlotStyle -> {Thick, Red}]],
   {u, 0, 2 \[Pi]}, {v, -\[Pi]/2, \[Pi]}]