From: PT-PMR IT Dept Newsgroups: sci.math Subject: Re: Finding the Distance Date: 15 Feb 1995 12:32:16 GMT kegolf@iquest.net (Kevin Egolf) wrote: > > I am trying to get a little help on figuring out distances. Is there > some sort of function to get 'as the crow flies' distance between > two places if you have the latitude and longitude of both places? > Not so much a function as a formula - the 'cosine rule' of spherical trigonometry. If a spherical triangle (a shape bounded by 3 arcs of great circles subtending angles of a, b and c at the centre of the sphere) has angles A, B and C opposite the correspondingly lettered sides, then the formula is cos(a)=cos(b)cos(c)+sin(b)sin(c)cos(A) In your case you will take A as the angle at the North Pole (the difference in longitudes), b and c as the angles of the points from the pole (90deg - latitude) and a will be the distance you require. [sig deleted - djr] ============================================================================== Newsgroups: sci.math.num-analysis From: hbaker@netcom.com (Henry G. Baker) Subject: Draft FAQ: Great Circle Distances from latitudes & Longitudes Date: Fri, 12 May 1995 16:05:49 GMT Soliciting comments/reviews on the following: Draft FAQ: Computing "Great Circle Distances" from latitudes and longitudes. This problem (and its solutions) are more than 500 years old. Henry Baker (hbaker@netcom.com), May 1995. Problem: Compute the great circle distance (either angular or linear) between the two points whose coordinates are given by (latitude,Longitude) pairs. Specifically, given P1=(l1,L1), P2=(l2,L2), compute the angular distance d between these two points on the globe. The linear distance is then r*d, where r is the radius of the globe (assuming that d is given in radians). Plan: Transform coordinates into Cartesian coordinates, compute the "dot" product, and then take the arccosine of the result. We can either perform laborious trigonometric calculations, or we can utilize geometric insight to realize that the great circle distance doesn't depend upon the actual values of the longitudes L1,L2, but only on their _difference_ dL=L2-L1. We therefore translate the problem to P1'=(l1,0), P2'=(l2,L2-L1)=(l2,dL). P1': x1=cos(l1), y1=0, z1=sin(l1). P2': x2=cos(l2)cos(dL), y2=cos(l2)sin(dL), z3=sin(l2). P1'.P2' = cos(d) = cos(l1)cos(l2)cos(dL) + sin(l1)sin(l2). We are now mathematically (but not computationally) done, because we have the following expression for the angular distance we seek: d = arccos(cos(l1)cos(l2)cos(dL) + sin(l1)sin(l2)). Discussion: There is a problem with this formula, however. If dL is very small, cos(dL) is essentially 1, and arccos(1)=0. In other words, we lose all accuracy in the (common) case where dL and |l1-l2| are both small. The problem is that cosines and arccosines do not have the property that f(0)=0. What we really want is a formula based on sines/arcsines or tangents/arctangents, which _do_ have this property. We therefore use the following centuries-old half-angle trick: sin(dL/2)^2 = (1-cos(dL))/2, or equivalently, cos(dL) = 1-2sin(dL/2)^2. cos(d) = cos(l1)cos(l2)(1 - 2 sin(dL/2)^2) + sin(l1)sin(l2) = cos(l1)cos(l2) - 2 cos(l1)cos(l2)sin(dL/2)^2 + sin(l1)sin(l2) (recognizing the cosine difference formula) = cos(l2-l1) - 2 cos(l1)cos(l2)sin(dL/2)^2 = cos(dl) - 2 cos(l1)cos(l2)sin(dL/2)^2 [dl = l1-l2] (using the half-angle trick again with dl) = 1 - 2 sin(dl/2)^2 - 2 cos(l1)cos(l2)sin(dL/2)^2 (using the half-angle trick a 3rd time with d) 1 - 2 sin(d/2)^2 = 1 - 2 sin(dl/2)^2 - 2 cos(l1)cos(l2)sin(dL/2)^2 sin(d/2)^2 = sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2 (This now looks a lot like a Pythagorean formula, except with half-sines and a correction coefficient of cos(l1)cos(l2).) sin(d/2) = sqrt(sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2) d = 2 asin( sqrt( sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2 ) ) This formula is now accurate when dl and/or dL are near zero. [As a check, if l1=l2=l, this formula simplifies to: d = 2 asin(min(1,cos(l)sin(dL/2))). If l=0, then d = dL. Furthermore, if dL=0, this formula simplifies to: d = 2 asin(sin(dl/2)) = 2 (dl/2) = dl.] Furthermore, sin()^2 >=0, and cos() >=0, since -pi/2 <= l1,l2 <= pi/2, so the argument to the square root is always >=0. It is also mathematically <=1, but errors may violate this by small amounts, so use min to guarantee that the argument to asin is in the range [0,1]. We finally have: d = 2 asin(min(1,sqrt(sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2))) Although in current subroutine libraries this appears to be a more expensive calculation, classical books tabulate the function versine(x) = vers(x) = 2 haversine(x) = 2 sin(x/2)^2 and its inverse. Therefore, if you're going to do this calculation a lot, you should also have these specialized subroutines for versine and its inverse. Since vers(x) = 2 sin(x/2)^2 = 1 - cos(x), its Taylor series is x^2/2! - x^4/4! + x^6/6! - x^8/8! + ... = cos(0) - cos(x), so versine is no more difficult to compute than cosine. Versine is more useful, however, because it is much more accurate near zero. Since accuracy is hard to get and easy to lose, the primitive trignometric functions found in your subroutine library should be _versine_ and _arcversine_ instead of _cosine_ and _arccosine_, and those (usually inappropriately) desiring cosine and arccosine can easily construct them as the macros: cos(x) = (1 - vers(x)) acos(x) = (pi/2 - asin(x)). [Also realize that the Earth is not really a sphere, but an oblate spheroid, so even these calculations are not accurate for the Earth.] [Historial note: The ancient mathematicians (astronomers and astrologers, mostly, who actually had to _compute numerically_ for a living instead of doing just symbolic algebra) used the function chord(x) = cd(x) = 2 sin(x/2) instead of sin(x). Using cd(x), the Pythagoreanness of the formula above is particularly striking: cd(d)^2 = cd(dl)^2 + cos(l1) cos(l2) cd(dL)^2 = cd(dl)^2 + (1-vers(l1)) (1-vers(l2)) cd(dL)^2 The Arab mathematicians circa 1000 AD used the "versed sine", or versine, in addition to the sine and chord. The cosine was almost never used, and apparently became popular much later as a result of Euler's identity: exp(iy) = cos(y) + i sin(y). In other words, cosines are important for _algebraic manipulation_, but not for _numerical calculation_. Apparently, those doing numerical calculations on computers still have much to learn from the ancient human computers. Reference: Smith, David E. History of Mathematics, Vol. II. Dover Publs, NY, 1925.] Copyright (c) 1995 by Henry Baker. All rights reserved. ============================================================================== From: sofiya@helium.gas.uug.arizona.edu (Sofiya B Vasina) Newsgroups: sci.math Subject: Re: great circle distance Date: 26 Jun 1995 05:20:21 GMT In article <3sj4h2$jjh$1@perth.DIALix.oz.au>, Neil Gerace wrote: >I'm after any simple formula to calculate the great circle distance >between any two points on the Earth's surface, given their lats and longs. gcd = R * arccos (cos(lat1)*cos(long1)*cos(lat2)*cos(long2) + cos(lat1)*sin(long1)*cos(lat2)*sin(long2) + sin(lat1)*sin(lat2)), where R is the radius of the earth. [sig deleted - djr] ============================================================================== From: tadam@megamed.com (Brian Adam) Newsgroups: sci.math,sci.geo.meteorology Subject: Re: Distance between two points on Earth Date: Mon, 20 Nov 1995 22:17:28 GMT dennis@emmi.physik.tu-berlin.de (D. S.) wrote: >Hello, >I am looking for an algorithm to calculate the distance between two >points on a sphere (in my case on Earth). I know that it has something >to do with spherical trigonometry but couldn't find a solution in my >books. >I would be happy with some kind of formula or algorithms in C, FORTRAN >or any other programming language. The input variables are longitude >and latitude of the two points. >As I need this for a meteorological program I post this also to >sci.geo.meteorology. >With kind regards, > Dennis >--- >Dennis Schulze Free University of Berlin >Fax/Voice: (+49 30) 793 49 51 Department of Meteorology > Email: dennis@bibo.met.fu-berlin.de >http://www.met.fu-berlin.de/~dennis/ #include Dug this up in from my old synoptic met class. It's not in any code, but perhaps you can write a code for it. Let a=lat of station 1 b=lat of station2 c=long of station1 d=long of station 2 Now x(radians)=(cos)-1 {cos(a)cos(b)cos(c-d) + sin(a)sin(b)} Then Distance=Radius(earth)*x(radians) Notes: 1. (cos)-1=inverse cosine (the only way I could type it) 2. for west longitudes use -(c) or -(d) in the computation 3. Note that x is in radians. Hope this helps. Brian Adam ============================================================================== From: Anthony Hugh Back Newsgroups: sci.math Subject: Re: how to determine distance between lat/long points? Date: Thu, 25 Jan 1996 16:01:58 GMT In article: ** Craig Cook ** writes: > > If I have the latitude and longitude of two points > > How do I determine the number of miles between them? > This can be done very easily by using the scalar product of two vectors to find the angle between those vectors. If the vectors are OA and OB where A and B are the two points on the surface of the earth and O is the centre of the earth. The scalar product gives OA*OB*cos(AOB) = R^2*cos(AOB) where R = radius of the earth. Having found angle AOB the distance between the points is R*(AOB) with AOB in radians. To find the scalar product we need the coordinates of the two points. Set up a three dimensional coordinate system with the x-axis in the longitudinal plane of OA and the xy plane containing the equator, the z-axis along the earth's axis. With this system, the coordinates of A will be Rcos(latA), 0, Rsin(latA) and the coords of B will be Rcos(latB)cos(lonB-lonA),Rcos(latB)sin(lonB-lonA),Rsin(latB) The scalar product is given by xA*xB + yA*yB + zA*zB = R^2cos(latA)cos(latB)cos(lonB-lonA)+ R^2sin(latA)sin(latB) Dividing out R^2 will give cos(AOB) cos(AOB)=cos(latA)cos(latB)cos(lonB-lonA)+sin(latA)sin(latB) This gives AOB, and the great circle distance between A and B will be R*(AOB) with AOB in radians. -- Anthony Hugh Back ============================================================================== From: bill@clyde.as.utexas.edu (William H. Jefferys) Newsgroups: sci.math Subject: Re: how to determine distance between lat/long points? Date: 3 Feb 1996 14:02:58 GMT In article , Craig Cook wrote: >If I have the latitude and longitude of two points >How do I determine the number of miles between them? Here's a method not using spherical trigonometry that is perhaps illuminating. Compute the vector pointing from the center of the (assumed spherical) Earth to each point. For point_1 this is x_1 = cos(lat_1)*cos(long_1) y_1 = cos(lat_1)*sin(long_1) z_1 = sin(lat_1) and similarly for point_2. Now, the cosine of the angle between point_1 and point_2 is the dot product of the two; and the sine of the angle between them is the modulus of the cross product between the two. In particular, the dot product gives you the spherical trig formula that the other two responses used, and one can get the result from the arc cosine of this. A point of practicality is that the arc cosine of the dot product loses accuracy if the two points are very close together (because the expansion of the cosine is 1 - theta^2/2 + ..., so computing the arc cosine of a small angle involves taking the square root of a very small number). In extreme cases, you can lose almost all significance doing this. So for algorithmic purposes, where one wants an accurate result, guaranteed, no matter what the inputs, it is better to compute _both_ the sine _and_ the cosine (as above) and then use the two-argument arctangent routine to calculate the angle (but use double precision when calculating the modulus). This will give maximum significant digits in all cases. The result is in angular measure, of course. To get linear distance you have to convert this to nautical miles. Things get more complicated if the Earth is not assumed spherical. Bill -- Bill Jefferys/Department of Astronomy/University of Texas/Austin, TX 78712 E-mail: bill@clyde.as.utexas.edu | URL: http://quasar.as.utexas.edu Finger for PGP Key: F7 11 FB 82 C6 21 D8 95 2E BD F7 6E 99 89 E1 82 -------------------------------------------------------------------------- ============================================================================== Newsgroups: sci.math From: hbaker@netcom.com (Henry Baker) Subject: Re: Long/Lat Distance Date: Mon, 24 Nov 1997 21:25:35 GMT In article <65c9mv$4r7@bgtnsc03.worldnet.att.net>, Steve Andrews <"steveandrews@worldnet.att.net"@worldnet.att.net> wrote: > I am looking for a formula to calculate the distance between two points > on the earth's surface given the longitudes and latitudes of those two > points. ftp://ftp.netcom.com/pub/hb/hbaker/FAQ-lat-long.txt