From dbca47f113b64ad87ba1391a8a43bdd7bdf669e0 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 13 Jul 2020 12:38:21 +0200
Subject: [PATCH] better handling of corner cases

copy and paste from python code
---
 src/quaternions.f90 | 2 ++
 src/rotations.f90   | 6 +++---
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/quaternions.f90 b/src/quaternions.f90
index c337bd681..a396f7f67 100644
--- a/src/quaternions.f90
+++ b/src/quaternions.f90
@@ -464,6 +464,8 @@ subroutine selfTest
   real(pReal), dimension(4) :: qu
   type(quaternion)          :: q, q_2
 
+  if(dNeq(abs(P),1.0_pReal))                           call IO_error(0,ext_msg='P not in {-1,+1}')
+
   call random_number(qu)
   qu = (qu-0.5_pReal) * 2.0_pReal
   q  = quaternion(qu)
diff --git a/src/rotations.f90 b/src/rotations.f90
index f95c54e98..7879523d2 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -581,7 +581,7 @@ pure function om2eu(om) result(eu)
   real(pReal),             dimension(3)   :: eu
   real(pReal)                             :: zeta
 
-  if (abs(om(3,3)) < 1.0_pReal) then
+  if(dNeq(abs(om(3,3)),1.0_pReal,1.e-9_pReal)) then
     zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal)
     eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), &
           acos(om(3,3)), &
@@ -589,8 +589,8 @@ pure function om2eu(om) result(eu)
   else
     eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ]
   end if
-
-  where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
+  where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal
+  where(eu<0.0_pReal)          eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI])
 
 end function om2eu