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]}]
How about a method which is both more accurate and faster? This is
provided by GeographicLib. Comparative timings (C++
implementations on a 2.66GHz Intel processor, using g++) are:
Vincenty direct: 1.11 us
GeographicLib::Geodesic::Direct: 0.88 us
GeographicLib::GeodesicLine::Position: 0.37 us
GeographicLib::GeodesicLine::ArcPosition: 0.31 us
The accuracy of Vincenty's formulas is about 0.1 mm, while the accuracy
of the GeographicLib algorithms about 0.01 um. Geodesic::Direct does a
straight solution of the direct problem. It's somewhat faster than
Vincenty because it's non-iterative and because it uses Clenshaw
summation to evaluate the trigonometric series. GeodesicLine::Position
allows you to calculate many points along a single geodesic about
2.4 times faster. If you merely want some
points on a geodesic which are approximately equally spaced (e.g., for plotting it), you can use
GeodesicLine::ArcPosition and shave a little extra time off the
computation. You can reduce the time still further by reducing the order of the
series used by GeographicLib from 6 to 3 by compiling with
-DGEOGRAPHICLIB_GEODESIC_ORDER=3
The accuracy is then 0.04 mm, i.e.,
comparable to, but slightly better than, Vincenty.
A cookbook recipe for solving the equivalent problem on a sphere is
given by the Wikipedia entry on great-circle navigation.
.
Best Answer
Questions like this can often be answered by solving triangles using laws of spherical trigonometry. You can look these up: they are referenced below. Normally, when you know three of the six parts of a spherical triangle, the laws let you find the other three in terms of the sines and cosines of the parts you know. The trick is to draw a useful triangle. When you're dealing with meridians, think of using a pole (which is common to all meridians) as one of the triangle's vertices. This suggests focusing on the triangles N-M-T1 or, equivalently, N-M-T2. Here are the details, a general solution, and a worked example.
Let the point in question be at latitude phi and the circle's radius be r (expressed as an angle, not as a distance. The angle in radians is the radial distance divided by the sphere's radius). We start by finding the difference between the local longitude and the extreme longitudes; let this difference be lambda.
Because the extreme meridians are tangent to the circle, angles M-T1-N and M-T2-N are right angles. This is nice because the sine and cosine of a right angle are easy to find and simple (they equal 1 and 0, respectively). The spherical Law of Sines says
Solving this gives
For example, let r = 45 degrees and phi = -30 degrees. The solution is lambda = 54.7356 degrees.
This much you already found out: it's the formula quoted in the question. Now we can apply the spherical Law of Cosines to relate the latitude phi' of the extremal points T1 and T2 to the radius r: this is what we are looking for.
The solution for the common latitude of the tangent points is
(for a circle centered at latitude phi of radius r radians.)
The remainder of the reply discusses and illustrates this result.
Notice in particular that when we start at the equator, phi = 0, easily giving phi' = 0 as expected. Otherwise, the sine of phi' is larger (in size) than the sine of phi, implying the solution is closer to the nearest pole, again as expected. In the example, phi' works out to -45 degrees.
The south pole is near the bottom of this pseudo-3D image of a sphere. The meridians start at the circle's center and go in 15 degree increments to either side, showing +-15, +-30, and +-45 degree longitudinal displacements, then stop at the limiting values of +-54.7356 degrees. The red dots are situated along these limiting meridians at latitudes of -45 degrees.
These formulas work for any circle that does not include either pole.
In this image, the circle's center (black dot) is just north of 60 degrees north latitude. Approximate circles of latitude are shown as dashed curves in 10 degree increments. The radius, equal almost to 30 degrees of arc, takes this circle almost to the north pole (where the meridians are converging at the top). The extremal meridians are therefore spread far apart (lambda is about 70 degrees) and, accordingly, the points of tangency are also close to the north pole, near 80 degrees latitude. This shows why the points of tangency (red dots) are usually closer to the nearest pole than the circle's center is.