r23619 - /branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor/frame_order.py -- June 03, 2014 - 11:24

Author: bugman
Date: Tue Jun  3 11:24:59 2014
New Revision: 23619

URL: http://svn.gna.org/viewcvs/relax?rev=23619&view=rev
Log:
Updated the rotor CaM test data frame_order.py script for the parameter 
reduction.

The rotor axis {theta, phi} polar angles have been replaced by the single 
axis alpha angle.  This
now matches the script for the 2nd rotor model.


Modified:
    
branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor/frame_order.py

Modified: 
branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor/frame_order.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor/frame_order.py?rev=23619&r1=23618&r2=23619&view=diff
==============================================================================
--- 
branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor/frame_order.py
    (original)
+++ 
branches/frame_order_cleanup/test_suite/shared_data/frame_order/cam/rotor/frame_order.py
    Tue Jun  3 11:24:59 2014
@@ -1,7 +1,46 @@
 # Script for optimising the rotor frame order test model of CaM.
 
 # Python module imports.
-from numpy import array
+from numpy import array, cross, float64, zeros
+from numpy.linalg import norm
+
+# relax module imports.
+from lib.frame_order.rotor_axis import create_rotor_axis_alpha
+from lib.geometry.lines import closest_point_ax
+from lib.geometry.coord_transform import spherical_to_cartesian
+from lib.geometry.rotations import reverse_euler_zyz
+from lib.geometry.vectors import vector_angle
+from pipe_control.structure.mass import pipe_centre_of_mass
+
+
+def alpha_angle(pivot=None, com=None, axis=None):
+    """Calculate and return the rotor alpha angle."""
+
+    # The CoM-pivot axis.
+    com_piv = com - pivot
+    com_piv /= norm(com_piv)
+
+    # The mu_xy vector.
+    z_axis = array([0, 0, 1], float64)
+    mu_xy = cross(z_axis, com_piv)
+    mu_xy /= norm(mu_xy)
+
+    # The alpha angle.
+    return vector_angle(mu_xy, axis, com_piv)
+
+
+def shift_pivot(pivot_orig=None, com=None, axis=None):
+    """Shift the pivot to the closest point on the rotor axis to the CoM.)"""
+
+    # The closest point.
+    pivot_new = closest_point_ax(line_pt=pivot_orig, axis=axis, point=com)
+
+    # Printout.
+    print("\n%-20s%s" % ("Original pivot:", pivot_orig))
+    print("%-20s%s" % ("New pivot:", pivot_new))
+
+    # Return the shifted pivot.
+    return pivot_new
 
 
 # The real parameter values.
@@ -10,6 +49,11 @@
 AXIS_THETA = 0.9600799785953431
 AXIS_PHI = 4.0322755062196229
 CONE_SIGMA_MAX = 30.0 / 360.0 * 2.0 * pi
+
+# Reconstruct the rotation axis.
+AXIS = zeros(3, float64)
+spherical_to_cartesian([1, AXIS_THETA, AXIS_PHI], AXIS)
+print("Rotation axis: %s" % AXIS)
 
 # Create the data pipe.
 pipe.create(pipe_name='frame order', pipe_type='frame order')
@@ -71,7 +115,7 @@
 frame_order.ref_domain('N')
 
 # Set the initial pivot point.
-pivot = array([ 37.254, 0.5, 16.7465])
+pivot = shift_pivot(pivot_orig=array([ 37.254, 0.5, 16.7465]), 
com=pipe_centre_of_mass(verbosity=0), axis=AXIS)
 frame_order.pivot(pivot, fix=True)
 
 # Set the paramagnetic centre.
@@ -87,8 +131,7 @@
 value.set(param='ave_pos_alpha', val=AVE_POS_ALPHA)
 value.set(param='ave_pos_beta', val=AVE_POS_BETA)
 value.set(param='ave_pos_gamma', val=AVE_POS_GAMMA)
-value.set(param='axis_theta', val=AXIS_THETA)
-value.set(param='axis_phi', val=AXIS_PHI)
+value.set(param='axis_alpha', val=alpha_angle(pivot=pivot, 
com=pipe_centre_of_mass(verbosity=0), axis=AXIS))
 value.set(param='cone_sigma_max', val=CONE_SIGMA_MAX)
 calc()
 
@@ -105,6 +148,17 @@
 # Optimise the pivot and model.
 frame_order.pivot(pivot, fix=False)
 minimise('simplex')
+
+# The distance from the optimised pivot and the rotation axis.
+opt_piv = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z])
+print("\n\nOptimised pivot displacement: %s" % norm(pivot - opt_piv))
+pt = closest_point_ax(line_pt=pivot, axis=AXIS, point=opt_piv)
+print("Distance from axis: %s\n" % norm(pt - opt_piv))
+
+# Recreate the axis and compare to the original.
+opt_axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=opt_piv, 
point=pipe_centre_of_mass(verbosity=0))
+print("Original axis:   %s" % AXIS)
+print("Optimised axis:  %s" % opt_axis)
 
 # Test Monte Carlo simulations.
 monte_carlo.setup(number=5)