Update rotation computation (#319)

* Update _rotations.py

Make sure that arccos is never given a value outside of its allowable range.
The previous version would occasionally yield a nan for stiff rods.

Co-authored-by: Seung Hyun Kim <skim0119@gmail.com>
Co-authored-by: Songyuan Cui <47127977+sy-cui@users.noreply.github.com>

Author: 3 people

elastica/_rotations.py

@@ -153,9 +153,13 @@ def _inv_rotate(director_collection: NDArray[np.float64]) -> NDArray[np.float64]
             )
         )
 
-        # TODO HARDCODED bugfix has to be changed. Remove 1e-14 tolerance
-        theta = arccos(0.5 * trace - 0.5 - 1e-10)
+        # Clip the trace to between -1 and 3.
+        # Any deviation beyond this is numerical error
+        trace = min(trace, 3.0)
+        trace = max(trace, -1.0)
+        theta = arccos(0.5 * trace - 0.5)
 
+        # TODO HARDCODED bugfix has to be changed. Remove 1e-14 tolerance
         vector_collection[0, k] *= -0.5 * theta / sin(theta + 1e-14)
         vector_collection[1, k] *= -0.5 * theta / sin(theta + 1e-14)
         vector_collection[2, k] *= -0.5 * theta / sin(theta + 1e-14)