r9325 - /1.3/maths_fns/rotation_matrix.py -- August 18, 2009 - 14:28

Author: bugman
Date: Tue Aug 18 14:28:08 2009
New Revision: 9325

URL: http://svn.gna.org/viewcvs/relax?rev=9325&view=rev
Log:
Bug fix for the quaternion_to_R() function.

The Wikipedia page was wrong (or at least misleading).  The calculation of 
the rotation matrix is
now correct.


Modified:
    1.3/maths_fns/rotation_matrix.py

Modified: 1.3/maths_fns/rotation_matrix.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/maths_fns/rotation_matrix.py?rev=9325&r1=9324&r2=9325&view=diff
==============================================================================
--- 1.3/maths_fns/rotation_matrix.py (original)
+++ 1.3/maths_fns/rotation_matrix.py Tue Aug 18 14:28:08 2009
@@ -36,7 +36,8 @@
         Q = |     2xy + 2zw      1 - 2x**2 - 2z**2      2yz - 2xw     |,
             |     2xz - 2yw          2yz + 2xw      1 - 2x**2 - 2y**2 |
 
-    and where the quaternion is defined as (x, y, z, w).
+    and where the quaternion is defined as q = (w, x, y, z).  This has been 
verified using Simo
+    Särkkä's "Notes on Quaternions" at http://www.lce.hut.fi/~ssarkka/.
 
 
     @param quat:    The quaternion.
@@ -45,28 +46,33 @@
     @type matrix:   numpy 3D, rank-2 array
     """
 
+    # Alias.
+    (w, x, y, z) = quat
+
     # Repetitive calculations.
-    q4_2 = quat[3]**2
-    q12 = quat[0] * quat[1]
-    q13 = quat[0] * quat[2]
-    q14 = quat[0] * quat[3]
-    q23 = quat[1] * quat[2]
-    q24 = quat[1] * quat[3]
-    q34 = quat[2] * quat[3]
+    x2 = 2.0 * x**2
+    y2 = 2.0 * y**2
+    z2 = 2.0 * z**2
+    xw = 2.0 * x*w
+    xy = 2.0 * x*y
+    xz = 2.0 * x*z
+    yw = 2.0 * y*w
+    yz = 2.0 * y*z
+    zw = 2.0 * z*w
 
     # The diagonal.
-    matrix[0, 0] = 2.0 * (quat[0]**2 + q4_2) - 1.0
-    matrix[1, 1] = 2.0 * (quat[1]**2 + q4_2) - 1.0
-    matrix[2, 2] = 2.0 * (quat[2]**2 + q4_2) - 1.0
-
-    # Off-diagonal.
-    matrix[0, 1] = 2.0 * (q12 - q34)
-    matrix[0, 2] = 2.0 * (q13 + q24)
-    matrix[1, 2] = 2.0 * (q23 - q14)
-
-    matrix[1, 0] = 2.0 * (q12 + q34)
-    matrix[2, 0] = 2.0 * (q13 - q24)
-    matrix[2, 1] = 2.0 * (q23 + q14)
+    matrix[0, 0] = 1.0 - y2 - z2
+    matrix[1, 1] = 1.0 - x2 - z2
+    matrix[2, 2] = 1.0 - x2 - y2
+
+    # The off-diagonal.
+    matrix[0, 1] = xy - zw
+    matrix[0, 2] = xz + yw
+    matrix[1, 2] = yz - xw
+
+    matrix[1, 0] = xy + zw
+    matrix[2, 0] = xz - yw
+    matrix[2, 1] = yz + xw
 
 
 def R_2vect(R, vector_orig, vector_fin):