diff --git a/python/damask/_crystal.py b/python/damask/_crystal.py
index 0b1c00458..077d2e13a 100644
--- a/python/damask/_crystal.py
+++ b/python/damask/_crystal.py
@@ -2,10 +2,11 @@ from typing import Union, Dict, List, Tuple
 
 import numpy as np
 
+from ._typehints import FloatSequence, CrystalFamily, CrystalLattice, CrystalKinematics
 from . import util
 from . import Rotation
 
-lattice_symmetries = {
+lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
                 'aP': 'triclinic',
 
                 'mP': 'monoclinic',
@@ -30,9 +31,9 @@ lattice_symmetries = {
 class Crystal():
     """Crystal lattice."""
 
-    def __init__(self,*,
-                 family = None,
-                 lattice = None,
+    def __init__(self, *,
+                 family: CrystalFamily = None,
+                 lattice: CrystalLattice = None,
                  a: float = None, b: float = None, c: float = None,
                  alpha: float = None, beta: float = None, gamma: float = None,
                  degrees: bool = False):
@@ -130,9 +131,8 @@ class Crystal():
             Crystal to check for equality.
 
         """
-        if not isinstance(other, Crystal):
-            return NotImplemented
-        return self.lattice == other.lattice and \
+        return NotImplemented if not isinstance(other, Crystal) else \
+               self.lattice == other.lattice and \
                self.parameters == other.parameters and \
                self.family == other.family
 
@@ -208,7 +208,7 @@ class Crystal():
         ...    }
 
         """
-        _basis  = {
+        _basis: Dict[CrystalFamily, Dict[str, np.ndarray]]  = {
             'cubic':    {'improper':np.array([ [-1.            ,  0.            ,  1. ],
                                                [ np.sqrt(2.)   , -np.sqrt(2.)   ,  0. ],
                                                [ 0.            ,  np.sqrt(3.)   ,  0. ] ]),
@@ -315,19 +315,19 @@ class Crystal():
                                               self.lattice[-1],None),dtype=float)
 
     def to_lattice(self, *,
-                   direction: np.ndarray = None,
-                   plane: np.ndarray = None) -> np.ndarray:
+                   direction: FloatSequence = None,
+                   plane: FloatSequence = None) -> np.ndarray:
         """
         Calculate lattice vector corresponding to crystal frame direction or plane normal.
 
         Parameters
         ----------
-        direction|plane : numpy.ndarray of shape (...,3)
+        direction|plane : numpy.ndarray, shape (...,3)
             Vector along direction or plane normal.
 
         Returns
         -------
-        Miller : numpy.ndarray of shape (...,3)
+        Miller : numpy.ndarray, shape (...,3)
             Lattice vector of direction or plane.
             Use util.scale_to_coprime to convert to (integer) Miller indices.
 
@@ -341,19 +341,19 @@ class Crystal():
 
 
     def to_frame(self, *,
-                 uvw: np.ndarray = None,
-                 hkl: np.ndarray = None) -> np.ndarray:
+                 uvw: FloatSequence = None,
+                 hkl: FloatSequence = None) -> np.ndarray:
         """
         Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
 
         Parameters
         ----------
-        uvw|hkl : numpy.ndarray of shape (...,3)
+        uvw|hkl : numpy.ndarray, shape (...,3)
             Miller indices of crystallographic direction or plane normal.
 
         Returns
         -------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Crystal frame vector along [uvw] direction or (hkl) plane normal.
 
         """
@@ -366,7 +366,7 @@ class Crystal():
 
 
     def kinematics(self,
-                   mode: str) -> Dict[str, List[np.ndarray]]:
+                   mode: CrystalKinematics) -> Dict[str, List[np.ndarray]]:
         """
         Return crystal kinematics systems.
 
@@ -381,7 +381,7 @@ class Crystal():
             Directions and planes of deformation mode families.
 
         """
-        _kinematics = {
+        _kinematics: Dict[CrystalLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
             'cF': {
                 'slip': [np.array([
                            [+0,+1,-1, +1,+1,+1],
@@ -626,7 +626,7 @@ class Crystal():
 
 
     def relation_operations(self,
-                            model: str) -> Tuple[str, Rotation]:
+                            model: str) -> Tuple[CrystalLattice, Rotation]:
         """
         Crystallographic orientation relationships for phase transformations.
 
@@ -658,7 +658,7 @@ class Crystal():
         https://doi.org/10.1016/j.actamat.2004.11.021
 
         """
-        _orientation_relationships = {
+        _orientation_relationships: Dict[str, Dict[CrystalLattice,np.ndarray]] = {
           'KS': {
             'cF' : np.array([
                 [[-1, 0, 1],[ 1, 1, 1]],
diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index b7a43e0fb..c7fbbbb85 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -1,8 +1,10 @@
 import inspect
 import copy
+from typing import Union, Callable, List, Dict, Any, Tuple, TypeVar
 
 import numpy as np
 
+from ._typehints import FloatSequence, IntSequence, CrystalFamily, CrystalLattice
 from . import Rotation
 from . import Crystal
 from . import util
@@ -33,6 +35,7 @@ _parameter_doc = \
 
        """
 
+MyType = TypeVar('MyType', bound='Orientation')
 
 class Orientation(Rotation,Crystal):
     """
@@ -93,12 +96,13 @@ class Orientation(Rotation,Crystal):
 
     @util.extend_docstring(_parameter_doc)
     def __init__(self,
-                 rotation = np.array([1.0,0.0,0.0,0.0]), *,
-                 family = None,
-                 lattice = None,
-                 a = None,b = None,c = None,
-                 alpha = None,beta = None,gamma = None,
-                 degrees = False):
+                 rotation: Union[FloatSequence, Rotation] = np.array([1.,0.,0.,0.]),
+                 *,
+                 family: CrystalFamily = None,
+                 lattice: CrystalLattice = None,
+                 a: float = None, b: float = None, c: float = None,
+                 alpha: float = None, beta: float = None, gamma: float = None,
+                 degrees: bool = False):
         """
         New orientation.
 
@@ -115,13 +119,13 @@ class Orientation(Rotation,Crystal):
                               a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees)
 
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         """Represent."""
         return '\n'.join([Crystal.__repr__(self),
                           Rotation.__repr__(self)])
 
-
-    def __copy__(self,rotation=None):
+    def __copy__(self: MyType,
+                 rotation: Union[FloatSequence, Rotation] = None) -> MyType:
         """Create deep copy."""
         dup = copy.deepcopy(self)
         if rotation is not None:
@@ -131,7 +135,9 @@ class Orientation(Rotation,Crystal):
     copy = __copy__
 
 
-    def __eq__(self,other):
+
+    def __eq__(self,
+               other: object) -> bool:
         """
         Equal to other.
 
@@ -141,12 +147,15 @@ class Orientation(Rotation,Crystal):
             Orientation to check for equality.
 
         """
+        if not isinstance(other, Orientation):
+            return NotImplemented
         matching_type = self.family == other.family and \
                         self.lattice == other.lattice and \
                         self.parameters == other.parameters
         return np.logical_and(matching_type,super(self.__class__,self.reduced).__eq__(other.reduced))
 
-    def __ne__(self,other):
+    def __ne__(self,
+               other: object) -> bool:
         """
         Not equal to other.
 
@@ -156,10 +165,14 @@ class Orientation(Rotation,Crystal):
             Orientation to check for equality.
 
         """
-        return np.logical_not(self==other)
+        return np.logical_not(self==other) if isinstance(other, Orientation) else NotImplemented
 
 
-    def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def isclose(self: MyType,
+                other: MyType,
+                rtol: float = 1e-5,
+                atol: float = 1e-8,
+                equal_nan: bool = True) -> bool:
         """
         Report where values are approximately equal to corresponding ones of other Orientation.
 
@@ -176,7 +189,7 @@ class Orientation(Rotation,Crystal):
 
         Returns
         -------
-        mask : numpy.ndarray bool
+        mask : numpy.ndarray of bool, shape (self.shape)
             Mask indicating where corresponding orientations are close.
 
         """
@@ -187,7 +200,11 @@ class Orientation(Rotation,Crystal):
 
 
 
-    def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def allclose(self: MyType,
+                 other: MyType,
+                 rtol: float = 1e-5,
+                 atol: float = 1e-8,
+                 equal_nan: bool = True) -> bool:
         """
         Test whether all values are approximately equal to corresponding ones of other Orientation.
 
@@ -208,10 +225,11 @@ class Orientation(Rotation,Crystal):
             Whether all values are close between both orientations.
 
         """
-        return np.all(self.isclose(other,rtol,atol,equal_nan))
+        return bool(np.all(self.isclose(other,rtol,atol,equal_nan)))
 
 
-    def __mul__(self,other):
+    def __mul__(self: MyType,
+                other: Union[Rotation, 'Orientation']) -> MyType:
         """
         Compose this orientation with other.
 
@@ -226,14 +244,15 @@ class Orientation(Rotation,Crystal):
             Compound rotation self*other, i.e. first other then self rotation.
 
         """
-        if isinstance(other,Orientation) or isinstance(other,Rotation):
-            return self.copy(rotation=Rotation.__mul__(self,Rotation(other.quaternion)))
+        if isinstance(other, (Orientation,Rotation)):
+            return self.copy(Rotation(self.quaternion)*Rotation(other.quaternion))
         else:
             raise TypeError('use "O@b", i.e. matmul, to apply Orientation "O" to object "b"')
 
 
     @staticmethod
-    def _split_kwargs(kwargs,target):
+    def _split_kwargs(kwargs: Dict[str, Any],
+                      target: Callable) -> Tuple[Dict[str, Any], ...]:
         """
         Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects.
 
@@ -252,7 +271,7 @@ class Orientation(Rotation,Crystal):
             Valid keyword arguments of Orientation object.
 
         """
-        kws = ()
+        kws: Tuple[Dict[str, Any], ...] = ()
         for t in (target,Orientation.__init__):
             kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},)
 
@@ -264,105 +283,108 @@ class Orientation(Rotation,Crystal):
 
 
     @classmethod
-    @util.extended_docstring(Rotation.from_random,_parameter_doc)
-    def from_random(cls,**kwargs):
+    @util.extended_docstring(Rotation.from_random, _parameter_doc)
+    def from_random(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random)
         return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_quaternion,_parameter_doc)
-    def from_quaternion(cls,**kwargs):
+    def from_quaternion(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion)
         return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc)
-    def from_Euler_angles(cls,**kwargs):
+    def from_Euler_angles(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles)
         return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_axis_angle,_parameter_doc)
-    def from_axis_angle(cls,**kwargs):
+    def from_axis_angle(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle)
         return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_basis,_parameter_doc)
-    def from_basis(cls,**kwargs):
+    def from_basis(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis)
         return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_matrix,_parameter_doc)
-    def from_matrix(cls,**kwargs):
+    def from_matrix(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix)
         return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc)
-    def from_Rodrigues_vector(cls,**kwargs):
+    def from_Rodrigues_vector(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector)
         return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_homochoric,_parameter_doc)
-    def from_homochoric(cls,**kwargs):
+    def from_homochoric(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric)
         return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_cubochoric,_parameter_doc)
-    def from_cubochoric(cls,**kwargs):
+    def from_cubochoric(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric)
         return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_spherical_component,_parameter_doc)
-    def from_spherical_component(cls,**kwargs):
+    def from_spherical_component(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component)
         return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extended_docstring(Rotation.from_fiber_component,_parameter_doc)
-    def from_fiber_component(cls,**kwargs):
+    def from_fiber_component(cls, **kwargs) -> 'Orientation':
         kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component)
         return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori)
 
 
     @classmethod
     @util.extend_docstring(_parameter_doc)
-    def from_directions(cls,uvw,hkl,**kwargs):
+    def from_directions(cls,
+                        uvw: FloatSequence,
+                        hkl: FloatSequence,
+                        **kwargs) -> 'Orientation':
         """
         Initialize orientation object from two crystallographic directions.
 
         Parameters
         ----------
-        uvw : list, numpy.ndarray of shape (...,3)
-            lattice direction aligned with lab frame x-direction.
-        hkl : list, numpy.ndarray of shape (...,3)
-            lattice plane normal aligned with lab frame z-direction.
+        uvw : numpy.ndarray, shape (...,3)
+            Lattice direction aligned with lab frame x-direction.
+        hkl : numpy.ndarray, shape (...,3)
+            Lattice plane normal aligned with lab frame z-direction.
 
         """
         o = cls(**kwargs)
         x = o.to_frame(uvw=uvw)
         z = o.to_frame(hkl=hkl)
         om = np.stack([x,np.cross(z,x),z],axis=-2)
-        return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
+        return o.copy(Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
 
 
     @property
-    def equivalent(self):
+    def equivalent(self: MyType) -> MyType:
         """
         Orientations that are symmetrically equivalent.
 
@@ -372,11 +394,11 @@ class Orientation(Rotation,Crystal):
         """
         sym_ops = self.symmetry_operations
         o = sym_ops.broadcast_to(sym_ops.shape+self.shape,mode='right')
-        return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
+        return self.copy(o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
 
 
     @property
-    def reduced(self):
+    def reduced(self: MyType) -> MyType:
         """Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry."""
         eq   = self.equivalent
         ok   = eq.in_FZ
@@ -387,13 +409,13 @@ class Orientation(Rotation,Crystal):
 
 
     @property
-    def in_FZ(self):
+    def in_FZ(self) -> Union[np.bool_, np.ndarray]:
         """
         Check whether orientation falls into fundamental zone of own symmetry.
 
         Returns
         -------
-        in : numpy.ndarray of bool, quaternion.shape
+        in : numpy.ndarray of bool, shape (self.shape)
             Whether Rodrigues-Frank vector falls into fundamental zone.
 
         Notes
@@ -431,13 +453,13 @@ class Orientation(Rotation,Crystal):
 
 
     @property
-    def in_disorientation_FZ(self):
+    def in_disorientation_FZ(self) -> np.ndarray:
         """
         Check whether orientation falls into fundamental zone of disorientations.
 
         Returns
         -------
-        in : numpy.ndarray of bool, quaternion.shape
+        in : numpy.ndarray of bool, shape (self.shape)
             Whether Rodrigues-Frank vector falls into disorientation FZ.
 
         References
@@ -471,8 +493,9 @@ class Orientation(Rotation,Crystal):
             else:
                 return np.ones_like(rho[...,0],dtype=bool)
 
-
-    def disorientation(self,other,return_operators=False):
+    def disorientation(self,
+                       other: 'Orientation',
+                       return_operators: bool = False) -> object:
         """
         Calculate disorientation between myself and given other orientation.
 
@@ -490,7 +513,7 @@ class Orientation(Rotation,Crystal):
         -------
         disorientation : Orientation
             Disorientation between self and other.
-        operators : numpy.ndarray int of shape (...,2), conditional
+        operators : numpy.ndarray of int, shape (...,2), conditional
             Index of symmetrically equivalent orientation that rotated vector to the SST.
 
         Notes
@@ -530,7 +553,7 @@ class Orientation(Rotation,Crystal):
             raise NotImplementedError('disorientation between different crystal families')
 
         blend = util.shapeblender(self.shape,other.shape)
-        s =  self.equivalent
+        s = self.equivalent
         o = other.equivalent
 
         s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right')
@@ -557,13 +580,15 @@ class Orientation(Rotation,Crystal):
                )
 
 
-    def average(self,weights=None,return_cloud=False):
+    def average(self,
+                weights: FloatSequence = None,
+                return_cloud: bool = False):
         """
         Return orientation average over last dimension.
 
         Parameters
         ----------
-        weights : numpy.ndarray, optional
+        weights : numpy.ndarray, shape (self.shape), optional
             Relative weights of orientations.
         return_cloud : bool, optional
             Return the set of symmetrically equivalent orientations that was used in averaging.
@@ -583,31 +608,30 @@ class Orientation(Rotation,Crystal):
 
         """
         eq = self.equivalent
-        m  = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,))
-                                          .broadcast_to(eq.shape))\
-               .as_axis_angle()[...,3]
+        m  = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,))             # type: ignore
+                                          .broadcast_to(eq.shape)).as_axis_angle()[...,3] # type: ignore
         r = Rotation(np.squeeze(np.take_along_axis(eq.quaternion,
                                                    np.argmin(m,axis=0)[np.newaxis,...,np.newaxis],
                                                    axis=0),
                                 axis=0))
-        return (
-                (self.copy(rotation=Rotation(r).average(weights)),
-                 self.copy(rotation=Rotation(r)))
-                if return_cloud else
-                self.copy(rotation=Rotation(r).average(weights))
+        return ((self.copy(Rotation(r).average(weights)),self.copy(Rotation(r))) if return_cloud else
+                self.copy(Rotation(r).average(weights))
                )
 
 
-    def to_SST(self,vector,proper=False,return_operators=False):
+    def to_SST(self,
+               vector: FloatSequence,
+               proper: bool = False,
+               return_operators: bool = False) -> np.ndarray:
         """
         Rotate vector to ensure it falls into (improper or proper) standard stereographic triangle of crystal symmetry.
 
         Parameters
         ----------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Lab frame vector to align with crystal frame direction.
             Shape of vector blends with shape of own rotation array.
-            For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
+            For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
         proper : bool, optional
             Consider only vectors with z >= 0, hence combine two neighboring SSTs.
             Defaults to False.
@@ -617,15 +641,18 @@ class Orientation(Rotation,Crystal):
 
         Returns
         -------
-        vector_SST : numpy.ndarray of shape (...,3)
+        vector_SST : numpy.ndarray, shape (...,3)
             Rotated vector falling into SST.
-        operators : numpy.ndarray int of shape (...), conditional
+        operators : numpy.ndarray of int, shape (...), conditional
             Index of symmetrically equivalent orientation that rotated vector to SST.
 
         """
+        vector_ = np.array(vector,float)
+        if vector_.shape[-1] != 3:
+            raise ValueError('input is not a field of three-dimensional vectors')
         eq  = self.equivalent
-        blend = util.shapeblender(eq.shape,np.array(vector).shape[:-1])
-        poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(np.array(vector),blend+(3,))
+        blend = util.shapeblender(eq.shape,vector_.shape[:-1])
+        poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(vector_,blend+(3,))
         ok    = self.in_SST(poles,proper=proper)
         ok   &= np.cumsum(ok,axis=0) == 1
         loc   = np.where(ok)
@@ -637,13 +664,15 @@ class Orientation(Rotation,Crystal):
                )
 
 
-    def in_SST(self,vector,proper=False):
+    def in_SST(self,
+               vector: FloatSequence,
+               proper: bool = False) -> Union[np.bool_, np.ndarray]:
         """
         Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry.
 
         Parameters
         ----------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Vector to check.
         proper : bool, optional
             Consider only vectors with z >= 0, hence combine two neighboring SSTs.
@@ -655,39 +684,43 @@ class Orientation(Rotation,Crystal):
             Whether vector falls into SST.
 
         """
-        if not isinstance(vector,np.ndarray) or vector.shape[-1] != 3:
+        vector_ = np.array(vector,float)
+        if vector_.shape[-1] != 3:
             raise ValueError('input is not a field of three-dimensional vectors')
 
         if self.standard_triangle is None:                                                          # direct exit for no symmetry
-            return  np.ones_like(vector[...,0],bool)
+            return np.ones_like(vector_[...,0],bool)
 
         if proper:
             components_proper   = np.around(np.einsum('...ji,...i',
-                                                      np.broadcast_to(self.standard_triangle['proper'], vector.shape+(3,)),
-                                                      vector), 12)
+                                                      np.broadcast_to(self.standard_triangle['proper'],   vector_.shape+(3,)),
+                                                      vector_), 12)
             components_improper = np.around(np.einsum('...ji,...i',
-                                                      np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
-                                                      vector), 12)
+                                                      np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
+                                                      vector_), 12)
             return   np.all(components_proper   >= 0.0,axis=-1) \
                    | np.all(components_improper >= 0.0,axis=-1)
         else:
             components = np.around(np.einsum('...ji,...i',
-                                             np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
-                                             np.block([vector[...,:2],np.abs(vector[...,2:3])])), 12)
+                                             np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
+                                             np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
 
             return np.all(components >= 0.0,axis=-1)
 
 
-    def IPF_color(self,vector,in_SST=True,proper=False):
+    def IPF_color(self,
+                  vector: FloatSequence,
+                  in_SST: bool = True,
+                  proper: bool = False) -> np.ndarray:
         """
         Map vector to RGB color within standard stereographic triangle of own symmetry.
 
         Parameters
         ----------
-        vector : numpy.ndarray of shape (...,3)
+        vector : numpy.ndarray, shape (...,3)
             Vector to colorize.
             Shape of vector blends with shape of own rotation array.
-            For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
+            For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
         in_SST : bool, optional
             Consider symmetrically equivalent orientations such that poles are located in SST.
             Defaults to True.
@@ -697,7 +730,7 @@ class Orientation(Rotation,Crystal):
 
         Returns
         -------
-        rgb : numpy.ndarray of shape (...,3)
+        rgb : numpy.ndarray, shape (...,3)
            RGB array of IPF colors.
 
         Examples
@@ -721,35 +754,35 @@ class Orientation(Rotation,Crystal):
 
         if proper:
             components_proper   = np.around(np.einsum('...ji,...i',
-                                                      np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
+                                                      np.broadcast_to(self.standard_triangle['proper'],   vector_.shape+(3,)),
                                                       vector_), 12)
             components_improper = np.around(np.einsum('...ji,...i',
                                                       np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
                                                       vector_), 12)
-            in_SST = np.all(components_proper   >= 0.0,axis=-1) \
-                   | np.all(components_improper >= 0.0,axis=-1)
-            components = np.where((in_SST & np.all(components_proper   >= 0.0,axis=-1))[...,np.newaxis],
+            in_SST_ = np.all(components_proper   >= 0.0,axis=-1) \
+                    | np.all(components_improper >= 0.0,axis=-1)
+            components = np.where((in_SST_ & np.all(components_proper   >= 0.0,axis=-1))[...,np.newaxis],
                                   components_proper,components_improper)
         else:
             components = np.around(np.einsum('...ji,...i',
                                              np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)),
                                              np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
 
-            in_SST = np.all(components >= 0.0,axis=-1)
+            in_SST_ = np.all(components >= 0.0,axis=-1)
 
         with np.errstate(invalid='ignore',divide='ignore'):
             rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5                # smoothen color ramps
             rgb = np.clip(rgb,0.,1.)                                                                # clip intensity
             rgb /= np.max(rgb,axis=-1,keepdims=True)                                                # normalize to (HS)V = 1
-        rgb[np.broadcast_to(~in_SST[...,np.newaxis],rgb.shape)] = 0.0
+        rgb[np.broadcast_to(~in_SST_[...,np.newaxis],rgb.shape)] = 0.0
 
         return rgb
 
 
     @property
-    def symmetry_operations(self):
+    def symmetry_operations(self) -> Rotation:
         """Symmetry operations as Rotations."""
-        _symmetry_operations = {
+        _symmetry_operations: Dict[CrystalFamily, List]  = {
             'cubic':         [
                               [ 1.0,            0.0,            0.0,            0.0            ],
                               [ 0.0,            1.0,            0.0,            0.0            ],
@@ -819,22 +852,25 @@ class Orientation(Rotation,Crystal):
 ####################################################################################################
     # functions that require lattice, not just family
 
-    def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False):
+    def to_pole(self, *,
+                uvw: FloatSequence = None,
+                hkl: FloatSequence = None,
+                with_symmetry: bool = False) -> np.ndarray:
         """
         Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl).
 
         Parameters
         ----------
-        uvw|hkl : numpy.ndarray of shape (...,3)
+        uvw|hkl : numpy.ndarray, shape (...,3)
             Miller indices of crystallographic direction or plane normal.
             Shape of vector blends with shape of own rotation array.
-            For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs.
+            For example, a rotation array, shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
         with_symmetry : bool, optional
             Calculate all N symmetrically equivalent vectors.
 
         Returns
         -------
-        vector : numpy.ndarray of shape (...,3) or (...,N,3)
+        vector : numpy.ndarray, shape (...,3) or (...,N,3)
             Lab frame vector (or vectors if with_symmetry) along [uvw] direction or (hkl) plane normal.
 
         """
@@ -846,23 +882,24 @@ class Orientation(Rotation,Crystal):
             blend += sym_ops.shape
             v = sym_ops.broadcast_to(shape) \
               @ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,))
-        return ~(self.broadcast_to(blend)) \
-               @ np.broadcast_to(v,blend+(3,))
+        return ~(self.broadcast_to(blend))@ np.broadcast_to(v,blend+(3,))
 
 
-    def Schmid(self,*,N_slip=None,N_twin=None):
+    def Schmid(self, *,
+               N_slip: IntSequence = None,
+               N_twin: IntSequence = None) -> np.ndarray:
         u"""
         Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems.
 
         Parameters
         ----------
-        N_slip|N_twin : iterable of int
+        N_slip|N_twin : '*' or iterable of int
             Number of deformation systems per family of the deformation system.
             Use '*' to select all.
 
         Returns
         -------
-        P : numpy.ndarray of shape (N,...,3,3)
+        P : numpy.ndarray, shape (N,...,3,3)
             Schmid matrix for each of the N deformation systems.
 
         Examples
@@ -887,6 +924,8 @@ class Orientation(Rotation,Crystal):
                             (self.kinematics('twin'),N_twin)
         if active == '*': active = [len(a) for a in kinematics['direction']]
 
+        if not active:
+            raise ValueError('Schmid matrix not defined')
         d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)]))
         p = self.to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)]))
         P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True),
@@ -897,7 +936,8 @@ class Orientation(Rotation,Crystal):
                @ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape)
 
 
-    def related(self,model):
+    def related(self: MyType,
+                model: str) -> MyType:
         """
         Orientations derived from the given relationship.
 
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 8e5f7f7f9..9c6bde2dc 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -1,7 +1,9 @@
 import copy
+from typing import Union, Sequence, Tuple, Literal, List, TypeVar
 
 import numpy as np
 
+from ._typehints import FloatSequence, IntSequence, NumpyRngSeed
 from . import tensor
 from . import util
 from . import grid_filters
@@ -13,6 +15,8 @@ _sc   = np.pi**(1./6.)/6.**(1./6.)
 _beta = np.pi**(5./6.)/6.**(1./6.)/2.
 _R1   = (3.*np.pi/4.)**(1./3.)
 
+MyType = TypeVar('MyType', bound='Rotation')
+
 class Rotation:
     u"""
     Rotation with functionality for conversion between different representations.
@@ -61,7 +65,8 @@ class Rotation:
 
     __slots__ = ['quaternion']
 
-    def __init__(self,rotation = np.array([1.0,0.0,0.0,0.0])):
+    def __init__(self,
+                 rotation: Union[FloatSequence, 'Rotation'] = np.array([1.0,0.0,0.0,0.0])):
         """
         New rotation.
 
@@ -73,6 +78,7 @@ class Rotation:
             Defaults to no rotation.
 
         """
+        self.quaternion: np.ndarray
         if isinstance(rotation,Rotation):
             self.quaternion = rotation.quaternion.copy()
         elif np.array(rotation).shape[-1] == 4:
@@ -81,13 +87,14 @@ class Rotation:
             raise TypeError('"rotation" is neither a Rotation nor a quaternion')
 
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         """Represent rotation as unit quaternion(s)."""
         return f'Quaternion{" " if self.quaternion.shape == (4,) else "s of shape "+str(self.quaternion.shape[:-1])+chr(10)}'\
                + str(self.quaternion)
 
 
-    def __copy__(self,rotation=None):
+    def __copy__(self: MyType,
+                 rotation: Union[FloatSequence, 'Rotation'] = None) -> MyType:
         """Create deep copy."""
         dup = copy.deepcopy(self)
         if rotation is not None:
@@ -97,14 +104,15 @@ class Rotation:
     copy = __copy__
 
 
-    def __getitem__(self,item):
+    def __getitem__(self,
+                    item: Union[Tuple[int], int, bool, np.bool_, np.ndarray]):
         """Return slice according to item."""
-        return self.copy() \
-               if self.shape == () else \
-               self.copy(rotation=self.quaternion[item+(slice(None),)] if isinstance(item,tuple) else self.quaternion[item])
+        return self.copy() if self.shape == () else \
+               self.copy(self.quaternion[item+(slice(None),)] if isinstance(item,tuple) else self.quaternion[item])
 
 
-    def __eq__(self,other):
+    def __eq__(self,
+               other: object) -> bool:
         """
         Equal to other.
 
@@ -114,11 +122,13 @@ class Rotation:
             Rotation to check for equality.
 
         """
-        return np.logical_or(np.all(self.quaternion ==      other.quaternion,axis=-1),
+        return NotImplemented if not isinstance(other, Rotation) else \
+               np.logical_or(np.all(self.quaternion ==      other.quaternion,axis=-1),
                              np.all(self.quaternion == -1.0*other.quaternion,axis=-1))
 
 
-    def __ne__(self,other):
+    def __ne__(self,
+               other: object) -> bool:
         """
         Not equal to other.
 
@@ -128,10 +138,13 @@ class Rotation:
             Rotation to check for inequality.
 
         """
-        return np.logical_not(self==other)
+        return np.logical_not(self==other) if isinstance(other, Rotation) else NotImplemented
 
-
-    def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def isclose(self: MyType,
+                other: MyType,
+                rtol: float = 1e-5,
+                atol: float = 1e-8,
+                equal_nan: bool = True) -> bool:
         """
         Report where values are approximately equal to corresponding ones of other Rotation.
 
@@ -148,7 +161,7 @@ class Rotation:
 
         Returns
         -------
-        mask : numpy.ndarray bool
+        mask : numpy.ndarray of bool, shape (self.shape)
             Mask indicating where corresponding rotations are close.
 
         """
@@ -158,7 +171,11 @@ class Rotation:
                              np.all(np.isclose(s,-1.0*o,rtol,atol,equal_nan),axis=-1))
 
 
-    def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def allclose(self: MyType,
+                 other: MyType,
+                 rtol: float = 1e-5,
+                 atol: float = 1e-8,
+                 equal_nan: bool = True) -> Union[np.bool_, bool]:
         """
         Test whether all values are approximately equal to corresponding ones of other Rotation.
 
@@ -188,27 +205,28 @@ class Rotation:
 
 
     @property
-    def size(self):
+    def size(self) -> int:
         return self.quaternion[...,0].size
 
     @property
-    def shape(self):
+    def shape(self) -> Tuple[int, ...]:
         return self.quaternion[...,0].shape
 
 
-    def __len__(self):
+    def __len__(self) -> int:
         """Length of leading/leftmost dimension of array."""
         return 0 if self.shape == () else self.shape[0]
 
 
-    def __invert__(self):
+    def __invert__(self: MyType) -> MyType:
         """Inverse rotation (backward rotation)."""
         dup = self.copy()
         dup.quaternion[...,1:] *= -1
         return dup
 
 
-    def __pow__(self,exp):
+    def __pow__(self: MyType,
+                exp: Union[float, int]) -> MyType:
         """
         Perform the rotation 'exp' times.
 
@@ -220,9 +238,10 @@ class Rotation:
         """
         phi = np.arccos(self.quaternion[...,0:1])
         p = self.quaternion[...,1:]/np.linalg.norm(self.quaternion[...,1:],axis=-1,keepdims=True)
-        return self.copy(rotation=Rotation(np.block([np.cos(exp*phi),np.sin(exp*phi)*p]))._standardize())
+        return self.copy(Rotation(np.block([np.cos(exp*phi),np.sin(exp*phi)*p]))._standardize())
 
-    def __ipow__(self,exp):
+    def __ipow__(self: MyType,
+                 exp: Union[float, int]) -> MyType:
         """
         Perform the rotation 'exp' times (in-place).
 
@@ -235,13 +254,14 @@ class Rotation:
         return self**exp
 
 
-    def __mul__(self,other):
+    def __mul__(self: MyType,
+                other: MyType) -> MyType:
         """
         Compose with other.
 
         Parameters
         ----------
-        other : Rotation of shape (self.shape)
+        other : Rotation, shape (self.shape)
             Rotation for composition.
 
         Returns
@@ -257,30 +277,32 @@ class Rotation:
             p_o = other.quaternion[...,1:]
             q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,)))
             p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
-            return Rotation(np.block([q,p]))._standardize()
+            return self.copy(Rotation(np.block([q,p]))._standardize())
         else:
             raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
 
-    def __imul__(self,other):
+    def __imul__(self: MyType,
+                 other: MyType) -> MyType:
         """
         Compose with other (in-place).
 
         Parameters
         ----------
-        other : Rotation of shape (self.shape)
+        other : Rotation, shape (self.shape)
             Rotation for composition.
 
         """
         return self*other
 
 
-    def __truediv__(self,other):
+    def __truediv__(self: MyType,
+                    other: MyType) -> MyType:
         """
         Compose with inverse of other.
 
         Parameters
         ----------
-        other : damask.Rotation of shape (self.shape)
+        other : damask.Rotation, shape (self.shape)
             Rotation to invert for composition.
 
         Returns
@@ -294,35 +316,37 @@ class Rotation:
         else:
             raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
 
-    def __itruediv__(self,other):
+    def __itruediv__(self: MyType,
+                     other: MyType) -> MyType:
         """
         Compose with inverse of other (in-place).
 
         Parameters
         ----------
-        other : Rotation of shape (self.shape)
+        other : Rotation, shape (self.shape)
             Rotation to invert for composition.
 
         """
         return self/other
 
 
-    def __matmul__(self,other):
+    def __matmul__(self,
+                   other: np.ndarray) -> np.ndarray:
         """
         Rotate vector, second order tensor, or fourth order tensor.
 
         Parameters
         ----------
-        other : numpy.ndarray of shape (...,3), (...,3,3), or (...,3,3,3,3)
+        other : numpy.ndarray, shape (...,3), (...,3,3), or (...,3,3,3,3)
             Vector or tensor on which to apply the rotation.
 
         Returns
         -------
-        rotated : numpy.ndarray of shape (...,3), (...,3,3), or (...,3,3,3,3)
+        rotated : numpy.ndarray, shape (...,3), (...,3,3), or (...,3,3,3,3)
             Rotated vector or tensor, i.e. transformed to frame defined by rotation.
 
         """
-        if isinstance(other,np.ndarray):
+        if isinstance(other, np.ndarray):
             if self.shape + (3,) == other.shape:
                 q_m = self.quaternion[...,0]
                 p_m = self.quaternion[...,1:]
@@ -342,7 +366,7 @@ class Rotation:
                 return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other)
             else:
                 raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors')
-        elif isinstance(other,Rotation):
+        elif isinstance(other, Rotation):
             raise TypeError('Use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"')
         else:
             raise TypeError(f'Cannot rotate {type(other)}')
@@ -350,42 +374,67 @@ class Rotation:
     apply = __matmul__
 
 
-    def _standardize(self):
+    def _standardize(self: MyType) -> MyType:
         """Standardize quaternion (ensure positive real hemisphere)."""
         self.quaternion[self.quaternion[...,0] < 0.0] *= -1
         return self
 
 
-    def append(self,other):
+    def append(self: MyType,
+               other: Union[MyType, List[MyType]]) -> MyType:
         """
         Extend array along first dimension with other array(s).
 
         Parameters
         ----------
-        other : damask.Rotation
+        other : (list of) damask.Rotation
 
         """
-        return self.copy(rotation=np.vstack(tuple(map(lambda x:x.quaternion,
-                                                      [self]+other if isinstance(other,list) else [self,other]))))
+        return self.copy(np.vstack(tuple(map(lambda x:x.quaternion,
+                                             [self]+other if isinstance(other,list) else [self,other]))))
 
 
-    def flatten(self,order = 'C'):
+    def flatten(self: MyType,
+                order: Literal['C','F','A'] = 'C') -> MyType:
         """
         Flatten array.
 
+        Parameters
+        ----------
+        order : {'C', 'F', 'A'}, optional
+            'C' flattens in row-major (C-style) order.
+            'F' flattens in column-major (Fortran-style) order.
+            'A' flattens in column-major order if object is Fortran contiguous in memory,
+            row-major order otherwise.
+            Defaults to 'C'.
+
         Returns
         -------
         flattened : damask.Rotation
             Rotation flattened to single dimension.
 
         """
-        return self.copy(rotation=self.quaternion.reshape((-1,4),order=order))
+        return self.copy(self.quaternion.reshape((-1,4),order=order))
 
 
-    def reshape(self,shape,order = 'C'):
+    def reshape(self: MyType,
+                shape: Union[int, Tuple[int, ...]],
+                order: Literal['C','F','A'] = 'C') -> MyType:
         """
         Reshape array.
 
+        Parameters
+        ----------
+        shape : int or tuple of ints
+            The new shape should be compatible with the original shape.
+            If an integer is supplied, then the result will be a 1-D array of that length.
+        order : {'C', 'F', 'A'}, optional
+            'C' flattens in row-major (C-style) order.
+            'F' flattens in column-major (Fortran-style) order.
+            'A' flattens in column-major order if object is Fortran contiguous in memory,
+            row-major order otherwise.
+            Defaults to 'C'.
+
         Returns
         -------
         reshaped : damask.Rotation
@@ -393,16 +442,18 @@ class Rotation:
 
         """
         if isinstance(shape,(int,np.integer)): shape = (shape,)
-        return self.copy(rotation=self.quaternion.reshape(tuple(shape)+(4,),order=order))
+        return self.copy(self.quaternion.reshape(tuple(shape)+(4,),order=order))
 
 
-    def broadcast_to(self,shape,mode = 'right'):
+    def broadcast_to(self: MyType,
+                     shape: Union[int, Tuple[int, ...]],
+                     mode: Literal['left', 'right'] = 'right') -> MyType:
         """
         Broadcast array.
 
         Parameters
         ----------
-        shape : tuple
+        shape : int or tuple of ints
             Shape of broadcasted array.
         mode : str, optional
             Where to preferentially locate missing dimensions.
@@ -415,17 +466,18 @@ class Rotation:
 
         """
         if isinstance(shape,(int,np.integer)): shape = (shape,)
-        return self.copy(rotation=np.broadcast_to(self.quaternion.reshape(util.shapeshifter(self.shape,shape,mode)+(4,)),
+        return self.copy(np.broadcast_to(self.quaternion.reshape(util.shapeshifter(self.shape,shape,mode)+(4,)),
                                                   shape+(4,)))
 
 
-    def average(self,weights = None):
+    def average(self: MyType,
+                weights: FloatSequence = None) -> MyType:
         """
         Average along last array dimension.
 
         Parameters
         ----------
-        weights : list of floats, optional
+        weights : numpy.ndarray, shape (self.shape), optional
             Relative weight of each rotation.
 
         Returns
@@ -443,22 +495,22 @@ class Rotation:
             """Intermediate representation supporting quaternion averaging."""
             return np.einsum('...i,...j',quat,quat)
 
-        if weights is None:
-            weights = np.ones(self.shape,dtype=float)
+        weights_ = np.ones(self.shape,dtype=float) if weights is None else np.array(weights,float)
 
-        eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights[...,np.newaxis,np.newaxis],axis=-3) \
-                                /np.sum(                      weights[...,np.newaxis,np.newaxis],axis=-3))
+        eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights_[...,np.newaxis,np.newaxis],axis=-3) \
+                                /np.sum(                      weights_[...,np.newaxis,np.newaxis],axis=-3))
 
-        return Rotation.from_quaternion(np.real(
-                                        np.squeeze(
-                                        np.take_along_axis(vec,
-                                                           eig.argmax(axis=-1)[...,np.newaxis,np.newaxis],
-                                                           axis=-1),
-                                        axis=-1)),
-                                        accept_homomorph = True)
+        return self.copy(Rotation.from_quaternion(np.real(
+                                                  np.squeeze(
+                                                  np.take_along_axis(vec,
+                                                                     eig.argmax(axis=-1)[...,np.newaxis,np.newaxis],
+                                                                     axis=-1),
+                                                  axis=-1)),
+                                                  accept_homomorph = True))
 
 
-    def misorientation(self,other):
+    def misorientation(self: MyType,
+                       other: MyType) -> MyType:
         """
         Calculate misorientation to other Rotation.
 
@@ -479,20 +531,20 @@ class Rotation:
     ################################################################################################
     # convert to different orientation representations (numpy arrays)
 
-    def as_quaternion(self):
+    def as_quaternion(self) -> np.ndarray:
         """
         Represent as unit quaternion.
 
         Returns
         -------
-        q : numpy.ndarray of shape (...,4)
+        q : numpy.ndarray, shape (...,4)
             Unit quaternion (q_0, q_1, q_2, q_3) in positive real hemisphere, i.e. ǀqǀ = 1, q_0 ≥ 0.
 
         """
         return self.quaternion.copy()
 
     def as_Euler_angles(self,
-                        degrees = False):
+                        degrees: bool = False) -> np.ndarray:
         """
         Represent as Bunge Euler angles.
 
@@ -503,7 +555,7 @@ class Rotation:
 
         Returns
         -------
-        phi : numpy.ndarray of shape (...,3)
+        phi : numpy.ndarray, shape (...,3)
             Bunge Euler angles (φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π])
             or (φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]) if degrees == True.
 
@@ -522,12 +574,11 @@ class Rotation:
 
         """
         eu = Rotation._qu2eu(self.quaternion)
-        if degrees: eu = np.degrees(eu)
-        return eu
+        return np.degrees(eu) if degrees else eu
 
     def as_axis_angle(self,
-                      degrees = False,
-                      pair = False):
+                      degrees: bool = False,
+                      pair: bool = False) -> Union[Tuple[np.ndarray, np.ndarray], np.ndarray]:
         """
         Represent as axis–angle pair.
 
@@ -540,7 +591,7 @@ class Rotation:
 
         Returns
         -------
-        axis_angle : numpy.ndarray of shape (...,4) or tuple ((...,3), (...)) if pair == True
+        axis_angle : numpy.ndarray, shape (...,4) or tuple ((...,3), (...)) if pair == True
             Axis and angle [n_1, n_2, n_3, ω] with ǀnǀ = 1 and ω ∈ [0,π]
             or ω ∈ [0,180] if degrees == True.
 
@@ -554,17 +605,17 @@ class Rotation:
         (array([0., 0., 1.]), array(0.))
 
         """
-        ax = Rotation._qu2ax(self.quaternion)
+        ax: np.ndarray = Rotation._qu2ax(self.quaternion)
         if degrees: ax[...,3] = np.degrees(ax[...,3])
         return (ax[...,:3],ax[...,3]) if pair else ax
 
-    def as_matrix(self):
+    def as_matrix(self) -> np.ndarray:
         """
         Represent as rotation matrix.
 
         Returns
         -------
-        R : numpy.ndarray of shape (...,3,3)
+        R : numpy.ndarray, shape (...,3,3)
             Rotation matrix R with det(R) = 1, R.T ∙ R = I.
 
         Examples
@@ -582,7 +633,7 @@ class Rotation:
         return Rotation._qu2om(self.quaternion)
 
     def as_Rodrigues_vector(self,
-                            compact = False):
+                            compact: bool = False) -> np.ndarray:
         """
         Represent as Rodrigues–Frank vector with separate axis and angle argument.
 
@@ -594,7 +645,7 @@ class Rotation:
 
         Returns
         -------
-        rho : numpy.ndarray of shape (...,4) or (...,3) if compact == True
+        rho : numpy.ndarray, shape (...,4) or (...,3) if compact == True
             Rodrigues–Frank vector [n_1, n_2, n_3, tan(ω/2)] with ǀnǀ = 1 and ω ∈ [0,π]
             or [n_1, n_2, n_3] with ǀnǀ = tan(ω/2) and ω ∈ [0,π] if compact == True.
 
@@ -615,13 +666,13 @@ class Rotation:
         else:
             return ro
 
-    def as_homochoric(self):
+    def as_homochoric(self) -> np.ndarray:
         """
         Represent as homochoric vector.
 
         Returns
         -------
-        h : numpy.ndarray of shape (...,3)
+        h : numpy.ndarray, shape (...,3)
             Homochoric vector (h_1, h_2, h_3) with ǀhǀ < (3/4*π)^(1/3).
 
         Examples
@@ -636,13 +687,13 @@ class Rotation:
         """
         return Rotation._qu2ho(self.quaternion)
 
-    def as_cubochoric(self):
+    def as_cubochoric(self) -> np.ndarray:
         """
         Represent as cubochoric vector.
 
         Returns
         -------
-        x : numpy.ndarray of shape (...,3)
+        x : numpy.ndarray, shape (...,3)
               Cubochoric vector (x_1, x_2, x_3) with max(x_i) < 1/2*π^(2/3).
 
         Examples
@@ -661,15 +712,15 @@ class Rotation:
     # Static constructors. The input data needs to follow the conventions, options allow to
     # relax the conventions.
     @staticmethod
-    def from_quaternion(q,
-                        accept_homomorph = False,
-                        P = -1):
+    def from_quaternion(q: Union[Sequence[FloatSequence], np.ndarray],
+                        accept_homomorph: bool = False,
+                        P: Literal[1, -1] = -1) -> 'Rotation':
         """
         Initialize from quaternion.
 
         Parameters
         ----------
-        q : numpy.ndarray of shape (...,4)
+        q : numpy.ndarray, shape (...,4)
             Unit quaternion (q_0, q_1, q_2, q_3) in positive real hemisphere, i.e. ǀqǀ = 1, q_0 ≥ 0.
         accept_homomorph : bool, optional
             Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere).
@@ -696,14 +747,14 @@ class Rotation:
         return Rotation(qu)
 
     @staticmethod
-    def from_Euler_angles(phi,
-                          degrees = False):
+    def from_Euler_angles(phi: np.ndarray,
+                          degrees: bool = False) -> 'Rotation':
         """
         Initialize from Bunge Euler angles.
 
         Parameters
         ----------
-        phi : numpy.ndarray of shape (...,3)
+        phi : numpy.ndarray, shape (...,3)
             Euler angles (φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π])
             or (φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]) if degrees == True.
         degrees : bool, optional
@@ -725,16 +776,16 @@ class Rotation:
         return Rotation(Rotation._eu2qu(eu))
 
     @staticmethod
-    def from_axis_angle(axis_angle,
-                        degrees = False,
-                        normalize = False,
-                        P = -1):
+    def from_axis_angle(axis_angle: np.ndarray,
+                        degrees:bool = False,
+                        normalize: bool = False,
+                        P: Literal[1, -1] = -1) -> 'Rotation':
         """
         Initialize from Axis angle pair.
 
         Parameters
         ----------
-        axis_angle : numpy.ndarray of shape (...,4)
+        axis_angle : numpy.ndarray, shape (...,4)
             Axis and angle (n_1, n_2, n_3, ω) with ǀnǀ = 1 and ω ∈ [0,π]
             or ω ∈ [0,180] if degrees == True.
         degrees : bool, optional
@@ -763,15 +814,15 @@ class Rotation:
         return Rotation(Rotation._ax2qu(ax))
 
     @staticmethod
-    def from_basis(basis,
-                   orthonormal = True,
-                   reciprocal = False):
+    def from_basis(basis: np.ndarray,
+                   orthonormal: bool = True,
+                   reciprocal: bool = False) -> 'Rotation':
         """
         Initialize from lattice basis vectors.
 
         Parameters
         ----------
-        basis : numpy.ndarray of shape (...,3,3)
+        basis : numpy.ndarray, shape (...,3,3)
             Three three-dimensional lattice basis vectors.
         orthonormal : bool, optional
             Basis is strictly orthonormal, i.e. is free of stretch components. Defaults to True.
@@ -799,28 +850,29 @@ class Rotation:
         return Rotation(Rotation._om2qu(om))
 
     @staticmethod
-    def from_matrix(R):
+    def from_matrix(R: np.ndarray) -> 'Rotation':
         """
         Initialize from rotation matrix.
 
         Parameters
         ----------
-        R : numpy.ndarray of shape (...,3,3)
+        R : numpy.ndarray, shape (...,3,3)
             Rotation matrix with det(R) = 1, R.T ∙ R = I.
 
         """
         return Rotation.from_basis(R)
 
     @staticmethod
-    def from_parallel(a,b):
+    def from_parallel(a: np.ndarray,
+                      b: np.ndarray ) -> 'Rotation':
         """
         Initialize from pairs of two orthogonal lattice basis vectors.
 
         Parameters
         ----------
-        a : numpy.ndarray of shape (...,2,3)
+        a : numpy.ndarray, shape (...,2,3)
             Two three-dimensional lattice vectors of first orthogonal basis.
-        b : numpy.ndarray of shape (...,2,3)
+        b : numpy.ndarray, shape (...,2,3)
             Corresponding three-dimensional lattice vectors of second basis.
 
         """
@@ -840,15 +892,15 @@ class Rotation:
 
 
     @staticmethod
-    def from_Rodrigues_vector(rho,
-                              normalize = False,
-                              P = -1):
+    def from_Rodrigues_vector(rho: np.ndarray,
+                              normalize: bool = False,
+                              P: Literal[1, -1] = -1) -> 'Rotation':
         """
         Initialize from Rodrigues–Frank vector (angle separated from axis).
 
         Parameters
         ----------
-        rho : numpy.ndarray of shape (...,4)
+        rho : numpy.ndarray, shape (...,4)
             Rodrigues–Frank vector (n_1, n_2, n_3, tan(ω/2)) with ǀnǀ = 1  and ω ∈ [0,π].
         normalize : bool, optional
             Allow ǀnǀ ≠ 1. Defaults to False.
@@ -872,14 +924,14 @@ class Rotation:
         return Rotation(Rotation._ro2qu(ro))
 
     @staticmethod
-    def from_homochoric(h,
-                        P = -1):
+    def from_homochoric(h: np.ndarray,
+                        P: Literal[1, -1] = -1) -> 'Rotation':
         """
         Initialize from homochoric vector.
 
         Parameters
         ----------
-        h : numpy.ndarray of shape (...,3)
+        h : numpy.ndarray, shape (...,3)
             Homochoric vector (h_1, h_2, h_3) with ǀhǀ < (3/4*π)^(1/3).
         P : int ∈ {-1,1}, optional
             Sign convention. Defaults to -1.
@@ -899,14 +951,14 @@ class Rotation:
         return Rotation(Rotation._ho2qu(ho))
 
     @staticmethod
-    def from_cubochoric(x,
-                        P = -1):
+    def from_cubochoric(x: np.ndarray,
+                        P: Literal[1, -1] = -1) -> 'Rotation':
         """
         Initialize from cubochoric vector.
 
         Parameters
         ----------
-        x : numpy.ndarray of shape (...,3)
+        x : numpy.ndarray, shape (...,3)
             Cubochoric vector (x_1, x_2, x_3) with max(x_i) < 1/2*π^(2/3).
         P : int ∈ {-1,1}, optional
             Sign convention. Defaults to -1.
@@ -927,8 +979,8 @@ class Rotation:
 
 
     @staticmethod
-    def from_random(shape = None,
-                    rng_seed = None):
+    def from_random(shape: Tuple[int, ...] = None,
+                    rng_seed: NumpyRngSeed = None) -> 'Rotation':
         """
         Initialize with random rotation.
 
@@ -944,7 +996,7 @@ class Rotation:
 
         """
         rng = np.random.default_rng(rng_seed)
-        r = rng.random(3 if shape is None else tuple(shape)+(3,) if hasattr(shape, '__iter__') else (shape,3))
+        r = rng.random(3 if shape is None else tuple(shape)+(3,) if hasattr(shape, '__iter__') else (shape,3)) #type: ignore
 
         A = np.sqrt(r[...,2])
         B = np.sqrt(1.0-r[...,2])
@@ -957,20 +1009,20 @@ class Rotation:
 
 
     @staticmethod
-    def from_ODF(weights,
-                 phi,
-                 N = 500,
-                 degrees = True,
-                 fractions = True,
-                 rng_seed = None):
+    def from_ODF(weights: np.ndarray,
+                 phi: np.ndarray,
+                 N: int = 500,
+                 degrees: bool = True,
+                 fractions: bool = True,
+                 rng_seed: NumpyRngSeed = None) -> 'Rotation':
         """
         Sample discrete values from a binned orientation distribution function (ODF).
 
         Parameters
         ----------
-        weights : numpy.ndarray of shape (n)
+        weights : numpy.ndarray, shape (n)
             Texture intensity values (probability density or volume fraction) at Euler space grid points.
-        phi : numpy.ndarray of shape (n,3)
+        phi : numpy.ndarray, shape (n,3)
             Grid coordinates in Euler space at which weights are defined.
         N : integer, optional
             Number of discrete orientations to be sampled from the given ODF.
@@ -986,14 +1038,14 @@ class Rotation:
 
         Returns
         -------
-        samples : damask.Rotation of shape (N)
-            Array of sampled rotations closely representing the input ODF.
+        samples : damask.Rotation, shape (N)
+            Array of sampled rotations that approximate the input ODF.
 
         Notes
         -----
         Due to the distortion of Euler space in the vicinity of ϕ = 0, probability densities, p, defined on
         grid points with ϕ = 0 will never result in reconstructed orientations as their dV/V = p dγ = p × 0.
-        Hence, it is recommended to transform any such dataset to cell centers that avoid grid points at ϕ = 0.
+        Hence, it is recommended to transform any such dataset to a cell-centered version, which avoids grid points at ϕ = 0.
 
         References
         ----------
@@ -1015,11 +1067,11 @@ class Rotation:
 
 
     @staticmethod
-    def from_spherical_component(center,
-                                 sigma,
-                                 N = 500,
-                                 degrees = True,
-                                 rng_seed = None):
+    def from_spherical_component(center: 'Rotation',
+                                 sigma: float,
+                                 N: int = 500,
+                                 degrees: bool = True,
+                                 rng_seed: NumpyRngSeed = None) -> 'Rotation':
         """
         Calculate set of rotations with Gaussian distribution around center.
 
@@ -1050,20 +1102,20 @@ class Rotation:
 
 
     @staticmethod
-    def from_fiber_component(alpha,
-                             beta,
-                             sigma = 0.0,
-                             N = 500,
-                             degrees = True,
-                             rng_seed = None):
+    def from_fiber_component(alpha: IntSequence,
+                             beta: IntSequence,
+                             sigma: float = 0.0,
+                             N: int = 500,
+                             degrees: bool = True,
+                             rng_seed: NumpyRngSeed = None):
         """
         Calculate set of rotations with Gaussian distribution around direction.
 
         Parameters
         ----------
-        alpha : numpy.ndarray of shape (2)
+        alpha : numpy.ndarray, shape (2)
             Polar coordinates (phi from x, theta from z) of fiber direction in crystal frame.
-        beta : numpy.ndarray of shape (2)
+        beta : numpy.ndarray, shape (2)
             Polar coordinates (phi from x, theta from z) of fiber direction in sample frame.
         sigma : float, optional
             Standard deviation of (Gaussian) misorientation distribution.
@@ -1078,7 +1130,8 @@ class Rotation:
 
         """
         rng = np.random.default_rng(rng_seed)
-        sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta)
+        sigma_,alpha_,beta_ = (np.radians(coordinate) for coordinate in (sigma,alpha,beta)) if degrees else \
+                              map(np.array, (sigma,alpha,beta))
 
         d_cr  = np.array([np.sin(alpha_[0])*np.cos(alpha_[1]), np.sin(alpha_[0])*np.sin(alpha_[1]), np.cos(alpha_[0])])
         d_lab = np.array([np.sin( beta_[0])*np.cos( beta_[1]), np.sin( beta_[0])*np.sin( beta_[1]), np.cos( beta_[0])])
@@ -1132,7 +1185,7 @@ class Rotation:
 ####################################################################################################
     #---------- Quaternion ----------
     @staticmethod
-    def _qu2om(qu):
+    def _qu2om(qu: np.ndarray) -> np.ndarray:
         qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2)
         om = np.block([qq + 2.0*qu[...,1:2]**2,
                        2.0*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]),
@@ -1147,7 +1200,7 @@ class Rotation:
         return om
 
     @staticmethod
-    def _qu2eu(qu):
+    def _qu2eu(qu: np.ndarray) -> np.ndarray:
         """Quaternion to Bunge Euler angles."""
         q02   = qu[...,0:1]*qu[...,2:3]
         q13   = qu[...,1:2]*qu[...,3:4]
@@ -1175,7 +1228,7 @@ class Rotation:
         return eu
 
     @staticmethod
-    def _qu2ax(qu):
+    def _qu2ax(qu: np.ndarray) -> np.ndarray:
         """
         Quaternion to axis–angle pair.
 
@@ -1191,7 +1244,7 @@ class Rotation:
         return ax
 
     @staticmethod
-    def _qu2ro(qu):
+    def _qu2ro(qu: np.ndarray) -> np.ndarray:
         """Quaternion to Rodrigues–Frank vector."""
         with np.errstate(invalid='ignore',divide='ignore'):
             s  = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True)
@@ -1205,7 +1258,7 @@ class Rotation:
         return ro
 
     @staticmethod
-    def _qu2ho(qu):
+    def _qu2ho(qu: np.ndarray) -> np.ndarray:
         """Quaternion to homochoric vector."""
         with np.errstate(invalid='ignore'):
             omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
@@ -1216,14 +1269,14 @@ class Rotation:
         return ho
 
     @staticmethod
-    def _qu2cu(qu):
+    def _qu2cu(qu: np.ndarray) -> np.ndarray:
         """Quaternion to cubochoric vector."""
         return Rotation._ho2cu(Rotation._qu2ho(qu))
 
 
     #---------- Rotation matrix ----------
     @staticmethod
-    def _om2qu(om):
+    def _om2qu(om: np.ndarray) -> np.ndarray:
         """
         Rotation matrix to quaternion.
 
@@ -1265,7 +1318,7 @@ class Rotation:
         return qu
 
     @staticmethod
-    def _om2eu(om):
+    def _om2eu(om: np.ndarray) -> np.ndarray:
         """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)
@@ -1284,7 +1337,7 @@ class Rotation:
         return eu
 
     @staticmethod
-    def _om2ax(om):
+    def _om2ax(om: np.ndarray) -> np.ndarray:
         """Rotation matrix to axis–angle pair."""
         diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2],
                                    om[...,2,0:1]-om[...,0,2:3],
@@ -1305,24 +1358,24 @@ class Rotation:
         return ax
 
     @staticmethod
-    def _om2ro(om):
+    def _om2ro(om: np.ndarray) -> np.ndarray:
         """Rotation matrix to Rodrigues–Frank vector."""
         return Rotation._eu2ro(Rotation._om2eu(om))
 
     @staticmethod
-    def _om2ho(om):
+    def _om2ho(om: np.ndarray) -> np.ndarray:
         """Rotation matrix to homochoric vector."""
         return Rotation._ax2ho(Rotation._om2ax(om))
 
     @staticmethod
-    def _om2cu(om):
+    def _om2cu(om: np.ndarray) -> np.ndarray:
         """Rotation matrix to cubochoric vector."""
         return Rotation._ho2cu(Rotation._om2ho(om))
 
 
     #---------- Bunge Euler angles ----------
     @staticmethod
-    def _eu2qu(eu):
+    def _eu2qu(eu: np.ndarray) -> np.ndarray:
         """Bunge Euler angles to quaternion."""
         ee = 0.5*eu
         cPhi = np.cos(ee[...,1:2])
@@ -1335,7 +1388,7 @@ class Rotation:
         return qu
 
     @staticmethod
-    def _eu2om(eu):
+    def _eu2om(eu: np.ndarray) -> np.ndarray:
         """Bunge Euler angles to rotation matrix."""
         c = np.cos(eu)
         s = np.sin(eu)
@@ -1353,7 +1406,7 @@ class Rotation:
         return om
 
     @staticmethod
-    def _eu2ax(eu):
+    def _eu2ax(eu: np.ndarray) -> np.ndarray:
         """Bunge Euler angles to axis–angle pair."""
         t = np.tan(eu[...,1:2]*0.5)
         sigma = 0.5*(eu[...,0:1]+eu[...,2:3])
@@ -1372,7 +1425,7 @@ class Rotation:
         return ax
 
     @staticmethod
-    def _eu2ro(eu):
+    def _eu2ro(eu: np.ndarray) -> np.ndarray:
         """Bunge Euler angles to Rodrigues–Frank vector."""
         ax = Rotation._eu2ax(eu)
         ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)])
@@ -1381,19 +1434,19 @@ class Rotation:
         return ro
 
     @staticmethod
-    def _eu2ho(eu):
+    def _eu2ho(eu: np.ndarray) -> np.ndarray:
         """Bunge Euler angles to homochoric vector."""
         return Rotation._ax2ho(Rotation._eu2ax(eu))
 
     @staticmethod
-    def _eu2cu(eu):
+    def _eu2cu(eu: np.ndarray) -> np.ndarray:
         """Bunge Euler angles to cubochoric vector."""
         return Rotation._ho2cu(Rotation._eu2ho(eu))
 
 
     #---------- Axis angle pair ----------
     @staticmethod
-    def _ax2qu(ax):
+    def _ax2qu(ax: np.ndarray) -> np.ndarray:
         """Axis–angle pair to quaternion."""
         c = np.cos(ax[...,3:4]*.5)
         s = np.sin(ax[...,3:4]*.5)
@@ -1401,7 +1454,7 @@ class Rotation:
         return qu
 
     @staticmethod
-    def _ax2om(ax):
+    def _ax2om(ax: np.ndarray) -> np.ndarray:
         """Axis-angle pair to rotation matrix."""
         c = np.cos(ax[...,3:4])
         s = np.sin(ax[...,3:4])
@@ -1418,12 +1471,12 @@ class Rotation:
         return om if _P < 0.0 else np.swapaxes(om,-1,-2)
 
     @staticmethod
-    def _ax2eu(ax):
+    def _ax2eu(ax: np.ndarray) -> np.ndarray:
         """Rotation matrix to Bunge Euler angles."""
         return Rotation._om2eu(Rotation._ax2om(ax))
 
     @staticmethod
-    def _ax2ro(ax):
+    def _ax2ro(ax: np.ndarray) -> np.ndarray:
         """Axis–angle pair to Rodrigues–Frank vector."""
         ro = np.block([ax[...,:3],
                        np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0),
@@ -1434,36 +1487,36 @@ class Rotation:
         return ro
 
     @staticmethod
-    def _ax2ho(ax):
+    def _ax2ho(ax: np.ndarray) -> np.ndarray:
         """Axis–angle pair to homochoric vector."""
         f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0)
         ho = ax[...,:3] * f
         return ho
 
     @staticmethod
-    def _ax2cu(ax):
+    def _ax2cu(ax: np.ndarray) -> np.ndarray:
         """Axis–angle pair to cubochoric vector."""
         return Rotation._ho2cu(Rotation._ax2ho(ax))
 
 
     #---------- Rodrigues-Frank vector ----------
     @staticmethod
-    def _ro2qu(ro):
+    def _ro2qu(ro: np.ndarray) -> np.ndarray:
         """Rodrigues–Frank vector to quaternion."""
         return Rotation._ax2qu(Rotation._ro2ax(ro))
 
     @staticmethod
-    def _ro2om(ro):
+    def _ro2om(ro: np.ndarray) -> np.ndarray:
         """Rodgrigues–Frank vector to rotation matrix."""
         return Rotation._ax2om(Rotation._ro2ax(ro))
 
     @staticmethod
-    def _ro2eu(ro):
+    def _ro2eu(ro: np.ndarray) -> np.ndarray:
         """Rodrigues–Frank vector to Bunge Euler angles."""
         return Rotation._om2eu(Rotation._ro2om(ro))
 
     @staticmethod
-    def _ro2ax(ro):
+    def _ro2ax(ro: np.ndarray) -> np.ndarray:
         """Rodrigues–Frank vector to axis–angle pair."""
         with np.errstate(invalid='ignore',divide='ignore'):
             ax = np.where(np.isfinite(ro[...,3:4]),
@@ -1473,7 +1526,7 @@ class Rotation:
         return ax
 
     @staticmethod
-    def _ro2ho(ro):
+    def _ro2ho(ro: np.ndarray) -> np.ndarray:
         """Rodrigues–Frank vector to homochoric vector."""
         f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi)
         ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-8,ro[...,0:3].shape),
@@ -1481,29 +1534,29 @@ class Rotation:
         return ho
 
     @staticmethod
-    def _ro2cu(ro):
+    def _ro2cu(ro: np.ndarray) -> np.ndarray:
         """Rodrigues–Frank vector to cubochoric vector."""
         return Rotation._ho2cu(Rotation._ro2ho(ro))
 
 
     #---------- Homochoric vector----------
     @staticmethod
-    def _ho2qu(ho):
+    def _ho2qu(ho: np.ndarray) -> np.ndarray:
         """Homochoric vector to quaternion."""
         return Rotation._ax2qu(Rotation._ho2ax(ho))
 
     @staticmethod
-    def _ho2om(ho):
+    def _ho2om(ho: np.ndarray) -> np.ndarray:
         """Homochoric vector to rotation matrix."""
         return Rotation._ax2om(Rotation._ho2ax(ho))
 
     @staticmethod
-    def _ho2eu(ho):
+    def _ho2eu(ho: np.ndarray) -> np.ndarray:
         """Homochoric vector to Bunge Euler angles."""
         return Rotation._ax2eu(Rotation._ho2ax(ho))
 
     @staticmethod
-    def _ho2ax(ho):
+    def _ho2ax(ho: np.ndarray) -> np.ndarray:
         """Homochoric vector to axis–angle pair."""
         tfit = np.array([+1.0000000000018852,      -0.5000000002194847,
                          -0.024999992127593126,    -0.003928701544781374,
@@ -1526,12 +1579,12 @@ class Rotation:
         return ax
 
     @staticmethod
-    def _ho2ro(ho):
+    def _ho2ro(ho: np.ndarray) -> np.ndarray:
         """Axis–angle pair to Rodrigues–Frank vector."""
         return Rotation._ax2ro(Rotation._ho2ax(ho))
 
     @staticmethod
-    def _ho2cu(ho):
+    def _ho2cu(ho: np.ndarray) -> np.ndarray:
         """
         Homochoric vector to cubochoric vector.
 
@@ -1571,32 +1624,32 @@ class Rotation:
 
     #---------- Cubochoric ----------
     @staticmethod
-    def _cu2qu(cu):
+    def _cu2qu(cu: np.ndarray) -> np.ndarray:
         """Cubochoric vector to quaternion."""
         return Rotation._ho2qu(Rotation._cu2ho(cu))
 
     @staticmethod
-    def _cu2om(cu):
+    def _cu2om(cu: np.ndarray) -> np.ndarray:
         """Cubochoric vector to rotation matrix."""
         return Rotation._ho2om(Rotation._cu2ho(cu))
 
     @staticmethod
-    def _cu2eu(cu):
+    def _cu2eu(cu: np.ndarray) -> np.ndarray:
         """Cubochoric vector to Bunge Euler angles."""
         return Rotation._ho2eu(Rotation._cu2ho(cu))
 
     @staticmethod
-    def _cu2ax(cu):
+    def _cu2ax(cu: np.ndarray) -> np.ndarray:
         """Cubochoric vector to axis–angle pair."""
         return Rotation._ho2ax(Rotation._cu2ho(cu))
 
     @staticmethod
-    def _cu2ro(cu):
+    def _cu2ro(cu: np.ndarray) -> np.ndarray:
         """Cubochoric vector to Rodrigues–Frank vector."""
         return Rotation._ho2ro(Rotation._cu2ho(cu))
 
     @staticmethod
-    def _cu2ho(cu):
+    def _cu2ho(cu: np.ndarray) -> np.ndarray:
         """
         Cubochoric vector to homochoric vector.
 
@@ -1639,7 +1692,8 @@ class Rotation:
 
 
     @staticmethod
-    def _get_pyramid_order(xyz,direction=None):
+    def _get_pyramid_order(xyz: np.ndarray,
+                           direction: Literal['forward', 'backward']) -> np.ndarray:
         """
         Get order of the coordinates.
 
diff --git a/python/damask/_table.py b/python/damask/_table.py
index df364d46e..c7b367248 100644
--- a/python/damask/_table.py
+++ b/python/damask/_table.py
@@ -194,7 +194,7 @@ class Table:
 
         Returns
         -------
-        mask : numpy.ndarray bool
+        mask : numpy.ndarray of bool
             Mask indicating where corresponding table values are close.
 
         """
diff --git a/python/damask/_typehints.py b/python/damask/_typehints.py
index 0b4a56a69..5fcf39a41 100644
--- a/python/damask/_typehints.py
+++ b/python/damask/_typehints.py
@@ -1,6 +1,6 @@
 """Functionality for typehints."""
 
-from typing import Sequence, Union, TextIO
+from typing import Sequence, Union, Literal, TextIO
 from pathlib import Path
 
 import numpy as np
@@ -9,6 +9,9 @@ import numpy as np
 FloatSequence = Union[np.ndarray,Sequence[float]]
 IntSequence = Union[np.ndarray,Sequence[int]]
 FileHandle = Union[TextIO, str, Path]
+CrystalFamily = Union[None,Literal['triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic']]
+CrystalLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
+CrystalKinematics = Literal['slip', 'twin']
 NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator]
 # BitGenerator does not exists in older numpy versions
 #NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator]
diff --git a/python/damask/util.py b/python/damask/util.py
index fc3f8ea8f..ab34df7c0 100644
--- a/python/damask/util.py
+++ b/python/damask/util.py
@@ -9,7 +9,7 @@ import re
 import fractions
 from collections import abc
 from functools import reduce
-from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, SupportsIndex, Sequence
+from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal
 from pathlib import Path
 
 import numpy as np
@@ -427,7 +427,7 @@ def hybrid_IA(dist: np.ndarray,
 def shapeshifter(fro: Tuple[int, ...],
                  to: Tuple[int, ...],
                  mode: Literal['left','right'] = 'left',
-                 keep_ones: bool = False) -> Sequence[SupportsIndex]:
+                 keep_ones: bool = False) -> Tuple[int, ...]:
     """
     Return dimensions that reshape 'fro' to become broadcastable to 'to'.
 
@@ -490,7 +490,7 @@ def shapeshifter(fro: Tuple[int, ...],
 
 
 def shapeblender(a: Tuple[int, ...],
-                 b: Tuple[int, ...]) -> Sequence[SupportsIndex]:
+                 b: Tuple[int, ...]) -> Tuple[int, ...]:
     """
     Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.