NavList:

   

Reply

On 2019-06-04 8:16, Tony Oz wrote:
> I'm so interested in this Cosine calculator issue because my TI-83+ was 
constantly throwing an exception for several days in March while doing plain 
arc-cosine of no-fancy Law of Cosines H_C computation at LAN. Investigation 
has shown that the input value to the arccos() was equal to 1,0000000000001. 
I had to guard that function with an "*If X<=1* ..." check.

Although Atan2(y, x) is more robust than arc cosine, it too can give
unexpected results. For example, suppose the result is nominally 0 - 180
because y is never negative. (For clarity I'll use degrees instead of
radians.) But if y is a calculated value and very near zero, a trivial
computational error might make it negative, thereby making Atan2()
return a value near -180 instead of near +180.

Another issue is what happens in the special case where y = 0 and x < 0.
Does Atan2() return +180 or -180? The answer is not clear on this page,
which just says the value returned is in the interval [-pi, +pi], or in
degree measure, -180 <= Atan2() <= +180. My reference book on the C++
language is equally ambiguous.

http://www.cplusplus.com/reference/cmath/atan2/

On the other hand, this page (converted to degree measure) says -180 <
Atan2() <= +180. In other words, it never returns -180.

https://docs.microsoft.com/en-us/dotnet/api/system.math.atan2?view=netframework-4.8

Yet another issue is what happens if x = 0 and y = 0. Does the
implementation return zero, or does it generate a domain error?

ANGULAR SEPARATION WITH ATAN2

With vectors and the Atan2() function it's possible to compute the
angular separation between two points on the celestial sphere with high
accuracy whether the angle is a microsecond of arc or tens of degrees.

First, if necessary convert coordinates of the points to vectors (xyz
format). Call them vec1 and vec2. Then compute their dot product and
also the magnitude of their cross product. In pseudocode:

x = Dot(vec1, vec2)
y = Magnitude(Cross(vec1, vec2))

Then the angular separation is Atan2(y, x). This depends on the fact
that the dot product equals the cosine of the angle, times the product
of the magnitudes of the two vectors. And the magnitude of the cross
product equals the sine of the angle, times the product of the magnitudes.

In effect, you get the cosine and sine of the angle, multiplied by the
same constant. Thus, a rectangular to polar conversion finds the angle
accurately regardless of its size. The result is always in the first or
second quadrant because the magnitude of a vector is never negative.
This is the algorithm in my SofaJpl astronomy toolbox DLL for Windows,
and the Lunar program (which depends on SofaJpl).

Some defensive coding is necessary to avoid the issues I mentioned
above. See the International Astronomical Union SOFA library function Sepp:

http://www.iausofa.org/2018_0130_C/sofa/sepp.html

It traps the degenerate case where both vectors are null, but not the
ambiguous case where y = 0 and x < 0. The corresponding routine in
SofaJpl is merely a shell for the SOFA routine, and it too lacks any
defensive coding for the ambiguous case. I've made a note to fix that in
the next release. This is not the first time that a discussion on an
unrelated topic caused me to look at my own code and find room for
improvement.

   

Reply