diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index e727c54ae..6c0dd2fa3 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -1,5 +1,6 @@
 import inspect
 import copy
+from typing import Union, Callable, Sequence, Dict, Any, Tuple, List
 
 import numpy as np
 
@@ -93,12 +94,12 @@ 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[Sequence[float], np.ndarray, Rotation] = np.array([1.,0.,0.,0.]), *,
+                 family: str = None,
+                 lattice: str = 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 +116,12 @@ 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,rotation: Union[Sequence[float], np.ndarray, Rotation] = None) -> "Orientation":
         """Create deep copy."""
         dup = copy.deepcopy(self)
         if rotation is not None:
@@ -131,7 +131,8 @@ class Orientation(Rotation,Crystal):
     copy = __copy__
 
 
-    def __eq__(self,other):
+
+    def __eq__(self, other: object) -> bool:
         """
         Equal to other.
 
@@ -141,12 +142,14 @@ class Orientation(Rotation,Crystal):
             Orientation to check for equality.
 
         """
+        if not isinstance(other, Orientation):
+            raise TypeError
         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 +159,16 @@ class Orientation(Rotation,Crystal):
             Orientation to check for equality.
 
         """
+        if not isinstance(other, Orientation):
+            raise TypeError
         return np.logical_not(self==other)
 
 
-    def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def isclose(self,
+                other: object,
+                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.
 
@@ -180,6 +189,8 @@ class Orientation(Rotation,Crystal):
             Mask indicating where corresponding orientations are close.
 
         """
+        if not isinstance(other, Orientation):
+            raise TypeError
         matching_type = self.family == other.family and \
                         self.lattice == other.lattice and \
                         self.parameters == other.parameters
@@ -187,7 +198,11 @@ class Orientation(Rotation,Crystal):
 
 
 
-    def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def allclose(self,
+                 other: object,
+                 rtol: float = 1e-5,
+                 atol: float = 1e-8,
+                 equal_nan: bool = True) -> Union[bool, np.bool_]:
         """
         Test whether all values are approximately equal to corresponding ones of other Orientation.
 
@@ -208,10 +223,12 @@ class Orientation(Rotation,Crystal):
             Whether all values are close between both orientations.
 
         """
+        if not isinstance(other, Orientation):
+            raise TypeError
         return np.all(self.isclose(other,rtol,atol,equal_nan))
 
 
-    def __mul__(self,other):
+    def __mul__(self, other: Union[Rotation, "Orientation"]) -> "Orientation":
         """
         Compose this orientation with other.
 
@@ -233,7 +250,7 @@ class Orientation(Rotation,Crystal):
 
 
     @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 +269,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,85 +281,88 @@ 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: Union[Sequence[float], np.ndarray],
+                        hkl: Union[Sequence[float], np.ndarray],
+                        **kwargs) -> "Orientation":
         """
         Initialize orientation object from two crystallographic directions.
 
@@ -362,7 +382,7 @@ class Orientation(Rotation,Crystal):
 
 
     @property
-    def equivalent(self):
+    def equivalent(self) -> "Orientation":
         """
         Orientations that are symmetrically equivalent.
 
@@ -376,7 +396,7 @@ class Orientation(Rotation,Crystal):
 
 
     @property
-    def reduced(self):
+    def reduced(self) -> "Orientation":
         """Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry."""
         eq   = self.equivalent
         ok   = eq.in_FZ
@@ -385,9 +405,8 @@ class Orientation(Rotation,Crystal):
         sort = 0 if len(loc) == 1 else np.lexsort(loc[:0:-1])
         return eq[ok][sort].reshape(self.shape)
 
-
     @property
-    def in_FZ(self):
+    def in_FZ(self) -> Union[np.bool_, np.ndarray]:
         """
         Check whether orientation falls into fundamental zone of own symmetry.
 
@@ -431,7 +450,7 @@ 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.
 
@@ -471,8 +490,7 @@ class Orientation(Rotation,Crystal):
             else:
                 return np.ones_like(rho[...,0],dtype=bool)
 
-
-    def disorientation(self,other,return_operators=False):
+    def disorientation(self, other, return_operators = False):
         """
         Calculate disorientation between myself and given other orientation.
 
@@ -518,8 +536,8 @@ class Orientation(Rotation,Crystal):
         if self.family != other.family:
             raise NotImplementedError('disorientation between different crystal families')
 
-        blend = util.shapeblender(self.shape,other.shape)
-        s =  self.equivalent
+        blend: Tuple[int, ...] = util.shapeblender(self.shape,other.shape)
+        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')
@@ -534,10 +552,9 @@ class Orientation(Rotation,Crystal):
         r = np.where(np.any(forward[...,np.newaxis],axis=(0,1),keepdims=True),
                      r_.quaternion,
                      _r.quaternion)
-        loc  = np.where(ok)
-        sort = 0 if len(loc) == 2 else np.lexsort(loc[:1:-1])
-        quat = r[ok][sort].reshape(blend+(4,))
-
+        loc: Tuple[float]  = np.where(ok)
+        sort: np.ndarray = 0 if len(loc) == 2 else np.lexsort(loc[:1:-1])
+        quat: np.ndarray = r[ok][sort].reshape(blend+(4,))
         return (
                 (self.copy(rotation=quat),
                  (np.vstack(loc[:2]).T)[sort].reshape(blend+(2,)))
@@ -546,7 +563,7 @@ class Orientation(Rotation,Crystal):
                )
 
 
-    def average(self,weights=None,return_cloud=False):
+    def average(self, weights = None, return_cloud = False):
         """
         Return orientation average over last dimension.
 
@@ -587,7 +604,10 @@ class Orientation(Rotation,Crystal):
                )
 
 
-    def to_SST(self,vector,proper=False,return_operators=False):
+    def to_SST(self,
+               vector: np.ndarray,
+               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.
 
@@ -626,7 +646,7 @@ class Orientation(Rotation,Crystal):
                )
 
 
-    def in_SST(self,vector,proper=False):
+    def in_SST(self, vector: np.ndarray, proper: bool = False) -> Union[np.bool_, np.ndarray]:
         """
         Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry.
 
@@ -667,7 +687,7 @@ class Orientation(Rotation,Crystal):
             return np.all(components >= 0.0,axis=-1)
 
 
-    def IPF_color(self,vector,in_SST=True,proper=False):
+    def IPF_color(self, vector: np.ndarray, in_SST: bool = True, proper: bool = False) -> np.ndarray:
         """
         Map vector to RGB color within standard stereographic triangle of own symmetry.
 
@@ -715,30 +735,30 @@ class Orientation(Rotation,Crystal):
             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) \
+            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 = 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[str, List[List]] = {
             'cubic':         [
                               [ 1.0,            0.0,            0.0,            0.0            ],
                               [ 0.0,            1.0,            0.0,            0.0            ],
@@ -808,7 +828,10 @@ 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: np.ndarray = None,
+                hkl: np.ndarray = None,
+                with_symmetry: bool = False) -> np.ndarray:
         """
         Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl).
 
@@ -839,7 +862,9 @@ class Orientation(Rotation,Crystal):
                @ np.broadcast_to(v,blend+(3,))
 
 
-    def Schmid(self,*,N_slip=None,N_twin=None):
+    def Schmid(self, *,
+               N_slip: Sequence[int] = None,
+               N_twin: Sequence[int] = None) -> np.ndarray:
         u"""
         Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems.
 
@@ -876,6 +901,7 @@ class Orientation(Rotation,Crystal):
                             (self.kinematics('twin'),N_twin)
         if active == '*': active = [len(a) for a in kinematics['direction']]
 
+        assert active
         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),
@@ -886,7 +912,7 @@ class Orientation(Rotation,Crystal):
                @ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape)
 
 
-    def related(self,model):
+    def related(self, model: str) -> "Orientation":
         """
         Orientations derived from the given relationship.
 
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index ac921d70a..0a4421090 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -6,6 +6,8 @@ from . import tensor
 from . import util
 from . import grid_filters
 
+from typing import Union, Sequence, Tuple, Literal, List
+
 _P = -1
 
 # parameters for conversion from/to cubochoric
@@ -61,7 +63,7 @@ class Rotation:
 
     __slots__ = ['quaternion']
 
-    def __init__(self,rotation = np.array([1.0,0.0,0.0,0.0])):
+    def __init__(self, rotation: Union[Sequence[float], np.ndarray, "Rotation"] = np.array([1.0,0.0,0.0,0.0])):
         """
         New rotation.
 
@@ -73,6 +75,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 +84,13 @@ 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, rotation: Union[Sequence[float], np.ndarray, "Rotation"] = None) -> "Rotation":
         """Create deep copy."""
         dup = copy.deepcopy(self)
         if rotation is not None:
@@ -97,14 +100,14 @@ 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])
 
 
-    def __eq__(self,other):
+    def __eq__(self, other: object) -> bool:
         """
         Equal to other.
 
@@ -114,11 +117,13 @@ class Rotation:
             Rotation to check for equality.
 
         """
+        if not isinstance(other, Rotation):
+            raise TypeError
         return 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 +133,12 @@ class Rotation:
             Rotation to check for inequality.
 
         """
+        if not isinstance(other, Rotation):
+            raise TypeError
         return np.logical_not(self==other)
 
 
-    def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True):
+    def isclose(self, other: "Rotation", 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.
 
@@ -158,7 +165,7 @@ 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, other: "Rotation", 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 +195,27 @@ 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) -> "Rotation":
         """Inverse rotation (backward rotation)."""
         dup = self.copy()
         dup.quaternion[...,1:] *= -1
         return dup
 
 
-    def __pow__(self,exp):
+    def __pow__(self, exp: int) -> "Rotation":
         """
         Perform the rotation 'exp' times.
 
@@ -222,7 +229,7 @@ class Rotation:
         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())
 
-    def __ipow__(self,exp):
+    def __ipow__(self, exp: int) -> "Rotation":
         """
         Perform the rotation 'exp' times (in-place).
 
@@ -235,7 +242,7 @@ class Rotation:
         return self**exp
 
 
-    def __mul__(self,other):
+    def __mul__(self, other: "Rotation") -> "Rotation":
         """
         Compose with other.
 
@@ -261,7 +268,7 @@ class Rotation:
         else:
             raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
 
-    def __imul__(self,other):
+    def __imul__(self, other: "Rotation") -> "Rotation":
         """
         Compose with other (in-place).
 
@@ -274,7 +281,7 @@ class Rotation:
         return self*other
 
 
-    def __truediv__(self,other):
+    def __truediv__(self, other: "Rotation") -> "Rotation":
         """
         Compose with inverse of other.
 
@@ -294,7 +301,7 @@ 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, other: "Rotation") -> "Rotation":
         """
         Compose with inverse of other (in-place).
 
@@ -307,7 +314,7 @@ class Rotation:
         return self/other
 
 
-    def __matmul__(self,other):
+    def __matmul__(self, other: np.ndarray) -> np.ndarray:
         """
         Rotate vector, second order tensor, or fourth order tensor.
 
@@ -350,13 +357,13 @@ class Rotation:
     apply = __matmul__
 
 
-    def _standardize(self):
+    def _standardize(self) -> "Rotation":
         """Standardize quaternion (ensure positive real hemisphere)."""
         self.quaternion[self.quaternion[...,0] < 0.0] *= -1
         return self
 
 
-    def append(self,other):
+    def append(self, other: Union["Rotation", List["Rotation"]]) -> "Rotation":
         """
         Extend array along first dimension with other array(s).
 
@@ -369,7 +376,7 @@ class Rotation:
                                                       [self]+other if isinstance(other,list) else [self,other]))))
 
 
-    def flatten(self,order = 'C'):
+    def flatten(self, order: Literal['C','F','A'] = 'C') -> "Rotation":
         """
         Flatten array.
 
@@ -382,7 +389,7 @@ class Rotation:
         return self.copy(rotation=self.quaternion.reshape((-1,4),order=order))
 
 
-    def reshape(self,shape,order = 'C'):
+    def reshape(self,shape: Union[int, Tuple[int, ...]], order: Literal['C','F','A'] = 'C') -> "Rotation":
         """
         Reshape array.
 
@@ -396,7 +403,7 @@ class Rotation:
         return self.copy(rotation=self.quaternion.reshape(tuple(shape)+(4,),order=order))
 
 
-    def broadcast_to(self,shape,mode = 'right'):
+    def broadcast_to(self, shape: Union[int, Tuple[int, ...]], mode: Literal["left", "right"] = 'right') -> "Rotation":
         """
         Broadcast array.
 
@@ -458,7 +465,7 @@ class Rotation:
                                         accept_homomorph = True)
 
 
-    def misorientation(self,other):
+    def misorientation(self, other: "Rotation") -> "Rotation":
         """
         Calculate misorientation to other Rotation.
 
@@ -479,7 +486,7 @@ class Rotation:
     ################################################################################################
     # convert to different orientation representations (numpy arrays)
 
-    def as_quaternion(self):
+    def as_quaternion(self) -> np.ndarray:
         """
         Represent as unit quaternion.
 
@@ -492,7 +499,7 @@ class Rotation:
         return self.quaternion.copy()
 
     def as_Euler_angles(self,
-                        degrees = False):
+                        degrees: bool = False) -> np.ndarray:
         """
         Represent as Bunge Euler angles.
 
@@ -526,8 +533,8 @@ class Rotation:
         return eu
 
     def as_axis_angle(self,
-                      degrees = False,
-                      pair = False):
+                      degrees: bool = False,
+                      pair: bool = False) -> Union[Tuple[float, ...], np.ndarray]:
         """
         Represent as axis–angle pair.
 
@@ -554,11 +561,11 @@ 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.
 
@@ -582,7 +589,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.
 
@@ -615,7 +622,7 @@ class Rotation:
         else:
             return ro
 
-    def as_homochoric(self):
+    def as_homochoric(self) -> np.ndarray:
         """
         Represent as homochoric vector.
 
@@ -636,7 +643,7 @@ class Rotation:
         """
         return Rotation._qu2ho(self.quaternion)
 
-    def as_cubochoric(self):
+    def as_cubochoric(self) -> np.ndarray:
         """
         Represent as cubochoric vector.
 
@@ -661,9 +668,9 @@ 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[Sequence[float]], np.ndarray],
+                        accept_homomorph: bool = False,
+                        P: Literal[1, -1] = -1) -> "Rotation":
         """
         Initialize from quaternion.
 
@@ -678,7 +685,7 @@ class Rotation:
             Sign convention. Defaults to -1.
 
         """
-        qu = np.array(q,dtype=float)
+        qu: np.ndarray = np.array(q,dtype=float)
         if qu.shape[:-2:-1] != (4,):
             raise ValueError('Invalid shape.')
         if abs(P) != 1:
@@ -696,8 +703,8 @@ 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.
 
@@ -725,10 +732,10 @@ 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.
 
@@ -763,9 +770,9 @@ 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.
 
@@ -799,7 +806,7 @@ class Rotation:
         return Rotation(Rotation._om2qu(om))
 
     @staticmethod
-    def from_matrix(R):
+    def from_matrix(R: np.ndarray) -> "Rotation":
         """
         Initialize from rotation matrix.
 
@@ -812,8 +819,9 @@ class Rotation:
         return Rotation.from_basis(R)
 
     @staticmethod
-    def from_parallel(a,b,
-                      **kwargs):
+    def from_parallel(a: np.ndarray,
+                      b: np.ndarray,
+                      **kwargs) -> "Rotation":
         """
         Initialize from pairs of two orthogonal lattice basis vectors.
 
@@ -841,9 +849,9 @@ 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).
 
@@ -873,8 +881,8 @@ 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.
 
@@ -900,8 +908,8 @@ 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.
 
@@ -928,8 +936,8 @@ class Rotation:
 
 
     @staticmethod
-    def from_random(shape = None,
-                    rng_seed = None):
+    def from_random(shape: Tuple[int, ...] = None,
+                    rng_seed: Union[None, int, Sequence[int], np.ndarray] = None) -> "Rotation":
         """
         Initialize with random rotation.
 
@@ -958,13 +966,13 @@ class Rotation:
 
 
     @staticmethod
-    def from_ODF(weights,
-                 phi,
-                 N = 500,
-                 degrees = True,
-                 fractions = True,
-                 rng_seed = None,
-                 **kwargs):
+    def from_ODF(weights: np.ndarray,
+                 phi: np.ndarray,
+                 N: int = 500,
+                 degrees: bool = True,
+                 fractions: bool = True,
+                 rng_seed: Union[None, int, Sequence[int], np.ndarray] = None,
+                 **kwargs) -> "Rotation":
         """
         Sample discrete values from a binned orientation distribution function (ODF).
 
@@ -1017,11 +1025,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: Union[None, int, Sequence[float], np.ndarray] = None) -> "Rotation":
         """
         Calculate set of rotations with Gaussian distribution around center.
 
@@ -1052,12 +1060,12 @@ class Rotation:
 
 
     @staticmethod
-    def from_fiber_component(alpha,
-                             beta,
-                             sigma = 0.0,
-                             N = 500,
-                             degrees = True,
-                             rng_seed = None):
+    def from_fiber_component(alpha: np.ndarray,
+                             beta: np.ndarray,
+                             sigma: float = 0.0,
+                             N: int = 500,
+                             degrees: bool = True,
+                             rng_seed: bool = None):
         """
         Calculate set of rotations with Gaussian distribution around direction.
 
@@ -1080,7 +1088,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_: np.ndarray; alpha_: np.ndarray; beta_: np.ndarray
+        sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta) # type: ignore
 
         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])])
@@ -1134,7 +1143,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]),
@@ -1149,7 +1158,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]
@@ -1177,7 +1186,7 @@ class Rotation:
         return eu
 
     @staticmethod
-    def _qu2ax(qu):
+    def _qu2ax(qu: np.ndarray) -> np.ndarray:
         """
         Quaternion to axis–angle pair.
 
@@ -1193,7 +1202,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)
@@ -1207,7 +1216,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))
@@ -1218,14 +1227,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.
 
@@ -1267,7 +1276,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)
@@ -1286,7 +1295,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],
@@ -1307,24 +1316,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])
@@ -1337,7 +1346,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)
@@ -1355,7 +1364,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])
@@ -1374,7 +1383,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)])
@@ -1383,19 +1392,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)
@@ -1403,7 +1412,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])
@@ -1420,12 +1429,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),
@@ -1436,36 +1445,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]),
@@ -1475,7 +1484,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),
@@ -1483,29 +1492,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,
@@ -1528,12 +1537,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.
 
@@ -1573,32 +1582,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.
 
@@ -1641,7 +1650,7 @@ 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.