diff --git a/CHANGES.rst b/CHANGES.rst index 9c20ead..8ae1420 100644 --- 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 index 29334ea..9a57c76 100644 --- 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 index 554a732..2db95cf 100644 --- 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) {