diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index dcc139b7c..4721143da 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -448,9 +448,12 @@ class Orientation(Rotation,Crystal):
             elif self.family == 'orthorhombic':
                 return (np.prod(1. >= rho_abs,axis=-1)).astype(bool)
             elif self.family == 'monoclinic':
-                return (1. >= rho_abs[...,1]).astype(bool)
+                return np.logical_or(   1. >= rho_abs[...,1],
+                                     np.isnan(rho_abs[...,1]))
+            elif self.family == 'triclinic':
+                return np.ones(rho_abs.shape[:-1]).astype(bool)
             else:
-                return np.all(np.isfinite(rho_abs),axis=-1)
+                raise TypeError(f'unknown symmetry "{self.family}"')
 
 
     @property
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 61529eb46..f93b3b012 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -1372,7 +1372,7 @@ class Rotation:
         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)<1e-12,
+        ax = np.where(np.abs(diag_delta)<1e-13,
                              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))
@@ -1581,14 +1581,13 @@ class Rotation:
     @staticmethod
     def _ho2ax(ho: np.ndarray) -> np.ndarray:
         """Homochoric vector to axis–angle pair."""
-        tfit = np.array([+1.0000000000018852,      -0.5000000002194847,
-                         -0.024999992127593126,    -0.003928701544781374,
-                         -0.0008152701535450438,   -0.0002009500426119712,
-                         -0.00002397986776071756,  -0.00008202868926605841,
-                         +0.00012448715042090092,  -0.0001749114214822577,
-                         +0.0001703481934140054,   -0.00012062065004116828,
-                         +0.000059719705868660826, -0.00001980756723965647,
-                         +0.000003953714684212874, -0.00000036555001439719544])
+        tfit = np.array([+0.9999999999999968,     -0.49999999999986866,     -0.025000000000632055,
+                         -0.003928571496460683,   -0.0008164666077062752,   -0.00019411896443261646,
+                         -0.00004985822229871769, -0.000014164962366386031, -1.9000248160936107e-6,
+                         -5.72184549898506e-6,    +7.772149920658778e-6,    -0.00001053483452909705,
+                         +9.528014229335313e-6,   -5.660288876265125e-6,    +1.2844901692764126e-6,
+                         +1.1255185726258763e-6,  -1.3834391419956455e-6,   +7.513691751164847e-7,
+                         -2.401996891720091e-7,   +4.386887017466388e-8,    -3.5917775353564864e-9])
         hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True)
         s = np.sum(tfit*hmag_squared**np.arange(len(tfit)),axis=-1,keepdims=True)
         with np.errstate(invalid='ignore'):
@@ -1679,7 +1678,7 @@ class Rotation:
 
         """
         with np.errstate(invalid='ignore',divide='ignore'):
-            # get pyramide and scale by grid parameter ratio
+            # get pyramid and scale by grid parameter ratio
             XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc
             order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1])
             q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \
diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py
index 5f9b83e70..2cb97d252 100644
--- a/python/tests/test_Orientation.py
+++ b/python/tests/test_Orientation.py
@@ -224,11 +224,11 @@ class TestOrientation:
 
     @pytest.mark.parametrize('family',crystal_families)
     def test_reduced_corner_cases(self,family):
-        # test whether there is always a sym-eq rotation that falls into the FZ
+        # test whether there is always exactly one sym-eq rotation that falls into the FZ
         N = np.random.randint(10,40)
         size = np.ones(3)*np.pi**(2./3.)
         grid = grid_filters.coordinates0_node([N+1,N+1,N+1],size,-size*.5)
-        evenly_distributed = Orientation.from_cubochoric(x=grid[:-2,:-2,:-2],family=family)
+        evenly_distributed = Orientation.from_cubochoric(x=grid,family=family)
         assert evenly_distributed.shape == evenly_distributed.reduced.shape
 
     @pytest.mark.parametrize('family',crystal_families)
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index c9dfc68b8..8b1bd0de8 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -301,14 +301,13 @@ def ro2ho(ro):
 #---------- Homochoric vector----------
 def ho2ax(ho):
     """Homochoric vector to axis angle pair."""
-    tfit = np.array([+1.0000000000018852,      -0.5000000002194847,
-                     -0.024999992127593126,    -0.003928701544781374,
-                     -0.0008152701535450438,   -0.0002009500426119712,
-                     -0.00002397986776071756,  -0.00008202868926605841,
-                     +0.00012448715042090092,  -0.0001749114214822577,
-                     +0.0001703481934140054,   -0.00012062065004116828,
-                     +0.000059719705868660826, -0.00001980756723965647,
-                     +0.000003953714684212874, -0.00000036555001439719544])
+    tfit = np.array([+0.9999999999999968,     -0.49999999999986866,     -0.025000000000632055,
+                     -0.003928571496460683,   -0.0008164666077062752,   -0.00019411896443261646,
+                     -0.00004985822229871769, -0.000014164962366386031, -1.9000248160936107e-6,
+                     -5.72184549898506e-6,    +7.772149920658778e-6,    -0.00001053483452909705,
+                     +9.528014229335313e-6,   -5.660288876265125e-6,    +1.2844901692764126e-6,
+                     +1.1255185726258763e-6,  -1.3834391419956455e-6,   +7.513691751164847e-7,
+                     -2.401996891720091e-7,   +4.386887017466388e-8,    -3.5917775353564864e-9])
     # normalize h and store the magnitude
     hmag_squared = np.sum(ho**2.)
     if iszero(hmag_squared):