diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index 007aecfd9..349c368fa 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -3,214 +3,8 @@
 import math
 import numpy as np
 from . import Lambert
-
-P = -1
-
-####################################################################################################
-class Quaternion:
-    u"""
-    Quaternion with basic operations
-    
-    q is the real part, p = (x, y, z) are the imaginary parts.
-    Defintion of multiplication depends on variable P, P ∉ {-1,1}.
-    """
-
-    def __init__(self,
-                 q    = 0.0,
-                 p    = np.zeros(3,dtype=float)):
-      """Initializes to identity unless specified"""
-      self.q = q
-      self.p = np.array(p)
-
-
-    def __copy__(self):
-      """Copy"""
-      return self.__class__(q=self.q,
-                            p=self.p.copy())
-
-    copy = __copy__
-
-
-    def __iter__(self):
-      """Components"""
-      return iter(self.asList())
-
-    def __repr__(self):
-      """Readable string"""
-      return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
-      
-
-    def __add__(self, other):
-      """Addition"""
-      if isinstance(other, Quaternion):
-        return self.__class__(q=self.q + other.q,
-                              p=self.p + other.p)
-      else:
-        return NotImplemented
-
-    def __iadd__(self, other):
-      """In-place addition"""
-      if isinstance(other, Quaternion):
-        self.q += other.q
-        self.p += other.p
-        return self
-      else:
-        return NotImplemented
-      
-    def __pos__(self):                                                                      
-      """Unary positive operator"""                                                                       
-      return self
-      
-      
-    def __sub__(self, other):
-      """Subtraction"""
-      if isinstance(other, Quaternion):
-        return self.__class__(q=self.q - other.q,
-                              p=self.p - other.p)
-      else:
-        return NotImplemented
-
-    def __isub__(self, other):
-      """In-place subtraction"""
-      if isinstance(other, Quaternion):
-        self.q -= other.q
-        self.p -= other.p
-        return self
-      else:
-        return NotImplemented
-        
-    def __neg__(self):                                                                      
-      """Unary positive operator"""  
-      self.q *= -1.0
-      self.p *= -1.0                                                                     
-      return self
-      
-      
-    def __mul__(self, other):
-      """Multiplication with quaternion or scalar"""
-      if isinstance(other, Quaternion):
-        return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
-                              p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
-      elif isinstance(other, (int, float)):
-        return self.__class__(q=self.q*other,
-                              p=self.p*other)
-      else:
-        return NotImplemented
-
-    def __imul__(self, other):
-      """In-place multiplication with quaternion or scalar"""
-      if isinstance(other, Quaternion):
-        self.q = self.q*other.q - np.dot(self.p,other.p)
-        self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
-        return self
-      elif isinstance(other, (int, float)):
-        self *= other
-        return self
-      else:
-        return NotImplemented
-
-
-    def __truediv__(self, other):
-      """Divsion with quaternion or scalar"""
-      if isinstance(other, Quaternion):
-        s = other.conjugate()/abs(other)**2.
-        return self.__class__(q=self.q * s,
-                              p=self.p * s)
-      elif isinstance(other, (int, float)):
-        self.q /= other
-        self.p /= other
-        return self
-      else:
-        return NotImplemented
-
-    def __itruediv__(self, other):
-      """In-place divsion with quaternion or scalar"""
-      if isinstance(other, Quaternion):
-        s = other.conjugate()/abs(other)**2.
-        self *= s
-        return self
-      elif isinstance(other, (int, float)):
-        self.q /= other
-        return self
-      else:
-        return NotImplemented
-
-
-    def __pow__(self, exponent):
-      """Power"""
-      if isinstance(exponent, (int, float)):
-        omega = np.acos(self.q)
-        return self.__class__(q=         np.cos(exponent*omega),
-                              p=self.p * np.sin(exponent*omega)/np.sin(omega))
-      else:
-        return NotImplemented
-
-    def __ipow__(self, exponent):
-      """In-place power"""
-      if isinstance(exponent, (int, float)):
-        omega = np.acos(self.q)
-        self.q  = np.cos(exponent*omega)
-        self.p *= np.sin(exponent*omega)/np.sin(omega)
-      else:
-        return NotImplemented
-      
-
-    def __abs__(self):
-      """Norm"""
-      return math.sqrt(self.q ** 2 + np.dot(self.p,self.p))
-
-    magnitude = __abs__
-
-
-    def __eq__(self,other):
-      """Equal (sufficiently close) to each other"""
-      return np.isclose(( self-other).magnitude(),0.0) \
-          or np.isclose((-self-other).magnitude(),0.0)
-
-    def __ne__(self,other):
-      """Not equal (sufficiently close) to each other"""
-      return not self.__eq__(other)
-
-
-    def asM(self):
-      """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)"""
-      return np.outer(self.asArray(),self.asArray())
-
-    def asArray(self):
-      """As numpy array"""
-      return np.array((self.q,self.p[0],self.p[1],self.p[2]))
- 
-    def asList(self):
-      return [self.q]+list(self.p)
-
-    def normalize(self):
-      d = self.magnitude()
-      if d > 0.0:
-          self.q /= d
-          self.p /= d
-      return self
-      
-    def normalized(self):
-      return self.copy().normalize()
-
-  
-    def conjugate(self):
-      self.p = -self.p
-      return self
-      
-    def conjugated(self):
-      return self.copy().conjugate()
-
-
-    def homomorph(self):
-      if self.q < 0.0:
-        self.q = -self.q
-        self.p = -self.p
-      return self
-      
-    def homomorphed(self):
-      return self.copy().homomorph()
-      
+from quaternion import Quaternion
+from quaternion import P as P
 
 
 ####################################################################################################
@@ -488,7 +282,7 @@ class Rotation:
 
       for i,(r,n) in enumerate(zip(rotations,weights)):
         M =          r.asM() * n if i == 0 \
-            else M + r.asM() * n                                                                      # noqa add (multiples) of this rotation to average noqa
+            else M + r.asM() * n                                                                    # noqa add (multiples) of this rotation to average noqa
       eig, vec = np.linalg.eig(M/N)
 
       return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py
new file mode 100644
index 000000000..7c98a6d1c
--- /dev/null
+++ b/python/damask/quaternion.py
@@ -0,0 +1,210 @@
+# -*- coding: UTF-8 no BOM -*-
+
+import numpy as np
+
+P = -1                                                                                              # convention (sed DOI:10.1088/0965-0393/23/8/083501)
+
+####################################################################################################
+class Quaternion:
+    u"""
+    Quaternion with basic operations
+    
+    q is the real part, p = (x, y, z) are the imaginary parts.
+    Defintion of multiplication depends on variable P, P ∈ {-1,1}.
+    """                                                                              
+
+    def __init__(self,
+                 q    = 0.0,
+                 p    = np.zeros(3,dtype=float)):
+      """Initializes to identity unless specified"""
+      self.q = q
+      self.p = np.array(p)
+
+
+    def __copy__(self):
+      """Copy"""
+      return self.__class__(q=self.q,
+                            p=self.p.copy())
+
+    copy = __copy__
+
+
+    def __iter__(self):
+      """Components"""
+      return iter(self.asList())
+
+    def __repr__(self):
+      """Readable string"""
+      return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
+      
+
+    def __add__(self, other):
+      """Addition"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q + other.q,
+                              p=self.p + other.p)
+      else:
+        return NotImplemented
+
+    def __iadd__(self, other):
+      """In-place addition"""
+      if isinstance(other, Quaternion):
+        self.q += other.q
+        self.p += other.p
+        return self
+      else:
+        return NotImplemented
+      
+    def __pos__(self):                                                                      
+      """Unary positive operator"""                                                                       
+      return self
+      
+      
+    def __sub__(self, other):
+      """Subtraction"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q - other.q,
+                              p=self.p - other.p)
+      else:
+        return NotImplemented
+
+    def __isub__(self, other):
+      """In-place subtraction"""
+      if isinstance(other, Quaternion):
+        self.q -= other.q
+        self.p -= other.p
+        return self
+      else:
+        return NotImplemented
+        
+    def __neg__(self):                                                                      
+      """Unary positive operator"""  
+      self.q *= -1.0
+      self.p *= -1.0                                                                     
+      return self
+      
+      
+    def __mul__(self, other):
+      """Multiplication with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
+                              p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
+      elif isinstance(other, (int, float)):
+        return self.__class__(q=self.q*other,
+                              p=self.p*other)
+      else:
+        return NotImplemented
+
+    def __imul__(self, other):
+      """In-place multiplication with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        self.q = self.q*other.q - np.dot(self.p,other.p)
+        self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
+        return self
+      elif isinstance(other, (int, float)):
+        self *= other
+        return self
+      else:
+        return NotImplemented
+
+
+    def __truediv__(self, other):
+      """Divsion with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        s = other.conjugate()/abs(other)**2.
+        return self.__class__(q=self.q * s,
+                              p=self.p * s)
+      elif isinstance(other, (int, float)):
+        self.q /= other
+        self.p /= other
+        return self
+      else:
+        return NotImplemented
+
+    def __itruediv__(self, other):
+      """In-place divsion with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        s = other.conjugate()/abs(other)**2.
+        self *= s
+        return self
+      elif isinstance(other, (int, float)):
+        self.q /= other
+        return self
+      else:
+        return NotImplemented
+
+
+    def __pow__(self, exponent):
+      """Power"""
+      if isinstance(exponent, (int, float)):
+        omega = np.acos(self.q)
+        return self.__class__(q=         np.cos(exponent*omega),
+                              p=self.p * np.sin(exponent*omega)/np.sin(omega))
+      else:
+        return NotImplemented
+
+    def __ipow__(self, exponent):
+      """In-place power"""
+      if isinstance(exponent, (int, float)):
+        omega = np.acos(self.q)
+        self.q  = np.cos(exponent*omega)
+        self.p *= np.sin(exponent*omega)/np.sin(omega)
+      else:
+        return NotImplemented
+      
+
+    def __abs__(self):
+      """Norm"""
+      return np.sqrt(self.q ** 2 + np.dot(self.p,self.p))
+
+    magnitude = __abs__
+
+
+    def __eq__(self,other):
+      """Equal (sufficiently close) to each other"""
+      return np.isclose(( self-other).magnitude(),0.0) \
+          or np.isclose((-self-other).magnitude(),0.0)
+
+    def __ne__(self,other):
+      """Not equal (sufficiently close) to each other"""
+      return not self.__eq__(other)
+
+
+    def asM(self):
+      """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)"""
+      return np.outer(self.asArray(),self.asArray())
+
+    def asArray(self):
+      """As numpy array"""
+      return np.array((self.q,self.p[0],self.p[1],self.p[2]))
+ 
+    def asList(self):
+      return [self.q]+list(self.p)
+
+    def normalize(self):
+      d = self.magnitude()
+      if d > 0.0:
+          self.q /= d
+          self.p /= d
+      return self
+      
+    def normalized(self):
+      return self.copy().normalize()
+
+  
+    def conjugate(self):
+      self.p = -self.p
+      return self
+      
+    def conjugated(self):
+      return self.copy().conjugate()
+
+
+    def homomorph(self):
+      if self.q < 0.0:
+        self.q = -self.q
+        self.p = -self.p
+      return self
+      
+    def homomorphed(self):
+      return self.copy().homomorph()