fix problem with C version of angle function (#223)

fix problem with C version of angle function
spacetelescope · Oct 21, 2022 · 1aca4c3 · 1aca4c3
1 parent 82a5a8c
commit 1aca4c3
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 8 deletions.
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -11,6 +11,9 @@ Release Notes
 
 - Set Python min version to 3.8. [#219]
 
+- Fix problem reported where angle of 180 degrees results in an
+  error with the C version of the code. [#223]
+
 
 1.2.22 (04-January-2022)
 ========================

diff --git a/spherical_geometry/great_circle_arc.py b/spherical_geometry/great_circle_arc.py
@@ -327,10 +327,8 @@ def angle(A, B, C, degrees=True):
 
         A, B, C = np.broadcast_arrays(A, B, C)
 
-        ABX = _fast_cross(A, B)
-        ABX = _cross_and_normalize(B, ABX)
-        BCX = _fast_cross(C, B)
-        BCX = _cross_and_normalize(B, BCX)
+        ABX = _cross_and_normalize(A, B)
+        BCX = _cross_and_normalize(C, B)
         X = _cross_and_normalize(ABX, BCX)
         diff = inner1d(B, X)
         inner = inner1d(ABX, BCX)

diff --git a/src/math_util.c b/src/math_util.c
@@ -691,13 +691,11 @@ DOUBLE_angle(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
         load_point_qd(ip2, is2, B);
         load_point_qd(ip3, is3, C);
 
-        cross_qd(A, B, TMP);
-        cross_qd(B, TMP, ABX);
+        cross_qd(A, B, ABX);
 
         if (normalize_qd(ABX, ABX)) return;
 
-        cross_qd(C, B, TMP);
-        cross_qd(B, TMP, BCX);
+        cross_qd(C, B, BCX);
 
         if (normalize_qd(BCX, BCX)) return;
 
@@ -708,6 +706,17 @@ DOUBLE_angle(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
         dot_qd(B, X, &diff);
         dot_qd(ABX, BCX, &inner);
 
+        /* The following threshold is currently arbitrary and
+        is based on observed errors in that value and several 
+        orders of magnitude larger than those. One day someone
+        should determine how qd is producing those errors, 
+        but for now...
+        */
+        if (fabs(inner.x[0]) == 1.0 && fabs(inner.x[1]) < 1e-60) {
+            inner.x[1] = 0.;
+            inner.x[2] = 0.;
+            inner.x[3] = 0.;
+        }
         if (inner.x[0] != inner.x[0] ||
             inner.x[0] < -1.0 ||
             inner.x[0] > 1.0) {