NavList:
A Community Devoted to the Preservation and Practice of Celestial Navigation and Other Methods of Traditional Wayfinding
Re: distance with atan2
From: George Huxtable
Date: 2005 Nov 12, 12:03 -0000
From: George Huxtable
Date: 2005 Nov 12, 12:03 -0000
Paul Hirose wrote- =========================================== One way to calculate distance (separation angle) on a sphere is to use vectors and the atan2 function or rectangular to polar conversion. Set x to the dot product of the two vectors, and y to the magnitude of their cross product. The separation angle is atan2(y, x), or you can get the angle by converting x and y to polar form. For example, Meeus calculates the separation angle between two stars: 14h 15m 39.7s +19°10′57″ 13h 25m 11.6s -11°09′41″ He uses spherical trig. Using the vector method, start by forming the rectanglar coordinates on the celestial sphere: -.783786 -.526988 .328578 -.914079 -.356354 -.193573 dot product = .840633 (x) cross product magnitude = .541606 (y) Converting (x, y) to polar form gives angle 32.7930°, which agrees with Meeus. An advantage of the vector method is its accuracy at all angles. A large portion of Meeus' chapter on separation angle discusses ways to avoid the inaccuracy of the traditional spherical trig formula near 0 or 180 degrees. ========================================= Comments from George. Paul often comes up with an interesting slant on things. I'm not fluent in vector manipulation, but perhaps I can see what he is getting at. I'll try to spell out some of the missing detail, and Paul can correct me if I've got it wrong. Three axes at right angles are set up, from the Earth's centre. The x axis is along the plane of the equator, to where it meets the Greenwich meridian., at lat = 0, long = 0. The y axis is also along the plane of the equator, toward a point 90 degrees East of Greenwich, at lat = 0, long = 90 degrees East. The z axis is toward the North pole. The radius of the Earth (assumed spherical), is taken to be 1 unit. Then a point A on the Earth's surface, at latA, longA, will have these three coordinates in space- Ax = cos latA cos longA Ay = cos latA sin longA Az = sin latA Paul takes an example from Meeus that uses star sky-angles of Right-Ascension (RA), and dec. but to my mind coordinates on the Earth's surface of lat and long are easier to visualise, and the mathematics is identical when the Earth is spherical. Note that here we are working in what's called a right-handed coordinate system, beloved of mathematicians, which assumes East is positive, but which may not be most appropriate for us navigators, and some modifications may be called for later. If you take Paul's example, but put long instead of RA, and lat instead of Dec, you can chech that his position A, at RA of 14h 15m 39.7 s, would translate, at 15d per hour, into longA = 213.9154 degrees of long (Easterly). Dec of +19d 10' 57" would become latA = +19.1825 degrees. You can work out for yourself the corresponding values for position B. Now, from that lat and long, you can deduce the three coordinates of A, (Ax, Ay, And Az), to be -.783786, -.526988, and +.327578, just as Paul gives them above. These are the components of a "vector" A, but you shouldn't let that word intimidate you. You can check if you've got the answer right, by squaring those three numbers, adding them, and taking the square-root. The result should be the length of the vector A. Because A stretches from the Earth's centre to a point on its surface, and we defined the Earth's surface to be 1 unit, then the length of A should be exactly 1. You can check it. And similarly, those for B, (Bx, By, Bz) become -.914079, -.356354, -.193573, as quoted by Paul. Then Paul works out the "dot-product" of the vectors A and B. In my learning days we called it the scalar product. I haven't remembered how to do these tricks, but have turned to Joos, "Theoretical Physics", a familiar student textbook from 50 years ago. The dot-product is given by AxBx + AyBy +AzBz. Paul has calculated this to be .840633, and I agree. Next he works out the length of the "cross-product", or vector product, which is a bit more complicated. First, you have to work out 3 terms as follows- (AyBz - AzBy), (AzBx - Ax Bz), and (AxBy - AyBx) Then square each term, add those three squares, and take the square-root of the result. If you've done all that right, I expect you will end up with the same result as Paul, of .541606. Finally, to get the distance between A and B, you have to work out the arctan of the ratio between the dot-product and the cross-product. This is best done using an ATAN2 or POL function which will put the answer into the right quadrant, as Paul has explained. In degrees, his result is 32.7930. We can multiply by 60 to turn it into nautical miles, and get 1967.6 miles as the great-circle distance between A and B. Note that in using the ATAN2 function, he has labelled the two input quantities x and y. That's simply because it's how the computer would label those inputs, those labels having no connection with the initial x and y axes. ================================ You may think all this is a long-winded way of going about obtaining the great-circle distance from A to B, and indeed it is. However, the procedures involved are all part of a mathematician's box of standard tools, and a computer can do those tasks in a twinkling. Also, as Paul explains, there are some advantages in precision, when the distance is very small, for which the standard formula is not well-suited.. It's certainly a different way of handling spherical problems to what we are used to. It leaves one or two questions, which Paul may be able to help with. If a POL function is used to provide an angle output, as the A of a (R,A) pair, is the resulting R of any value, or must it be discarded? That's a similar question to what Bill has been asking. Can that vector method also give the initial great-circle course, leaving A towards B? If so, do any mods need to be made to conform with common navigator's notations (LHA positive westwards, azimuth clockwise from North)? George. contact George Huxtable at george@huxtable.u-net.com or at +44 1865 820222 (from UK, 01865 820222) or at 1 Sandy Lane, Southmoor, Abingdon, Oxon OX13 5HX, UK.