using new Orientation class

.gitlab-ci.yml

@@ -207,7 +207,9 @@ Post_ParaviewRelated:
 
 Post_OrientationConversion:
   stage: postprocessing
-  script: OrientationConversion/test.py
+  script:
+    - OrientationConversion/test.py
+    - OrientationConversion/test2.py
   except:
     - master
     - release

PRIVATE

@@ -1 +1 @@
-Subproject commit f0090997df817f0a0b5a480a60e81929875b1010
+Subproject commit 8deb37dd4526fb5e1425fe1d2360508d01b6ac3e

processing/post/addIPFcolor.py

@@ -41,6 +41,10 @@ parser.set_defaults(pole = (0.0,0.0,1.0),
 
 (options, filenames) = parser.parse_args()
 
+# damask.Orientation requires Bravais lattice, but we are only interested in symmetry
+symmetry2lattice={'cubic':'bcc','hexagonal':'hex','tetragonal':'bct'}
+lattice = symmetry2lattice[options.symmetry]
+
 pole = np.array(options.pole)
 pole /= np.linalg.norm(pole)
 
@@ -78,8 +82,8 @@ for name in filenames:
 
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
+    o = damask.Orientation(np.array(list(map(float,table.data[column:column+4]))),
-                           symmetry   = options.symmetry).reduced()
+                           lattice   = lattice).reduced()
 
     table.data_append(o.IPFcolor(pole))
     outputAlive = table.data_write()                                                                # output processed line

processing/post/addPole.py

@@ -75,9 +75,9 @@ for name in filenames:
 # ------------------------------------------ process data ------------------------------------------
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
+    o = damask.Rotation(np.array(list(map(float,table.data[column:column+4]))))
 
-    rotatedPole = o.quaternion*pole                                                                 # rotate pole according to crystal orientation
+    rotatedPole = o*pole                                                                 # rotate pole according to crystal orientation
     (x,y) = rotatedPole[0:2]/(1.+abs(pole[2]))                                                      # stereographic projection
 
     table.data_append([np.sqrt(x*x+y*y),np.arctan2(y,x)] if options.polar else [x,y])               # cartesian coordinates

processing/post/rotateData.py

@@ -31,7 +31,7 @@ parser.add_option('--degrees',
                   action = 'store_true',
                   help = 'angles are given in degrees')
 
-parser.set_defaults(rotation = (0.,1.,0.,0.),                                                       # no rotation about 1,0,0
+parser.set_defaults(rotation = (0.,1.,1.,1.),                                                       # no rotation about 1,1,1
                     degrees = False,
                    )
                     

processing/pre/geom_fromTable.py

@@ -52,11 +52,8 @@ parser.add_option('--crystallite',
 
 
 parser.set_defaults(symmetry       = [damask.Symmetry.lattices[-1]],
-                    tolerance      = 0.0,
-                    degrees        = False,
                     homogenization = 1,
                     crystallite    = 1,
-                    verbose        = False,
                     pos            = 'pos',
                    )
 

python/damask/__init__.py

@@ -13,7 +13,7 @@ from .asciitable  import ASCIItable       # noqa
     
 from .config      import Material         # noqa
 from .colormaps   import Colormap, Color  # noqa
-from .orientation import Quaternion, Symmetry, Rotation, Orientation # noqa
+from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa
 
 #from .block       import Block           # only one class
 from .result      import Result           # noqa

python/damask/orientation.py

@@ -1,13 +1,13 @@
 # -*- coding: UTF-8 no BOM -*-
 
-import math,os
+import math
 import numpy as np
 from . import Lambert
 
 P = -1
 
 ####################################################################################################
-class Quaternion2:
+class Quaternion:
     u"""
     Quaternion with basic operations
     
@@ -50,7 +50,7 @@ class Quaternion2:
 
     def __add__(self, other):
       """Addition"""
-      if isinstance(other, Quaternion2):
+      if isinstance(other, Quaternion):
         return self.__class__(q=self.q + other.q,
                               p=self.p + other.p)
       else:
@@ -58,7 +58,7 @@ class Quaternion2:
 
     def __iadd__(self, other):
       """In-place addition"""
-      if isinstance(other, Quaternion2):
+      if isinstance(other, Quaternion):
         self.q += other.q
         self.p += other.p
         return self
@@ -72,7 +72,7 @@ class Quaternion2:
       
     def __sub__(self, other):
       """Subtraction"""
-      if isinstance(other, Quaternion2):
+      if isinstance(other, Quaternion):
         return self.__class__(q=self.q - other.q,
                               p=self.p - other.p)
       else:
@@ -80,7 +80,7 @@ class Quaternion2:
 
     def __isub__(self, other):
       """In-place subtraction"""
-      if isinstance(other, Quaternion2):
+      if isinstance(other, Quaternion):
         self.q -= other.q
         self.p -= other.p
         return self
@@ -96,7 +96,7 @@ class Quaternion2:
       
     def __mul__(self, other):
       """Multiplication with quaternion or scalar"""
-      if isinstance(other, Quaternion2):
+      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)):
@@ -107,7 +107,7 @@ class Quaternion2:
 
     def __imul__(self, other):
       """In-place multiplication with quaternion or scalar"""
-      if isinstance(other, Quaternion2):
+      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
@@ -120,7 +120,7 @@ class Quaternion2:
 
     def __truediv__(self, other):
       """Divsion with quaternion or scalar"""
-      if isinstance(other, Quaternion2):
+      if isinstance(other, Quaternion):
         s = other.conjugate()/abs(other)**2.
         return self.__class__(q=self.q * s,
                               p=self.p * s)
@@ -133,7 +133,7 @@ class Quaternion2:
 
     def __itruediv__(self, other):
       """In-place divsion with quaternion or scalar"""
-      if isinstance(other, Quaternion2):
+      if isinstance(other, Quaternion):
         s = other.conjugate()/abs(other)**2.
         self *= s
         return self
@@ -215,7 +215,8 @@ class Rotation:
     u"""
     Orientation stored as unit quaternion.
     
-    All methods and naming conventions based on Rowenhorst_etal2015
+    Following: D Rowenhorst et al. Consistent representations of and conversions between 3D rotations
+               10.1088/0965-0393/23/8/083501
     Convention 1: coordinate frames are right-handed
     Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
                   when viewing from the end point of the rotation axis towards the origin
@@ -242,10 +243,10 @@ class Rotation:
       If a quaternion is given, it needs to comply with the convection. Use .fromQuaternion
       to check the input.
       """
-      if isinstance(quaternion,Quaternion2):
+      if isinstance(quaternion,Quaternion):
         self.quaternion = quaternion.copy()
       else:
-        self.quaternion = Quaternion2(q=quaternion[0],p=quaternion[1:4])
+        self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4])
         self.quaternion.homomorph() # ToDo: Needed?
 
     def __repr__(self):
@@ -455,373 +456,6 @@ class Rotation:
       return self.__class__(other.quaternion*self.quaternion.conjugated())
       
 
-# ******************************************************************************************
-class Quaternion:
-    u"""
-    Orientation represented as unit quaternion.
-    
-    All methods and naming conventions based on Rowenhorst_etal2015
-    Convention 1: coordinate frames are right-handed
-    Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
-                  when viewing from the end point of the rotation axis towards the origin
-    Convention 3: rotations will be interpreted in the passive sense
-    Convention 4: Euler angle triplets are implemented using the Bunge convention,
-                  with the angular ranges as [0, 2π],[0, π],[0, 2π]
-    Convention 5: the rotation angle ω is limited to the interval [0, π]
-    Convention 6: P = -1 (as default)
-    
-    w is the real part, (x, y, z) are the imaginary parts.
-
-    Vector "a" (defined in coordinate system "A") is passively rotated
-               resulting in new coordinates "b" when expressed in system "B".
-    b = Q * a
-    b = np.dot(Q.asMatrix(),a)
-    """
-
-    def __init__(self,
-                 quat = None,
-                 q    =  1.0,
-                 p    = np.zeros(3,dtype=float)):
-      """Initializes to identity unless specified"""
-      self.q = quat[0] if quat is not None else q
-      self.p = np.array(quat[1:4]) if quat is not None else p
-      self.homomorph()
-
-    def __iter__(self):
-      """Components"""
-      return iter(self.asList())
-
-    def __copy__(self):
-      """Copy"""
-      return self.__class__(q=self.q,p=self.p.copy())
-
-    copy = __copy__
-
-    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 __pow__(self, exponent):
-      """Power"""
-      omega = math.acos(self.q)
-      return self.__class__(q=         math.cos(exponent*omega),
-                            p=self.p * math.sin(exponent*omega)/math.sin(omega))
-
-    def __ipow__(self, exponent):
-      """In-place power"""
-      omega = math.acos(self.q)
-      self.q  = math.cos(exponent*omega)
-      self.p *= math.sin(exponent*omega)/math.sin(omega)
-      return self
-
-    def __mul__(self, other):
-      """Multiplication"""
-      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
-      P = -1.0
-      try:                                                           # 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))
-      except: pass
-      try:                                                           # vector (perform passive rotation)
-          ( x, y, z) = self.p
-          (Vx,Vy,Vz) = other[0:3]
-          A = self.q*self.q - np.dot(self.p,self.p)
-          B = 2.0 * (x*Vx + y*Vy + z*Vz)
-          C = 2.0 * P*self.q
-
-          return np.array([
-            A*Vx + B*x + C*(y*Vz - z*Vy),
-            A*Vy + B*y + C*(z*Vx - x*Vz),
-            A*Vz + B*z + C*(x*Vy - y*Vx),
-            ])
-      except: pass
-      try:                                                           # scalar
-          return self.__class__(q=self.q*other,
-                                p=self.p*other)
-      except:
-          return self.copy()
-
-    def __imul__(self, other):
-      """In-place multiplication"""
-      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
-      P = -1.0
-      try:                                                           # 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)
-      except: pass
-      return self
-
-    def __div__(self, other):
-      """Division"""
-      if isinstance(other, (int,float)):
-        return self.__class__(q=self.q / other,
-                              p=self.p / other)
-      else:
-          return NotImplemented
-
-    def __idiv__(self, other):
-      """In-place division"""
-      if isinstance(other, (int,float)):
-          self.q /= other
-          self.p /= other
-      return self
-
-    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
-
-    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
-
-    def __neg__(self):
-      """Additive inverse"""
-      self.q = -self.q
-      self.p = -self.p
-      return self
-
-    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 __cmp__(self,other):
-      """Linear ordering"""
-      return (1 if np.linalg.norm(self.asRodrigues()) > np.linalg.norm(other.asRodrigues()) else 0) \
-           - (1 if np.linalg.norm(self.asRodrigues()) < np.linalg.norm(other.asRodrigues()) else 0)
-
-    def magnitude_squared(self):
-      return self.q ** 2 + np.dot(self.p,self.p)
-
-    def normalize(self):
-      d = self.magnitude()
-      if d > 0.0:
-          self.q /= d
-          self.p /= d
-      return self
-
-    def conjugate(self):
-      self.p = -self.p
-      return self
-
-    def homomorph(self):
-      if self.q < 0.0:
-        self.q = -self.q
-        self.p = -self.p
-      return self
-
-    def normalized(self):
-      return self.copy().normalize()
-
-    def conjugated(self):
-      return self.copy().conjugate()
-
-    def homomorphed(self):
-      return self.copy().homomorph()
-
-    def asList(self):
-      return [self.q]+list(self.p)
-
-    def asM(self):                                                   # to find Averaging Quaternions (see F. Landis Markley et al.)
-      return np.outer(self.asList(),self.asList())
-      
-    def asMatrix(self):
-      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
-      P = -1.0
-      qbarhalf = 0.5*(self.q**2 - np.dot(self.p,self.p))
-      return 2.0*np.array(
-        [[        qbarhalf + self.p[0]**2          ,
-           self.p[0]*self.p[1] -P* self.q*self.p[2],
-           self.p[0]*self.p[2] +P* self.q*self.p[1] ],
-         [ self.p[0]*self.p[1] +P* self.q*self.p[2],
-                  qbarhalf + self.p[1]**2          ,
-           self.p[1]*self.p[2] -P* self.q*self.p[0] ],
-         [ self.p[0]*self.p[2] -P* self.q*self.p[1],
-           self.p[1]*self.p[2] +P* self.q*self.p[0],
-                  qbarhalf + self.p[2]**2           ],
-        ])
-
-    def asAngleAxis(self,
-                    degrees = False,
-                    flat = False):
-
-      angle = 2.0*math.acos(self.q)
-
-      if np.isclose(angle,0.0):
-        angle = 0.0
-        axis  = np.array([0.0,0.0,1.0])
-      elif np.isclose(self.q,0.0):
-        angle = math.pi
-        axis  = self.p
-      else:
-        axis  = np.sign(self.q)*self.p/np.linalg.norm(self.p)
-
-      angle = np.degrees(angle) if degrees else angle
-
-      return np.hstack((angle,axis)) if flat else (angle,axis)
-
-    def asRodrigues(self):
-      return np.inf*np.ones(3) if np.isclose(self.q,0.0) else self.p/self.q
-
-
-#    # Static constructors
-    @classmethod
-    def fromIdentity(cls):
-      return cls()
-
-
-    @classmethod
-    def fromRandom(cls,randomSeed = None):
-      import binascii
-      if randomSeed is None:
-        randomSeed = int(binascii.hexlify(os.urandom(4)),16)
-      np.random.seed(randomSeed)
-      r = np.random.random(3)
-      A = math.sqrt(max(0.0,r[2]))
-      B = math.sqrt(max(0.0,1.0-r[2]))
-      w = math.cos(2.0*math.pi*r[0])*A
-      x = math.sin(2.0*math.pi*r[1])*B
-      y = math.cos(2.0*math.pi*r[1])*B
-      z = math.sin(2.0*math.pi*r[0])*A
-      return cls(quat=[w,x,y,z])
-
-
-    @classmethod
-    def fromRodrigues(cls, rodrigues):
-      if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues)
-      norm = np.linalg.norm(rodrigues)
-      halfangle = math.atan(norm)
-      s = math.sin(halfangle)
-      c = math.cos(halfangle)
-      return cls(q=c,p=s*rodrigues/norm)
-
-
-    @classmethod
-    def fromAngleAxis(cls,
-                      angle,
-                      axis,
-                      degrees = False):
-      if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype=float)
-      axis = axis.astype(float)/np.linalg.norm(axis)
-      angle = np.radians(angle) if degrees else angle
-      s = math.sin(0.5 * angle)
-      c = math.cos(0.5 * angle)
-      return cls(q=c,p=axis*s)
-
-
-    @classmethod
-    def fromEulers(cls,
-                   eulers,
-                   degrees = False):
-      if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype=float)
-      eulers = np.radians(eulers) if degrees else eulers
-
-      sigma = 0.5*(eulers[0]+eulers[2])
-      delta = 0.5*(eulers[0]-eulers[2])
-      c = np.cos(0.5*eulers[1])
-      s = np.sin(0.5*eulers[1])
-
-      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
-      P = -1.0
-      w =       c * np.cos(sigma) 
-      x =  -P * s * np.cos(delta) 
-      y =  -P * s * np.sin(delta)
-      z =  -P * c * np.sin(sigma) 
-      return cls(quat=[w,x,y,z])
-
-
-# Modified Method to calculate Quaternion from Orientation Matrix,
-# Source: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
-
-    @classmethod
-    def fromMatrix(cls, m):
-      if m.shape != (3,3) and np.prod(m.shape) == 9:
-        m = m.reshape(3,3)
-
-      # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
-      P = -1.0
-      w =   0.5*math.sqrt(max(0.0,1.0+m[0,0]+m[1,1]+m[2,2]))
-      x = P*0.5*math.sqrt(max(0.0,1.0+m[0,0]-m[1,1]-m[2,2]))
-      y = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]+m[1,1]-m[2,2]))
-      z = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]-m[1,1]+m[2,2]))
-
-      x *= -1 if m[2,1] < m[1,2] else 1
-      y *= -1 if m[0,2] < m[2,0] else 1
-      z *= -1 if m[1,0] < m[0,1] else 1
-
-      return cls(quat=np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2))
-
-
-    @classmethod
-    def new_interpolate(cls, q1, q2, t):
-        """
-        Interpolation
- 
-        See http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf
-        for (another?) way to interpolate quaternions.
-        """
-        assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion)
-        Q = cls()
-
-        costheta = q1.q*q2.q + np.dot(q1.p,q2.p)
-        if costheta < 0.:
-            costheta = -costheta
-            q1 = q1.conjugated()
-        elif costheta > 1.:
-            costheta = 1.
-
-        theta = math.acos(costheta)
-        if abs(theta) < 0.01:
-            Q.q = q2.q
-            Q.p = q2.p
-            return Q
-
-        sintheta = math.sqrt(1.0 - costheta * costheta)
-        if abs(sintheta) < 0.01:
-            Q.q = (q1.q + q2.q) * 0.5
-            Q.p = (q1.p + q2.p) * 0.5
-            return Q
-
-        ratio1 = math.sin((1.0 - t) * theta) / sintheta
-        ratio2 = math.sin(       t  * theta) / sintheta
-
-        Q.q = q1.q * ratio1 + q2.q * ratio2
-        Q.p = q1.p * ratio1 + q2.p * ratio2
-        return Q
-
-
 # ******************************************************************************************
 class Symmetry:
   """
@@ -932,26 +566,16 @@ class Symmetry:
                     [ 1.0,0.0,0.0,0.0 ],
                   ]
 
-    return list(map(Quaternion,
+    return np.array(symQuats)
-                    np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))]))
-    
-    
-  def equivalentQuaternions(self,
-                            quaternion,
-                            who = []):
-    """List of symmetrically equivalent quaternions based on own symmetry."""
-    return [q*quaternion for q in self.symmetryQuats(who)]
     
 
   def inFZ(self,R):
-    """Check whether given Rodrigues vector falls into fundamental zone of own symmetry."""
+    """
-    if isinstance(R, Quaternion): R = R.asRodrigues()                                               # translate accidentally passed quaternion
+    Check whether given Rodrigues vector falls into fundamental zone of own symmetry.
-# fundamental zone in Rodrigues space is point symmetric around origin
     
-    if R.shape[0]==4: # transition old (length not stored separately) to new
+    Fundamental zone in Rodrigues space is point symmetric around origin.
+    """
     Rabs = abs(R[0:3]*R[3])
-    else:
-      Rabs = abs(R)
 
     if self.lattice == 'cubic':
       return     math.sqrt(2.0)-1.0 >= Rabs[0] \
@@ -1367,27 +991,30 @@ class Lattice:
 
     relationship = models[model]
 
-    r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()),
+    r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()),                    # target lattice
          'rotations':[] }
 
     myPlane_id    = relationship['mapping'][self.lattice]
     otherPlane_id = (myPlane_id+1)%2
     myDir_id      = myPlane_id +2
     otherDir_id   = otherPlane_id +2
+
     for miller in np.hstack((relationship['planes'],relationship['directions'])):
       myPlane     = miller[myPlane_id]/    np.linalg.norm(miller[myPlane_id])
       myDir       = miller[myDir_id]/      np.linalg.norm(miller[myDir_id])
+      myMatrix    = np.array([myDir,np.cross(myPlane,myDir),myPlane]).T
+
       otherPlane  = miller[otherPlane_id]/ np.linalg.norm(miller[otherPlane_id])
       otherDir    = miller[otherDir_id]/   np.linalg.norm(miller[otherDir_id])
-      
-      myMatrix     = np.array([myDir,np.cross(myPlane,myDir),myPlane]).T
       otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane]).T
+
       r['rotations'].append(Rotation.fromMatrix(np.dot(otherMatrix,myMatrix.T)))
+
     return r
 
 
 
-class Orientation2:
+class Orientation:
   """
   Crystallographic orientation
   
@@ -1431,27 +1058,23 @@ class Orientation2:
     
     for i,sA in enumerate(mySymEqs):
       for j,sB in enumerate(otherSymEqs):
-        theQ = sB.rotation*mis*sA.rotation.inversed()
+        r = sB.rotation*mis*sA.rotation.inversed()
         for k in range(2):
-          theQ.inversed()
+          r.inversed()
-          breaker = self.lattice.symmetry.inFZ(theQ.asRodrigues()) \
+          breaker = self.lattice.symmetry.inFZ(r.asRodrigues()) \
-                    and (not SST or other.lattice.symmetry.inDisorientationSST(theQ.asRodrigues()))
+                    and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues()))
           if breaker: break
         if breaker: break
       if breaker: break
 
-# disorientation, own sym, other sym, self-->other: True, self<--other: False
+    return r
-    return theQ
 
   def inFZ(self):
     return self.lattice.symmetry.inFZ(self.rotation.asRodrigues())
   
   def equivalentOrientations(self):
     """List of orientations which are symmetrically equivalent"""
-    q = self.lattice.symmetry.symmetryQuats()
+    return [self.__class__(q*self.rotation.quaternion,self.lattice) for q in self.lattice.symmetry.symmetryQuats()]
-    q2 = [Quaternion2(q=a.asList()[0],p=a.asList()[1:4]) for a in q] # convert Quaternion to Quaternion2
-    x = [self.__class__(q3*self.rotation.quaternion,self.lattice) for q3 in q2]
-    return x
     
   def relatedOrientations(self,model):
     """List of orientations related by the given orientation relationship"""
@@ -1465,176 +1088,68 @@ class Orientation2:
 
     return self.__class__(me.rotation,self.lattice)
     
-# ******************************************************************************************
-class Orientation:
-
-  __slots__ = ['quaternion','symmetry']
-
-  def __init__(self,
-               quaternion = Quaternion.fromIdentity(),
-               Rodrigues  = None,
-               Eulers     = None,
-               random     = False,                                                                  # integer to have a fixed seed or True for real random
-               symmetry   = None,
-               degrees    = False,
-              ):
-    if random:                                                                                      # produce random orientation
-      if isinstance(random, bool ):
-        self.quaternion = Quaternion.fromRandom()
-      else:
-        self.quaternion = Quaternion.fromRandom(randomSeed=random)
-    elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,):                                   # based on given Euler angles
-      self.quaternion = Quaternion.fromEulers(Eulers,degrees=degrees)
-    elif isinstance(Rodrigues, np.ndarray) and Rodrigues.shape == (3,):                             # based on given Rodrigues vector
-      self.quaternion = Quaternion.fromRodrigues(Rodrigues)
-    elif isinstance(quaternion, Quaternion):                                                        # based on given quaternion
-      self.quaternion = quaternion.homomorphed()
-    elif (isinstance(quaternion, np.ndarray) and quaternion.shape == (4,)) or \
-         (isinstance(quaternion, list)       and len(quaternion)  ==  4  ):                         # based on given quaternion-like array
-      self.quaternion = Quaternion(quat=quaternion).homomorphed()
-
-    self.symmetry = Symmetry(symmetry)
-
-  def __copy__(self):
-    """Copy"""
-    return self.__class__(quaternion=self.quaternion,symmetry=self.symmetry.lattice)
-
-  copy = __copy__
-
-
-  def __repr__(self):
-    """Value as all implemented representations"""
-    return '\n'.join([
-              'Symmetry: {}'.format(self.symmetry),
-              'Quaternion: {}'.format(self.quaternion),
-              'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
-              ])
-
-  def asRodrigues(self):
-    return self.quaternion.asRodrigues()
-  rodrigues = property(asRodrigues)
-
-  def asAngleAxis(self,
-                  degrees = False,
-                  flat = False):
-    return self.quaternion.asAngleAxis(degrees,flat)
-
-  def asMatrix(self):
-    return self.quaternion.asMatrix()
-  matrix = property(asMatrix)
-
-  def inFZ(self):
-    return self.symmetry.inFZ(self.quaternion.asRodrigues())
-
-  def equivalentQuaternions(self,
-                            who = []):
-    return self.symmetry.equivalentQuaternions(self.quaternion,who)
-
-  def equivalentOrientations(self,
-                             who = []):
-    return [Orientation(quaternion = q, symmetry = self.symmetry.lattice) for q in self.equivalentQuaternions(who)]
-
-  def reduced(self):
-    """Transform orientation to fall into fundamental zone according to symmetry"""
-    for me in self.symmetry.equivalentQuaternions(self.quaternion):
-      if self.symmetry.inFZ(me.asRodrigues()): break
-
-    return Orientation(quaternion=me,symmetry=self.symmetry.lattice)
-
-
-  def disorientation(self,
-                     other,
-                     SST = True):
-    """
-    Disorientation between myself and given other orientation.
-
-    Rotation axis falls into SST if SST == True.
-    (Currently requires same symmetry for both orientations.
-     Look into A. Heinz and P. Neumann 1991 for cases with differing sym.)
-    """
-    if self.symmetry != other.symmetry:
-      raise NotImplementedError('disorientation between different symmetry classes not supported yet.')
-
-    misQ = other.quaternion*self.quaternion.conjugated()
-    mySymQs    =  self.symmetry.symmetryQuats() if SST else self.symmetry.symmetryQuats()[:1]       # take all or only first sym operation
-    otherSymQs = other.symmetry.symmetryQuats()
-    
-    for i,sA in enumerate(mySymQs):
-      for j,sB in enumerate(otherSymQs):
-        theQ = sB*misQ*sA.conjugated()
-        for k in range(2):
-          theQ.conjugate()
-          breaker = self.symmetry.inFZ(theQ) \
-                    and (not SST or other.symmetry.inDisorientationSST(theQ))
-          if breaker: break
-        if breaker: break
-      if breaker: break
-
-# disorientation, own sym, other sym, self-->other: True, self<--other: False
-    return (Orientation(quaternion = theQ,symmetry = self.symmetry.lattice),
-            i,j, k == 1)                                                                             
-
-
   def inversePole(self,
                   axis,
                   proper = False,
                   SST = True):
     """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)"""
     if SST:                                                                                         # pole requested to be within SST
-      for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)):                   # test all symmetric equivalent quaternions
+      for i,o in enumerate(self.equivalentOrientations()):                                          # test all symmetric equivalent quaternions
-        pole = q*axis                                                                               # align crystal direction to axis
+        pole = o.rotation*axis                                                                      # align crystal direction to axis
-        if self.symmetry.inSST(pole,proper): break                                                  # found SST version
+        if self.lattice.symmetry.inSST(pole,proper): break                                          # found SST version
     else:
-      pole = self.quaternion*axis                                                                   # align crystal direction to axis
+      pole = self.rotation*axis                                                                     # align crystal direction to axis
 
     return (pole,i if SST else 0)
     
+    
   def IPFcolor(self,axis):
     """TSL color of inverse pole figure for given axis"""
     color = np.zeros(3,'d')
 
-    for q in self.symmetry.equivalentQuaternions(self.quaternion):
+    for o in self.equivalentOrientations():
-      pole = q*axis                                                                                 # align crystal direction to axis
+      pole = o.rotation*axis                                                                        # align crystal direction to axis
-      inSST,color = self.symmetry.inSST(pole,color=True)
+      inSST,color = self.lattice.symmetry.inSST(pole,color=True)
       if inSST: break
 
     return color                                                                        
 
-  @classmethod
-  def average(cls,
-              orientations,
-              multiplicity = []):
-    """
-    Average orientation
 
-    ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman.
+  # @classmethod
-         Averaging Quaternions,
+  # def average(cls,
-         Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197.
+              # orientations,
-         doi: 10.2514/1.28949
+              # multiplicity = []):
-    usage:
+    # """
-         a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal')
+    # Average orientation
-         b = Orientation(Eulers=np.radians([20, 0, 0]),  symmetry='hexagonal')
-         avg = Orientation.average([a,b])
-    """
-    if not all(isinstance(item, Orientation) for item in orientations):
-      raise TypeError("Only instances of Orientation can be averaged.")
 
-    N = len(orientations)
+    # ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman.
-    if multiplicity == [] or not multiplicity:
+         # Averaging Quaternions,
-      multiplicity = np.ones(N,dtype='i')
+         # Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197.
+         # doi: 10.2514/1.28949
+    # usage:
+         # a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal')
+         # b = Orientation(Eulers=np.radians([20, 0, 0]),  symmetry='hexagonal')
+         # avg = Orientation.average([a,b])
+    # """
+    # if not all(isinstance(item, Orientation) for item in orientations):
+      # raise TypeError("Only instances of Orientation can be averaged.")
 
-    reference = orientations[0]                                                                     # take first as reference
+    # N = len(orientations)
-    for i,(o,n) in enumerate(zip(orientations,multiplicity)):
+    # if multiplicity == [] or not multiplicity:
-      closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0]             # select sym orientation with lowest misorientation
+      # multiplicity = np.ones(N,dtype='i')
-      M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n              # noqa add (multiples) of this orientation to average noqa
-    eig, vec = np.linalg.eig(M/N)
 
-    return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])),
+    # reference = orientations[0]                                                                     # take first as reference
-                       symmetry = reference.symmetry.lattice)
+    # for i,(o,n) in enumerate(zip(orientations,multiplicity)):
+      # closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0]             # select sym orientation with lowest misorientation
+      # M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n              # noqa add (multiples) of this orientation to average noqa
+    # eig, vec = np.linalg.eig(M/N)
+
+    # return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])),
+                       # symmetry = reference.symmetry.lattice)
 
 
 ####################################################################################################
-# Code below available according to the followin conditions on https://github.com/MarDiehl/3Drotations
+# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations
 ####################################################################################################
 # Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
 # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
@@ -1769,7 +1284,6 @@ def qu2eu(qu):
     eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0,   0.0]) if iszero(q12) else \
          np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2),         np.pi, 0.0])
   else:
-    #chiInv = 1.0/chi         ToDo: needed for what?
     eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ),
                    np.arctan2( 2.0*chi, q03-q12 ), 
                    np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )])