Skip to content

Commit

Permalink
Update _rotations.py
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
nmnaughton committed Jan 5, 2024
1 parent cdf6ec2 commit 300e381
Showing 1 changed file with 5 additions and 2 deletions.
7 changes: 5 additions & 2 deletions elastica/_rotations.py
Original file line number Diff line number Diff line change
Expand Up @@ -146,9 +146,12 @@ def _inv_rotate(director_collection):
)
)

# 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 -3 and 3. It any deviation beyond this is numerical error
trace = min(trace,3.0)
trace = max(trace,-3.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)
Expand Down

0 comments on commit 300e381

Please sign in to comment.