From c8d774f456b17506874324f49a0dd1631195a6ec Mon Sep 17 00:00:00 2001
From: Harsha	Phukan <phykanha@egr.msu.edu>
Date: Fri, 22 May 2015 20:04:05 +0000
Subject: [PATCH] adopted math from Martin Baker to calculate quaternion from
 matrix representation.

---
 lib/damask/orientation.py | 46 ++++++++++++++++++++-------------------
 1 file changed, 24 insertions(+), 22 deletions(-)

diff --git a/lib/damask/orientation.py b/lib/damask/orientation.py
index c66db9c82..d88165b6d 100644
--- a/lib/damask/orientation.py
+++ b/lib/damask/orientation.py
@@ -428,50 +428,52 @@ class Quaternion:
         return cls([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[0,0] + m[1,1] + m[2,2] > 0.00000001:
-        t = m[0,0] + m[1,1] + m[2,2] + 1.0
-        s = 0.5/math.sqrt(t)
+      tr=m[0,0]+m[1,1]+m[2,2]
+      if tr > 0.00000001:
+        s = math.sqrt(tr + 1.0)*2.0
 
         return cls(
-          [ s*t,
-            (m[1,2] - m[2,1])*s,
-            (m[2,0] - m[0,2])*s,
-            (m[0,1] - m[1,0])*s,
+          [ s*0.25,
+            (m[2,1] - m[1,2])/s,
+            (m[0,2] - m[2,0])/s,
+            (m[1,0] - m[0,1])/s,
           ])
 
       elif m[0,0] > m[1,1] and m[0,0] > m[2,2]:
         t = m[0,0] - m[1,1] - m[2,2] + 1.0
-        s = 0.5/math.sqrt(t)
+        s = 2.0*math.sqrt(t)
 
         return cls(
-          [ (m[1,2] - m[2,1])*s,
-            s*t,
-            (m[0,1] + m[1,0])*s,
-            (m[2,0] + m[0,2])*s,
+          [ (m[2,1] - m[1,2])/s,
+            s*0.25,
+            (m[0,1] + m[1,0])/s,
+            (m[2,0] + m[0,2])/s,
           ])
 
       elif m[1,1] > m[2,2]:
         t = -m[0,0] + m[1,1] - m[2,2] + 1.0
-        s = 0.5/math.sqrt(t)
+        s = 2.0*math.sqrt(t)
 
         return cls(
-          [ (m[2,0] - m[0,2])*s,
-            (m[0,1] + m[1,0])*s,
-            s*t,
-            (m[1,2] + m[2,1])*s,
+          [ (m[0,2] - m[2,0])/s,
+            (m[0,1] + m[1,0])/s,
+            s*0.25,
+            (m[1,2] + m[2,1])/s,
           ])
 
       else:
         t = -m[0,0] - m[1,1] + m[2,2] + 1.0
-        s = 0.5/math.sqrt(t)
+        s = 2.0*math.sqrt(t)
 
         return cls(
-          [ (m[0,1] - m[1,0])*s,
-            (m[2,0] + m[0,2])*s,
-            (m[1,2] + m[2,1])*s,
-            s*t,
+          [ (m[1,0] - m[0,1])/s,
+            (m[2,0] + m[0,2])/s,
+            (m[1,2] + m[2,1])/s,
+            s*0.25,
           ])