From 1aed224c3b0a249fe33babe6eb0839da8d0bbe9a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 6 Apr 2019 11:23:31 +0200 Subject: [PATCH] numerically more stable avoids division by zero --- src/rotations.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index b899adacb..7c29b03bc 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -777,12 +777,12 @@ pure function qu2ax(qu) result(ax) real(pReal) :: omega, s - omega = 2.0 * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)) ! if the angle equals zero, then we return the rotation axis as [001] - if (dEq0(omega)) then - ax = [ 0.0, 0.0, 1.0, 0.0 ] + if (dEq0(sqrt(qu%x**2+qu%y**2+qu%z**2))) then + ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] elseif (dNeq0(qu%w)) then s = sign(1.0_pReal,qu%w)/sqrt(qu%x**2+qu%y**2+qu%z**2) + omega = 2.0_pReal * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)) ax = [ qu%x*s, qu%y*s, qu%z*s, omega ] else ax = [ qu%x, qu%y, qu%z, PI ]