Author: bugman
Date: Tue Aug 3 14:44:56 2010
New Revision: 11389
URL: http://svn.gna.org/viewcvs/relax?rev=11389&view=rev
Log:
Created the pseudo-ellipse, free rotor target function.
This model is not yet available.
Modified:
1.3/maths_fns/frame_order.py
1.3/maths_fns/frame_order_matrix_ops.py
Modified: 1.3/maths_fns/frame_order.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/maths_fns/frame_order.py?rev=11389&r1=11388&r2=11389&view=diff
==============================================================================
--- 1.3/maths_fns/frame_order.py (original)
+++ 1.3/maths_fns/frame_order.py Tue Aug 3 14:44:56 2010
@@ -31,7 +31,7 @@
from generic_fns.frame_order import print_frame_order_2nd_degree
from maths_fns.alignment_tensor import to_5D, to_tensor
from maths_fns.chi2 import chi2
-from maths_fns.frame_order_matrix_ops import compile_2nd_matrix_iso_cone,
compile_2nd_matrix_pseudo_ellipse,
compile_2nd_matrix_pseudo_ellipse_torsionless, reduce_alignment_tensor
+from maths_fns.frame_order_matrix_ops import compile_2nd_matrix_iso_cone,
compile_2nd_matrix_pseudo_ellipse,
compile_2nd_matrix_pseudo_ellipse_free_rotor,
compile_2nd_matrix_pseudo_ellipse_torsionless, reduce_alignment_tensor
from maths_fns.rotation_matrix import euler_to_R_zyz as euler_to_R
from relax_errors import RelaxError
@@ -330,10 +330,10 @@
return chi2(self.red_tensors, self.red_tensors_bc, self.red_errors)
- def func_pseudo_ellipse_torsionless(self, params):
- """Target function for torsionless pseudo-elliptic cone model
optimisation using alignment tensors.
-
- @param params: The vector of parameter values {alpha, beta, gamma,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the
first 3 are the average position rotation Euler angles, the next 3 are the
Euler angles defining the eigenframe, and the last 2 are the torsionless
pseudo-elliptic cone geometric parameters.
+ def func_pseudo_ellipse_free_rotor(self, params):
+ """Target function for free_rotor pseudo-elliptic cone model
optimisation using alignment tensors.
+
+ @param params: The vector of parameter values {alpha, beta, gamma,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the
first 3 are the average position rotation Euler angles, the next 3 are the
Euler angles defining the eigenframe, and the last 2 are the free_rotor
pseudo-elliptic cone geometric parameters.
@type params: list of float
@return: The chi-squared or SSE value.
@rtype: float
@@ -343,7 +343,7 @@
alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma,
cone_theta_x, cone_theta_y = params
# Generate the 2nd degree Frame Order super matrix.
- frame_order_2nd =
compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, self.rot,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_)
+ frame_order_2nd =
compile_2nd_matrix_pseudo_ellipse_free_rotor(self.frame_order_2nd, self.rot,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_)
# Average position rotation.
euler_to_R(alpha, beta, gamma, self.rot)
@@ -373,3 +373,48 @@
# Return the chi-squared value.
return chi2(self.red_tensors, self.red_tensors_bc, self.red_errors)
+
+
+ def func_pseudo_ellipse_torsionless(self, params):
+ """Target function for torsionless pseudo-elliptic cone model
optimisation using alignment tensors.
+
+ @param params: The vector of parameter values {alpha, beta, gamma,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_y} where the
first 3 are the average position rotation Euler angles, the next 3 are the
Euler angles defining the eigenframe, and the last 2 are the torsionless
pseudo-elliptic cone geometric parameters.
+ @type params: list of float
+ @return: The chi-squared or SSE value.
+ @rtype: float
+ """
+
+ # Unpack the parameters.
+ alpha, beta, gamma, eigen_alpha, eigen_beta, eigen_gamma,
cone_theta_x, cone_theta_y = params
+
+ # Generate the 2nd degree Frame Order super matrix.
+ frame_order_2nd =
compile_2nd_matrix_pseudo_ellipse_torsionless(self.frame_order_2nd, self.rot,
eigen_alpha, eigen_beta, eigen_gamma, cone_theta_x, cone_theta_)
+
+ # Average position rotation.
+ euler_to_R(alpha, beta, gamma, self.rot)
+
+ # Back calculate the reduced tensors.
+ for i in range(self.num_tensors):
+ # Tensor indices.
+ index1 = i*5
+ index2 = i*5+5
+
+ # Reduce the tensor.
+ reduce_alignment_tensor(frame_order_2nd,
self.full_tensors[index1:index2], self.red_tensors_bc[index1:index2])
+
+ # Convert the tensor to 3D, rank-2 form.
+ to_tensor(self.tensor_3D, self.red_tensors_bc[index1:index2])
+
+ # Rotate the tensor (normal R.X.RT rotation).
+ if self.full_in_ref_frame[i]:
+ tensor_3D = dot(self.rot, dot(self.tensor_3D,
transpose(self.rot)))
+
+ # Rotate the tensor (inverse RT.X.R rotation).
+ else:
+ tensor_3D = dot(transpose(self.rot), dot(self.tensor_3D,
self.rot))
+
+ # Convert the tensor back to 5D, rank-1 form.
+ to_5D(self.red_tensors_bc[index1:index2], tensor_3D)
+
+ # Return the chi-squared value.
+ return chi2(self.red_tensors, self.red_tensors_bc, self.red_errors)
Modified: 1.3/maths_fns/frame_order_matrix_ops.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/maths_fns/frame_order_matrix_ops.py?rev=11389&r1=11388&r2=11389&view=diff
==============================================================================
--- 1.3/maths_fns/frame_order_matrix_ops.py (original)
+++ 1.3/maths_fns/frame_order_matrix_ops.py Tue Aug 3 14:44:56 2010
@@ -173,8 +173,8 @@
return matrix
-def compile_2nd_matrix_pseudo_ellipse_torsionless(matrix, R, eigen_alpha,
eigen_beta, eigen_gamma, theta_x, theta_y):
- """Generate the 2nd degree Frame Order matrix for the torsionless pseudo
ellipse.
+def compile_2nd_matrix_pseudo_ellipse_free_rotor(matrix, R, eigen_alpha,
eigen_beta, eigen_gamma, theta_x, theta_y):
+ """Generate the 2nd degree Frame Order matrix for the free rotor pseudo
ellipse.
@param matrix: The Frame Order matrix, 2nd degree to be populated.
@type matrix: numpy 9D, rank-2 array
@@ -193,6 +193,61 @@
"""
# The surface area normalisation factor.
+ fact = 0.125 / pec(theta_x, theta_y)
+
+ # Diagonal.
+ matrix[0, 0] = fact * quad(part_int_daeg2_pseudo_ellipse_free_rotor_00,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
+ matrix[1, 1] = matrix[3, 3] = fact *
quad(part_int_daeg2_pseudo_ellipse_free_rotor_11, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0]
+ matrix[4, 4] = fact * quad(part_int_daeg2_pseudo_ellipse_free_rotor_44,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
+ matrix[8, 8] = fact * quad(part_int_daeg2_pseudo_ellipse_free_rotor_88,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
+
+ # Off diagonal set 1.
+ matrix[0, 4] = matrix[4, 0] = matrix[0, 0]
+ matrix[0, 8] = matrix[8, 0] = fact *
quad(part_int_daeg2_pseudo_ellipse_free_rotor_08, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0]
+ matrix[4, 8] = matrix[8, 4] = fact *
quad(part_int_daeg2_pseudo_ellipse_free_rotor_48, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0]
+
+ # Off diagonal set 2.
+ matrix[1, 3] = matrix[3, 1] = -matrix[1, 1]
+
+ # Average position rotation.
+ euler_to_R_zyz(eigen_alpha, eigen_beta, eigen_gamma, R)
+
+ # The outer product of R.
+ R_kron = kron_prod(R, R)
+
+ # Perform the T23 transpose to obtain the Kronecker product matrix!
+ transpose_23(matrix)
+
+ # Rotate.
+ matrix = dot(R_kron, dot(matrix, transpose(R_kron)))
+
+ # Perform T23 again to return back.
+ transpose_23(matrix)
+
+ # Return the matrix.
+ return matrix
+
+
+def compile_2nd_matrix_pseudo_ellipse_torsionless(matrix, R, eigen_alpha,
eigen_beta, eigen_gamma, theta_x, theta_y):
+ """Generate the 2nd degree Frame Order matrix for the torsionless pseudo
ellipse.
+
+ @param matrix: The Frame Order matrix, 2nd degree to be populated.
+ @type matrix: numpy 9D, rank-2 array
+ @param R: The rotation matrix to be populated.
+ @type R: numpy 3D, rank-2 array
+ @param eigen_alpha: The eigenframe rotation alpha Euler angle.
+ @type eigen_alpha: float
+ @param eigen_beta: The eigenframe rotation beta Euler angle.
+ @type eigen_beta: float
+ @param eigen_gamma: The eigenframe rotation gamma Euler angle.
+ @type eigen_gamma: float
+ @param theta_x: The cone opening angle along x.
+ @type theta_x: float
+ @param theta_y: The cone opening angle along y.
+ @type theta_y: float
+ """
+
+ # The surface area normalisation factor.
fact = 0.25 / pec(theta_x, theta_y)
# Diagonal.
@@ -723,6 +778,124 @@
return 2.0/3.0 * smax * (1.0 - cos(tmax)**3)
+def part_int_daeg2_pseudo_ellipse_free_rotor_00(phi, x, y):
+ """The theta-sigma partial integral of the 2nd degree Frame Order matrix
element 11 for the free rotor pseudo ellipse.
+
+ @param phi: The azimuthal tilt-torsion angle.
+ @type phi: float
+ @param x: The cone opening angle along x.
+ @type x: float
+ @param y: The cone opening angle along y.
+ @type y: float
+ @return: The theta-sigma partial integral.
+ @rtype: float
+ """
+
+ # Theta max.
+ tmax = tmax_pseudo_ellipse(phi, x, y)
+
+ # The theta-sigma integral.
+ return cos(phi)**2*sin(2.0*tmax) + (4.0 - 2.0*cos(phi)**2)*tmax
+
+
+def part_int_daeg2_pseudo_ellipse_free_rotor_08(phi, x, y):
+ """The theta-sigma partial integral of the 2nd degree Frame Order matrix
element 19 for the free rotor pseudo ellipse.
+
+ @param phi: The azimuthal tilt-torsion angle.
+ @type phi: float
+ @param x: The cone opening angle along x.
+ @type x: float
+ @param y: The cone opening angle along y.
+ @type y: float
+ @return: The theta-sigma partial integral.
+ @rtype: float
+ """
+
+ # Theta max.
+ tmax = tmax_pseudo_ellipse(phi, x, y)
+
+ # The theta-sigma integral.
+ return 4.0*cos(phi)**2*tmax - 2.0*cos(phi)**2*sin(2.0*tmax)
+
+
+def part_int_daeg2_pseudo_ellipse_free_rotor_11(phi, x, y):
+ """The theta-sigma partial integral of the 2nd degree Frame Order matrix
element 22 for the free rotor pseudo ellipse.
+
+ @param phi: The azimuthal tilt-torsion angle.
+ @type phi: float
+ @param x: The cone opening angle along x.
+ @type x: float
+ @param y: The cone opening angle along y.
+ @type y: float
+ @return: The theta-sigma partial integral.
+ @rtype: float
+ """
+
+ # Theta max.
+ tmax = tmax_pseudo_ellipse(phi, x, y)
+
+ # The theta-sigma integral.
+ return 4.0*sin(tmax)
+
+
+def part_int_daeg2_pseudo_ellipse_free_rotor_44(phi, x, y):
+ """The theta-sigma partial integral of the 2nd degree Frame Order matrix
element 55 for the free rotor pseudo ellipse.
+
+ @param phi: The azimuthal tilt-torsion angle.
+ @type phi: float
+ @param x: The cone opening angle along x.
+ @type x: float
+ @param y: The cone opening angle along y.
+ @type y: float
+ @return: The theta-sigma partial integral.
+ @rtype: float
+ """
+
+ # Theta max.
+ tmax = tmax_pseudo_ellipse(phi, x, y)
+
+ # The theta-sigma integral.
+ return sin(phi)**2*sin(2.0*tmax) + (4.0 - 2.0*sin(phi)**2)*tmax
+
+
+def part_int_daeg2_pseudo_ellipse_free_rotor_48(phi, x, y):
+ """The theta-sigma partial integral of the 2nd degree Frame Order matrix
element 59 for the free rotor pseudo ellipse.
+
+ @param phi: The azimuthal tilt-torsion angle.
+ @type phi: float
+ @param x: The cone opening angle along x.
+ @type x: float
+ @param y: The cone opening angle along y.
+ @type y: float
+ @return: The theta-sigma partial integral.
+ @rtype: float
+ """
+
+ # Theta max.
+ tmax = tmax_pseudo_ellipse(phi, x, y)
+
+ # The theta-sigma integral.
+ return 4.0*sin(phi)**2*tmax - 2.0*sin(phi)**2*sin(2.0*tmax)
+
+
+def part_int_daeg2_pseudo_ellipse_free_rotor_88(phi, x, y):
+ """The theta-sigma partial integral of the 2nd degree Frame Order matrix
element 99 for the free rotor pseudo ellipse.
+
+ @param phi: The azimuthal tilt-torsion angle.
+ @type phi: float
+ @param x: The cone opening angle along x.
+ @type x: float
+ @param y: The cone opening angle along y.
+ @type y: float
+ @return: The theta-sigma partial integral.
+ @rtype: float
+ """
+
+ # Theta max.
+ tmax = tmax_pseudo_ellipse(phi, x, y)
+
+ # The theta-sigma integral.
+ return 2.0*sin(2.0*tmax) + 4.0*tmax
def part_int_daeg2_pseudo_ellipse_torsionless_00(phi, x, y):
"""The theta partial integral of the 2nd degree Frame Order matrix
element 11 for the torsionless pseudo ellipse.
r11389 - in /1.3/maths_fns: frame_order.py frame_order_matrix_ops.py