removed Rodrigues class; removed quaternion "inverse" method; adopted Rowenhorst algorithm for angle--axis; use math.isclose(...,0.0) to check for ==0.0

lib/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, Rodrigues, Symmetry, Orientation # noqa
+from .orientation import Quaternion, Symmetry, Orientation # noqa
 
 #from .block       import Block           # only one class
 from .result      import Result           # noqa

lib/damask/orientation.py

@@ -7,24 +7,6 @@
 import math,os
 import numpy as np
 
-# ******************************************************************************************
-class Rodrigues:
-
-    def __init__(self, vector = np.zeros(3)):
-      self.vector = vector
-
-    def asQuaternion(self):
-      norm = np.linalg.norm(self.vector)
-      halfAngle = np.arctan(norm)
-      return Quaternion(np.cos(halfAngle),np.sin(halfAngle)*self.vector/norm)
-
-    def asAngleAxis(self):
-      norm = np.linalg.norm(self.vector)
-      halfAngle = np.arctan(norm)
-      return (2.0*halfAngle,self.vector/norm)
-
-
-
 # ******************************************************************************************
 class Quaternion:
     u"""
@@ -178,12 +160,13 @@ class Quaternion:
     magnitude = __abs__
 
     def __eq__(self,other):
-      """Equal at e-8 precision"""
+      """Equal (sufficiently close) to each other"""
-      return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8
+      return math.isclose(( self-other).magnitude(),0.0) \
+          or math.isclose((-self-other).magnitude(),0.0)
 
     def __ne__(self,other):
-      """Not equal at e-8 precision"""
+      """Not equal (sufficiently close) to each other"""
-      return not self.__eq__(self,other)
+      return not self.__eq__(other)
 
     def __cmp__(self,other):
       """Linear ordering"""
@@ -193,11 +176,6 @@ class Quaternion:
     def magnitude_squared(self):
       return self.q ** 2 + np.dot(self.p,self.p)
 
-    def identity(self):
-      self.q = 1.
-      self.p = np.zeros(3,dtype=float)
-      return self
-
     def normalize(self):
       d = self.magnitude()
       if d > 0.0:
@@ -209,13 +187,6 @@ class Quaternion:
       self.p = -self.p
       return self
 
-    def inverse(self):
-      d = self.magnitude()
-      if d > 0.0:
-        self.conjugate()
-        self /= d
-      return self
-
     def homomorph(self):
       if self.q < 0.0:
         self.q = -self.q
@@ -228,9 +199,6 @@ class Quaternion:
     def conjugated(self):
       return self.copy().conjugate()
 
-    def inversed(self):
-      return self.copy().inverse()
-
     def homomorphed(self):
       return self.copy().homomorph()
 
@@ -259,27 +227,24 @@ class Quaternion:
     def asAngleAxis(self,
                     degrees = False,
                     flat = False):
-      if self.q > 1.:
-          self.normalize()
 
-      s = math.sqrt(1. - self.q**2)
+      angle = 2.0*math.acos(self.q)
-      x = 2.*self.q**2 - 1.
-      y = 2.*self.q * s
 
-      angle = math.atan2(y,x)
+      if math.isclose(angle,0.0):
-      if angle < 0.0:
+        angle = 0.0
-        angle *= -1.
+        axis  = np.array([0.0,0.0,1.0])
-        s     *= -1.
+      elif math.isclose(self.q,0.0):
-
+        angle = math.pi
-      if flat:
+        axis  = self.p
-        return np.hstack((np.degrees(angle) if degrees else angle,
-                          np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s)))
       else:
-        return (np.degrees(angle) if degrees else angle,
+        axis  = np.sign(self.q)*self.p/np.linalg.norm(self.p)
-                np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s))
+
+      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 self.q == 0.0 else self.p/self.q
+      return np.inf*np.ones(3) if math.isclose(self.q,0.0) else self.p/self.q
 
     def asEulers(self,
                  degrees = False):
@@ -290,9 +255,9 @@ class Quaternion:
       q12 = self.p[0]**2 + self.p[1]**2
       chi = np.sqrt(q03*q12)
       
-      if abs(chi) < 1e-10 and abs(q12) < 1e-10:
+      if math.isclose(chi,0.0) and math.isclose(q12,0.0):
         eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0])
-      elif abs(chi) < 1e-10 and abs(q03) < 1e-10:
+      elif math.isclose(chi,0.0) and math.isclose(q03,0.0):
         eulers = np.array([math.atan2( 2  *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0])
       else:
         eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi),
@@ -456,11 +421,11 @@ class Symmetry:
 
 
   def __eq__(self, other):
-    """Equal"""
+    """Equal to other"""
     return self.lattice == other.lattice
 
   def __neq__(self, other):
-    """Not equal"""
+    """Not equal to other"""
     return not self.__eq__(other)
 
   def __cmp__(self,other):
@@ -549,7 +514,7 @@ class Symmetry:
 
   def inFZ(self,R):
     """Check whether given Rodrigues vector falls into fundamental zone of own symmetry."""
-    if isinstance(R, Quaternion): R = R.asRodrigues()                                               # translate accidentially passed quaternion
+    if isinstance(R, Quaternion): R = R.asRodrigues()                                               # translate accidentally passed quaternion
 # fundamental zone in Rodrigues space is point symmetric around origin
     R = abs(R)                                                                                      
     if self.lattice == 'cubic':