diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 97f494947..4756e94c2 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -640,15 +640,25 @@ class Rotation:
                 i = np.where(np.isclose(w,1.0+0.0j))[0][0]
                 ax[0:3] = np.real(vr[0:3,i])
                 diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
-                ax[0:3] = np.where(np.abs(diagDelta)<1.e-6, ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta))
-            return ax
+                ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta))
         else:
             diag_delta = 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]
                                   ])
+            t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,))
             w,vr = np.linalg.eig(om)
-            # TODO ------------------
+            # 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.
+            ax = np.where(np.abs(diag_delta)<0,
+                          np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
+                          np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)) \
+                          * np.abs(diag_delta)*np.sign(-P*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]
+        return ax
+
 
     @staticmethod
     def om2ro(om):
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 875a05d99..6bb0f5160 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -80,7 +80,7 @@ def default():
 
     return [Rotation.fromQuaternion(s) for s in specials] + \
            [Rotation.fromQuaternion(s) for s in specials_scatter] + \
-           [Rotation.fromRandom() for _ in range(n-len(specials))]
+           [Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))]
 
 @pytest.fixture
 def reference_dir(reference_dir_base):
@@ -159,17 +159,22 @@ class TestRotation:
                                            Rotation.qu2ho])
     def test_quaternion_vectorization(self,default,conversion):
         qu = np.array([rot.asQuaternion() for rot in default])
+        #dev_null = conversion(qu.reshape(qu.shape[0]//2,-1,4))
         co = conversion(qu)
         for q,c in zip(qu,co):
-             assert np.allclose(conversion(q),c)
+            print(q,c)
+            assert np.allclose(conversion(q),c)
 
     @pytest.mark.parametrize('conversion',[Rotation.om2eu,
+                                          #Rotation.om2ax,
                                           ])
     def test_matrix_vectorization(self,default,conversion):
-        qu = np.array([rot.asMatrix() for rot in default])
-        co = conversion(qu)
-        for q,c in zip(qu,co):
-             assert np.allclose(conversion(q),c)
+        om = np.array([rot.asMatrix() for rot in default])
+        #dev_null = conversion(om.reshape(om.shape[0]//2,-1,3,3))
+        co = conversion(om)
+        for o,c in zip(om,co):
+            print(o,c)
+            assert np.allclose(conversion(o),c)
 
     @pytest.mark.parametrize('conversion',[Rotation.eu2qu,
                                            Rotation.eu2om,
@@ -177,10 +182,12 @@ class TestRotation:
                                            Rotation.eu2ro,
                                           ])
     def test_Euler_vectorization(self,default,conversion):
-        qu = np.array([rot.asEulers() for rot in default])
-        co = conversion(qu)
-        for q,c in zip(qu,co):
-             assert np.allclose(conversion(q),c)
+        eu = np.array([rot.asEulers() for rot in default])
+        #dev_null = conversion(eu.reshape(eu.shape[0]//2,-1,3))
+        co = conversion(eu)
+        for e,c in zip(eu,co):
+            print(e,c)
+            assert np.allclose(conversion(e),c)
 
     @pytest.mark.parametrize('conversion',[Rotation.ax2qu,
                                            Rotation.ax2om,
@@ -192,4 +199,5 @@ class TestRotation:
         dev_null = conversion(ax.reshape(ax.shape[0]//2,-1,4))
         co = conversion(ax)
         for a,c in zip(ax,co):
-             assert np.allclose(conversion(a),c)
+            print(a,c)
+            assert np.allclose(conversion(a),c)