silver@csu-cs.UUCP (07/04/83)
I'm working on a program (gcdist) which produces tables of distances and/or headings between places, given latitudes and longitudes. (I'd be glad to post a copy when I'm done, by the way.) It's finished EXCEPT for one stone wall I ran into. The VNR Concise Encyclopedia of Mathematics, 1975, section 12.2 explains how to compute headings, but they don't give a clear algorithm for disambiguating the results of asin(). You get two supplementary headings (angles beta and pi - beta) from asin (sin (b) * sin (alpha) / sin (a)), where beta and alpha are internal angles (alpha is the difference in longitude) and b and a are opposite sides (b is north pole to other place; a is the great circle distance between places). Since "the greater angle is opposite the greater side", you are supposed to take the value of beta that is "greater or smaller than the angle alpha, according as the side b is greater or smaller than the side a." However, what do you do when BOTH values are either greater or smaller? For example, suppose the places are at 0 N, 0 E and 1 N, 1 E. Here, alpha is very small (1 degree) while the values for beta are 44.996 deg and 135.004 deg. The former is obviously correct, but how can the program know that? I give up! After a week of thinking about it, I admit defeat and call for help. Any math wizards know the answer? Please mail me, and thanks in advance! Alan Silverstein, Hewlett-Packard Fort Collins Systems Division, Colorado ucbvax!hplabs!hpfcld!ajs, 303-226-3800 x3053, N 40 31'31" W 105 00'43"
bill@utastro.UUCP (07/09/83)
Here is the best way to solve a spherical triangle. Let the sides
and angles be labelled as below:
/\
/ C\
/ \
b / \ a
/ \
/ \
/ A B \
----------------
c
Capital letters are the angles and lower case letters the opposite
sides. In your case, angle C is at the north pole, so that
C = longitude difference between the two points. Sides a and b
are the colatitudes of the two points of interest (colatitude =
pi/2 - latitude). You want angle A (the bearing from
A to B) and side c (the distance from A to B).
The following three quantities can be calculated from the known
quantities a, b and C:
x = sin c cos A = cos a sin b - sin a cos b cos C
y = sin c sin A = sin a sin C
z = cos c = cos a cos b + cos b cos a cos C
You should calculate all three quantities and then do the following:
The ratio of the first two quantities
is tan A; but don't use the ordinary arctangent function! Instead,
use the relation angle A = atan2(y,x), which gives you the angle
*resolved for quadrant*. By definition, sin c >= 0, because all
sides of the spherical triangle are less than pi, so we don't
have to know sin c at this point.
Now compute w = sin c = sqrt (x*X + y*y); as remarked above, it is
non-negative. Finally, compute side c = atan2 (w, z). For the reasons
given below, you should *not* calculate c from the acos of z!
This method has several advantages:
1) It automatically resolves the quadrant issue;
2) It is accurate even near A = pi/2 or c = 0;
the other approach (using asin and acos)
loses accuracy because the sin and cos have
extrema near pi/2 and 0, respectively.
3) It works for southern latitudes as well. Keep
point C at the north pole; the math will still
work out.
The same method is the method of choice on a pocket calculator, as
well. Most of the scientific calculators have a vector-to-polar
key that is the analog of the function atan2 (y, x).
Bill Jefferys
Astronomy Dept
University of Texas
Austin TX 78712
(...ucbvax!nbires!ut-ngp!utastro!bill)
(...decvax!eagle!ut-ngp!utastro!bill)
( utastro!bill@utexas-11)