basic projection question

[Gerald I. Evenden]

>Date: Wed, 7 Jul 93 02:37:21 -0500
>Message-Id: <9307070737.AA12583 at bushland.ecn.purdue.edu>
>From: Darrell McCauley <mccauley at ecn.purdue.edu>
>To: grassp-list at max.cecer.army.mil
>Subject: basic projection question
>
>Trying to make my programs more lat-long friendly...
>
>  G_begin_distance_calculations();
>
>I need to calculate the angle of a line made between two
>points (x1,y1) and (x2,y2).

First, you need to define what you mean by a "line" between two point
when dealing with geographic coordinate: a rhumb line (constant
azimuth) or geodesic (shortest distance).

>The distance between the points (the hypotenuse in
>geometric/planimetric terms) is
>
>  d = G_distance (x1, y1, x2, y2);
>
>In a planimetric projection, the angle (in degrees) would be:

What do you mean by "planimetric projection?"  Projections have nothing
to do with rhumb or geodesic distances---if that's what we are dealing
with.

>  angle = 180.0 / 3.14159 * acos ( (x1-x2) /d )
>
>but with lat-long I don't believe that this is the case
>(because the distance between two longitudes is different,
>depending upon what latitude you are at).  Would the 
>correct way to compute this angle be:
>
>  angle = 180.0 / 3.14159 
>        * acos ((G_distance(x1,y1,x2,y1)+G_distance(x1,y2,x2,y2)) /2.0 
>        / d);
>
>(that is, averaging the two distances at lat1 and lat2 to come
>up with the numerator for the inverse cosine quotient). 
>This approach seems sort of naive but I cannot think of
>anything better at the moment.

You have totally lost me.

>Is this correct approach? Is there a library function that I am
>missing which would do what I want?
>
>--Darrell

I am not sure that I understand your problem.  But you seem to be needing
the formulary for either a rhumbline or geodesic and I suspect the latter.
The program geod distributed with the PROJ.4 system will give you
forward and inverse geodesic computations with forward and back azimuths
at each point and great circle distances for both spherical and elliptical
Earths.  The elliptical formulary is too complex to post but the spherical
formulae are basically:

distance = 2*R*asin(sqrt(sin^2((lat2-lat1)/2) + cos(lat1)*cos(lat2)*
	sin^2((lon2-lon1)/2)))

Azi1-2 = atan2(cos(lat2)*sin(lon2-lon1), cos(lat1)*sin(lat2)-
	sin(lat1)*cos(lat2)*cos(lon2-lon1))

where lat1,lon1 and lat2,lon2 are geographic coordinates of points,
R is Earth's radius.  To get back azimuth from point 2 to point 1,
use second equation and with reversed coodinantes.  For more details
and verification of the above see Snyder's "Map Projections---A Working
Manual", USGS Prof.Paper 1395, p. 30.  For rhumb line computations, see
p. 46 of same publication.

Gerald (Jerry) I. Evenden   Internet: gie at charon.er.usgs.gov
voice: (508)563-6766          Postal: P.O. Box 1027
  fax: (508)457-2310                  N.Falmouth, MA 02556-1027