r11021 - /1.3/generic_fns/diffusion_tensor.py -- March 18, 2010 - 15:09

Author: bugman
Date: Thu Mar 18 15:09:00 2010
New Revision: 11021

URL: http://svn.gna.org/viewcvs/relax?rev=11021&view=rev
Log:
Final bug fix for the ellipsoid() function for the (Dxx, Dyy, Dzz, Dxy, Dxz, 
Dyz) parameter set.

The 2304 unique Euler angle values for a rank-2 tensor with no skew (12 axis 
orders, 2 frames of
ref., 2 directions (forwards + reverse), 4 rotational symmetries (pi 
rotations about each
eigenvector), 2 handednesses (left + right), and 6 eigenvalue orders) are now 
all properly handled.
The resultant Euler angle values are a single set of the 2304, matching that 
already used in relax
(zyz, rotated axis system, reverse direction, rotational symmetry collapse to 
0 <= alpha, beta,
gamma <= pi, right handedness, and Dx <= Dy <= Dz).


Modified:
    1.3/generic_fns/diffusion_tensor.py

Modified: 1.3/generic_fns/diffusion_tensor.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/diffusion_tensor.py?rev=11021&r1=11020&r2=11021&view=diff
==============================================================================
--- 1.3/generic_fns/diffusion_tensor.py (original)
+++ 1.3/generic_fns/diffusion_tensor.py Thu Mar 18 15:09:00 2010
@@ -26,8 +26,8 @@
 # Python module imports.
 from copy import deepcopy
 from math import cos, pi, sin
-from numpy import float64, zeros
-from numpy.linalg import eig
+from numpy import cross, float64, transpose, zeros
+from numpy.linalg import eig, norm
 from re import search
 
 # relax module imports.
@@ -376,14 +376,32 @@
         reorder = zeros(3, int)
         Di_sort = sorted(Di)
         Di = Di.tolist()
-        R_sort = zeros((3, 3), float64)
+        R_new = zeros((3, 3), float64)
 
         # Reorder columns.
         for i in range(3):
-            R_sort[:, i] = R[:, Di.index(Di_sort[i])]
-
-        # Euler angles.
-        alpha, beta, gamma = R_to_euler_zyz(R_sort)
+            R_new[:, i] = R[:, Di.index(Di_sort[i])]
+
+        # Switch from the left handed to right handed universes (if needed).
+        if norm(cross(R_new[:, 0], R_new[:, 1]) - R_new[:, 2]) > 1e-7:
+            R_new[:, 2] = -R_new[:, 2]
+        
+        # Reverse the rotation.
+        R_new = transpose(R_new)
+
+        # Euler angles (reverse rotation in the rotated axis system).
+        gamma, beta, alpha = R_to_euler_zyz(R_new)
+
+        # Collapse the pi axis rotation symmetries.
+        if alpha >= pi:
+            alpha = alpha - pi
+        if gamma >= pi:
+            alpha = pi - alpha
+            beta = pi - beta
+            gamma = gamma - pi
+        if beta >= pi:
+            alpha = pi - alpha
+            beta = beta - pi
 
         # Set the parameters.
         set(value=[Di_sort[0], Di_sort[1], Di_sort[2]], param=['Dx', 'Dy', 
'Dz'])