Author: bugman
Date: Wed Aug 4 22:42:48 2010
New Revision: 11414
URL: http://svn.gna.org/viewcvs/relax?rev=11414&view=rev
Log:
Fixes for the torsionless frame order matrix equations.
The previous equations were plain wrong!
Modified:
1.3/maths_fns/frame_order_matrix_ops.py
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=11414&r1=11413&r2=11414&view=diff
==============================================================================
--- 1.3/maths_fns/frame_order_matrix_ops.py (original)
+++ 1.3/maths_fns/frame_order_matrix_ops.py Wed Aug 4 22:42:48 2010
@@ -254,23 +254,23 @@
"""
# The surface area normalisation factor.
- fact = 0.25 / pec(theta_x, theta_y)
+ fact = 1.0 / (6.0 * pec(theta_x, theta_y))
# Diagonal.
- matrix[0, 0] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_00,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
- matrix[1, 1] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_11,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
- matrix[2, 2] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_22,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
+ matrix[0, 0] = fact * (6.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_00, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
+ matrix[1, 1] = fact * (2.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_11, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
+ matrix[2, 2] = fact * (5.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_22, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
matrix[3, 3] = matrix[1, 1]
- matrix[4, 4] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_44,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
- matrix[5, 5] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_55,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
+ matrix[4, 4] = fact * (6.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_44, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
+ matrix[5, 5] = fact * (5.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_55, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
matrix[6, 6] = matrix[2, 2]
matrix[7, 7] = matrix[5, 5]
matrix[8, 8] = fact * quad(part_int_daeg2_pseudo_ellipse_torsionless_88,
-pi, pi, args=(theta_x, theta_y), full_output=1)[0]
# Off diagonal set 1.
- matrix[0, 4] = matrix[4, 0] = fact *
quad(part_int_daeg2_pseudo_ellipse_torsionless_04, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0]
- matrix[0, 8] = matrix[8, 0] = fact *
quad(part_int_daeg2_pseudo_ellipse_torsionless_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_torsionless_48, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0]
+ matrix[0, 4] = matrix[4, 0] = fact * (2.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_04, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
+ matrix[0, 8] = matrix[8, 0] = fact * (4.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_08, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
+ matrix[4, 8] = matrix[8, 4] = fact * (4.0*pi +
quad(part_int_daeg2_pseudo_ellipse_torsionless_48, -pi, pi, args=(theta_x,
theta_y), full_output=1)[0])
# Off diagonal set 2.
matrix[1, 3] = matrix[3, 1] = matrix[0, 4]
@@ -936,6 +936,8 @@
# 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 for
the torsionless pseudo-ellipse.
@@ -953,7 +955,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return (cos(phi)**4*sin(2.0*tmax) +
8.0*cos(phi)**2*sin(phi)**2*sin(tmax) + (4.0*sin(phi)**4 +
2.0*cos(phi)**4)*tmax)
+ return (2*cos(phi)**4*cos(tmax) +
6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - (6*sin(phi)**4 +
2*cos(phi)**4)*cos(tmax)
def part_int_daeg2_pseudo_ellipse_torsionless_04(phi, x, y):
@@ -973,7 +975,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return (cos(phi)**2*sin(phi)**2*sin(2.0*tmax) -
8.0*cos(phi)**2*sin(phi)**2*sin(tmax) + 6*cos(phi)**2*sin(phi)**2*tmax)
+ return (2*cos(phi)**2*sin(phi)**2*cos(tmax) -
6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - 8*cos(phi)**2*sin(phi)**2*cos(tmax)
def part_int_daeg2_pseudo_ellipse_torsionless_08(phi, x, y):
@@ -993,7 +995,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return - (cos(phi)**2*sin(2.0*tmax) - 2.0*cos(phi)**2*tmax)
+ return 2*cos(phi)**2*cos(tmax)**3 - 6*cos(phi)**2*cos(tmax)
def part_int_daeg2_pseudo_ellipse_torsionless_11(phi, x, y):
@@ -1013,7 +1015,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return (cos(phi)**2*sin(phi)**2*sin(2.0*tmax) + (4.0*sin(phi)**4 +
4.0*cos(phi)**4)*sin(tmax) + 6*cos(phi)**2*sin(phi)**2*tmax)
+ return (2*cos(phi)**2*sin(phi)**2*cos(tmax) + 3*sin(phi)**4 +
3*cos(phi)**4)*sin(tmax)**2 - 8*cos(phi)**2*sin(phi)**2*cos(tmax)
def part_int_daeg2_pseudo_ellipse_torsionless_22(phi, x, y):
@@ -1033,7 +1035,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return (cos(phi)**2*sin(2.0*tmax) + 4.0*sin(phi)**2*sin(tmax) +
2.0*cos(phi)**2*tmax)
+ return (2*sin(phi)**2 - 2)*cos(tmax)**3 - 3*sin(phi)**2*cos(tmax)**2
def part_int_daeg2_pseudo_ellipse_torsionless_44(phi, x, y):
@@ -1053,7 +1055,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return (sin(phi)**4*sin(2.0*tmax) +
8.0*cos(phi)**2*sin(phi)**2*sin(tmax) + (2.0*sin(phi)**4 +
4.0*cos(phi)**4)*tmax)
+ return (2*sin(phi)**4*cos(tmax) +
6*cos(phi)**2*sin(phi)**2)*sin(tmax)**2 - (2*sin(phi)**4 +
6*cos(phi)**4)*cos(tmax)
def part_int_daeg2_pseudo_ellipse_torsionless_48(phi, x, y):
@@ -1073,7 +1075,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return - (sin(phi)**2*sin(2.0*tmax) - 2.0*sin(phi)**2*tmax)
+ return 2*sin(phi)**2*cos(tmax)**3 - 6*sin(phi)**2*cos(tmax)
def part_int_daeg2_pseudo_ellipse_torsionless_55(phi, x, y):
@@ -1093,7 +1095,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return sin(phi)**2*sin(2.0*tmax) + 4.0*cos(phi)**2*sin(tmax) +
2.0*sin(phi)**2*tmax
+ return (2*cos(phi)**2 - 2)*cos(tmax)**3 - 3*cos(phi)**2*cos(tmax)**2
def part_int_daeg2_pseudo_ellipse_torsionless_88(phi, x, y):
@@ -1113,7 +1115,7 @@
tmax = tmax_pseudo_ellipse(phi, x, y)
# The theta integral.
- return sin(2.0*tmax) + 2.0*tmax
+ return 2 - 2*cos(tmax)**3
def populate_1st_eigenframe_iso_cone(matrix, angle):
r11414 - /1.3/maths_fns/frame_order_matrix_ops.py