diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 167552b68..ee497e72b 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -669,7 +669,7 @@ class Rotation:
         """Rotation matrix to Bunge-Euler angles."""
         with np.errstate(invalid='ignore',divide='ignore'):
             zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2)
-            eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4),
+            eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-9),
                           np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]),
                                     np.pi*0.5*(1-om[...,2,2:3]),
                                     np.zeros(om.shape[:-2]+(1,)),
@@ -679,30 +679,30 @@ class Rotation:
                                     np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta)
                                    ])
                           )
-        eu[np.abs(eu)<1.e-6] = 0.0
+        eu[np.abs(eu)<1.e-8] = 0.0
         eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
         return eu
 
     @staticmethod
     def om2ax(om):
         """Rotation matrix to axis angle pair."""
+        return Rotation.qu2ax(Rotation.om2qu(om)) # HOTFIX
         diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2],
                                    om[...,2,0:1]-om[...,0,2:3],
                                    om[...,0,1:2]-om[...,1,0:1]
                                  ])
-        diag_delta[np.abs(diag_delta)<1.e-6] = 1.0
         t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,))
         w,vr = np.linalg.eig(om)
         # mask duplicated real eigenvalues
         w[np.isclose(w[...,0],1.0+0.0j),1:] = 0.
         w[np.isclose(w[...,1],1.0+0.0j),2:] = 0.
         vr = np.swapaxes(vr,-1,-2)
-        ax = np.where(np.abs(diag_delta)<0,
+        ax = np.where(np.abs(diag_delta)<1e-12,
                              np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
                       np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \
                       *np.sign(diag_delta))
         ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))])
-        ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0]
+        ax[np.abs(ax[...,3])<1.e-8] = [ 0.0, 0.0, 1.0, 0.0]
         return ax
 
     @staticmethod
@@ -804,6 +804,7 @@ class Rotation:
     @staticmethod
     def ax2om(ax):
         """Axis angle pair to rotation matrix."""
+        return Rotation.qu2om(Rotation.ax2qu(ax))  # HOTFIX
         c = np.cos(ax[...,3:4])
         s = np.sin(ax[...,3:4])
         omc = 1. -c
diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py
index 5219f3260..ab503f402 100644
--- a/python/tests/rotation_conversion.py
+++ b/python/tests/rotation_conversion.py
@@ -134,25 +134,26 @@ def om2qu(a):
 
 def om2eu(om):
     """Rotation matrix to Bunge-Euler angles."""
-    if not np.isclose(np.abs(om[2,2]),1.0,1.e-4):
+    if not np.isclose(np.abs(om[2,2]),1.0,1.e-9):
         zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
         eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
                        np.arccos(om[2,2]),
                        np.arctan2(om[0,2]*zeta, om[1,2]*zeta)])
     else:
         eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0])                    # following the paper, not the reference implementation
-    eu[np.abs(eu)<1.e-6] = 0.0
+    eu[np.abs(eu)<1.e-8] = 0.0
     eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
     return eu
 
 def om2ax(om):
     """Rotation matrix to axis angle pair."""
+    return qu2ax(om2qu(om)) # HOTFIX
     ax=np.empty(4)
 
     # first get the rotation angle
     t = 0.5*(om.trace() -1.0)
     ax[3] = np.arccos(np.clip(t,-1.0,1.0))
-    if np.abs(ax[3])<1.e-6:
+    if np.abs(ax[3])<1.e-8:
         ax = np.array([ 0.0, 0.0, 1.0, 0.0])
     else:
         w,vr = np.linalg.eig(om)
@@ -160,8 +161,7 @@ def om2ax(om):
         i = np.where(np.isclose(w,1.0+0.0j))[0][0]
         ax[0:3] = np.real(vr[0:3,i])
         diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
-        diagDelta[np.abs(diagDelta)<1.e-6] = 1.0
-        ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta))
+        ax[0:3] = np.where(np.abs(diagDelta)<1e-12, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta))
     return ax
 
 #---------- Bunge-Euler angles ----------
@@ -229,6 +229,7 @@ def ax2qu(ax):
 
 def ax2om(ax):
     """Axis angle pair to rotation matrix."""
+    return qu2om(ax2qu(ax)) # HOTFIX
     c = np.cos(ax[3])
     s = np.sin(ax[3])
     omc = 1.0-c
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 4dc98a797..adf507e49 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -107,10 +107,10 @@ class TestRotation:
             print(m,o,rot.as_quaternion())
             assert ok and np.isclose(np.linalg.norm(o),1.0)
 
-    @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om)])
-                                                 #(Rotation.om2eu,Rotation.eu2om),
-                                                 #(Rotation.om2ax,Rotation.ax2om),
-                                                 #(Rotation.om2ro,Rotation.ro2om),
+    @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om),
+                                                 (Rotation.om2eu,Rotation.eu2om),
+                                                 (Rotation.om2ax,Rotation.ax2om),
+                                                 (Rotation.om2ro,Rotation.ro2om)])
                                                  #(Rotation.om2ho,Rotation.ho2om),
                                                  #(Rotation.om2cu,Rotation.cu2om)])
     def test_matrix_internal(self,default,forward,backward):
@@ -121,12 +121,12 @@ class TestRotation:
             print(m,o,rot.as_quaternion())
             assert ok and np.isclose(np.linalg.det(o),1.0)
 
-    @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu)])
-                                                #(Rotation.eu2om,Rotation.om2eu),
-                                                #(Rotation.eu2ax,Rotation.ax2eu),
-                                                #(Rotation.eu2ro,Rotation.ro2eu),
-                                                #(Rotation.eu2ho,Rotation.ho2eu),
-                                                #(Rotation.eu2cu,Rotation.cu2eu)])
+    @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu),
+                                                 (Rotation.eu2om,Rotation.om2eu),
+                                                 (Rotation.eu2ax,Rotation.ax2eu),
+                                                 (Rotation.eu2ro,Rotation.ro2eu),
+                                                 (Rotation.eu2ho,Rotation.ho2eu),
+                                                 (Rotation.eu2cu,Rotation.cu2eu)])
     def test_Eulers_internal(self,default,forward,backward):
         for rot in default:
             m = rot.as_Eulers()
@@ -140,6 +140,22 @@ class TestRotation:
             print(m,o,rot.as_quaternion())
             assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all()
 
+    @pytest.mark.parametrize('forward,backward',[(Rotation.ax2qu,Rotation.qu2ax),
+                                                 (Rotation.ax2om,Rotation.om2ax),
+                                                 (Rotation.ax2eu,Rotation.eu2ax),
+                                                 (Rotation.ax2ro,Rotation.ro2ax),
+                                                 (Rotation.ax2ho,Rotation.ho2ax),
+                                                 (Rotation.ax2cu,Rotation.cu2ax)])
+    def test_axis_angle_internal(self,default,forward,backward):
+        for rot in default:
+            m = rot.as_axis_angle()
+            o = backward(forward(m))
+            ok = np.allclose(m,o,atol=atol)
+            if np.isclose(m[3],np.pi,atol=atol):
+                ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
+            print(m,o,rot.as_quaternion())
+            assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7
+
 
     @pytest.mark.parametrize('degrees',[True,False])
     def test_Eulers(self,default,degrees):
@@ -177,7 +193,7 @@ class TestRotation:
             if np.isclose(m[3],np.pi,atol=atol):
                 ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
             print(m,o,rot.as_quaternion())
-            assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9
+            assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7
 
     @pytest.mark.parametrize('P',[1,-1])
     @pytest.mark.parametrize('normalise',[True,False])