r8893 - /1.3/data/align_tensor.py -- February 27, 2009 - 18:58

Author: bugman
Date: Fri Feb 27 18:58:24 2009
New Revision: 8893

URL: http://svn.gna.org/viewcvs/relax?rev=8893&view=rev
Log:
The diagonalised alignment tensor is now being created (and is being done 
correctly).


Modified:
    1.3/data/align_tensor.py

Modified: 1.3/data/align_tensor.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/data/align_tensor.py?rev=8893&r1=8892&r2=8893&view=diff
==============================================================================
--- 1.3/data/align_tensor.py (original)
+++ 1.3/data/align_tensor.py Fri Feb 27 18:58:24 2009
@@ -24,6 +24,7 @@
 from re import search
 from math import cos, sin
 from numpy import dot, float64, identity, transpose, zeros
+from numpy.linalg import eigvals
 from types import ListType
 
 # relax module imports.
@@ -505,8 +506,8 @@
     return tensor
 
 
-def calc_tensor_diag(rotation, tensor):
-    """Function for calculating the diagonalised alignment tensor.
+def calc_tensor_diag(tensor):
+    """Calculate the diagonalised alignment tensor.
 
     The diagonalised alignment tensor is defined as::
 
@@ -514,10 +515,8 @@
         tensor  =  |  0   Ayy'  0  |.
                    |  0    0   Azz'|
 
-    The diagonalised alignment tensor is calculated using the tensor and the 
rotation matrix
-    through the equation::
-
-        R^T . tensor_diag . R.
+    The diagonalised alignment tensor is calculated by eigenvalue 
decomposition.
+
 
     @param rotation:    The rotation matrix.
     @type rotation:     numpy array ((3, 3), float64)
@@ -527,8 +526,28 @@
     @rtype:             numpy array ((3, 3), float64)
     """
 
-    # Rotation (R^T . tensor_diag . R).
-    return dot(transpose(rotation), dot(tensor_diag, rotation))
+    # The eigenvalues.
+    vals = eigvals(tensor)
+
+    # Find the x < y < z indices.
+    abs_vals = abs(vals).tolist()
+    Axx_index = abs_vals.index(min(abs_vals))
+    Azz_index = abs_vals.index(max(abs_vals))
+    last_index = range(3)
+    last_index.pop(Axx_index)
+    last_index.pop(Azz_index)
+    Ayy_index = last_index[0]
+
+    # Empty tensor.
+    tensor_diag = zeros((3, 3), float64)
+
+    # Fill the elements.
+    tensor_diag[0, 0] = vals[Axx_index]
+    tensor_diag[1, 1] = vals[Ayy_index]
+    tensor_diag[2, 2] = vals[Azz_index]
+
+    # Return the tensor.
+    return tensor_diag
 
 
 def dependency_generator():
@@ -561,10 +580,10 @@
     yield ('Axx_unit',      ['alpha', 'beta', 'gamma'],                     
['alpha', 'beta', 'gamma'])
     yield ('Ayy_unit',      ['alpha', 'beta', 'gamma'],                     
['alpha', 'beta', 'gamma'])
     yield ('Azz_unit',      ['alpha', 'beta'],                              
['alpha', 'beta'])
-    yield ('tensor_diag',   ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'],            
['tensor', 'rotation'])
     yield ('rotation',      ['alpha', 'beta', 'gamma'],                     
['Axx_unit', 'Ayy_unit', 'Azz_unit'])
     yield ('tensor',        ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'],            
['Axx', 'Ayy', 'Azz', 'Axy', 'Axz', 'Ayz'])
     yield ('tensor_5D',     ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'],            
['Axx', 'Ayy', 'Azz', 'Axy', 'Axz', 'Ayz'])
+    yield ('tensor_diag',   ['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'],            
['tensor'])