r26427 - /branches/frame_order_cleanup/lib/frame_order/pseudo_ellipse.py -- November 06, 2014 - 15:48

Author: bugman
Date: Thu Nov  6 15:48:03 2014
New Revision: 26427

URL: http://svn.gna.org/viewcvs/relax?rev=26427&view=rev
Log:
Simplifications and fixes for the 1st degree frame order matrix calculation 
for the pseudo-ellipse.

The compile_1st_matrix_pseudo_ellipse() function of the 
lib.frame_order.pseudo_ellipse module has
been significantly simplified by shifting a lot of maths outside of the 
quadratic integration.


Modified:
    branches/frame_order_cleanup/lib/frame_order/pseudo_ellipse.py

Modified: branches/frame_order_cleanup/lib/frame_order/pseudo_ellipse.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/lib/frame_order/pseudo_ellipse.py?rev=26427&r1=26426&r2=26427&view=diff
==============================================================================
--- branches/frame_order_cleanup/lib/frame_order/pseudo_ellipse.py      
(original)
+++ branches/frame_order_cleanup/lib/frame_order/pseudo_ellipse.py      Thu 
Nov  6 15:48:03 2014
@@ -54,7 +54,7 @@
     """
 
     # The surface area normalisation factor.
-    fact = 2.0 * sigma_max * pec(theta_x, theta_y)
+    fact = 2.0 * pec(theta_x, theta_y)
 
     # Invert.
     if fact == 0.0:
@@ -62,10 +62,13 @@
     else:
         fact = 1.0 / fact
 
+    # Sinc value.
+    sinc_smax = sinc(sigma_max/pi)
+
     # Numerical integration of phi of each element.
-    matrix[0, 0] = fact * quad(part_int_daeg1_pseudo_ellipse_xx, -pi, pi, 
args=(theta_x, theta_y, sigma_max), full_output=1)[0]
-    matrix[1, 1] = fact * quad(part_int_daeg1_pseudo_ellipse_yy, -pi, pi, 
args=(theta_x, theta_y, sigma_max), full_output=1)[0]
-    matrix[2, 2] = fact * quad(part_int_daeg1_pseudo_ellipse_zz, -pi, pi, 
args=(theta_x, theta_y, sigma_max), full_output=1)[0]
+    matrix[0, 0] = fact * sinc_smax * (2.0*pi + 
quad(part_int_daeg1_pseudo_ellipse_xx, -pi, pi, args=(theta_x, theta_y, 
sigma_max), full_output=1)[0])
+    matrix[1, 1] = fact * sinc_smax * (2.0*pi + 
quad(part_int_daeg1_pseudo_ellipse_yy, -pi, pi, args=(theta_x, theta_y, 
sigma_max), full_output=1)[0])
+    matrix[2, 2] = 0.5 * fact * quad(part_int_daeg1_pseudo_ellipse_zz, -pi, 
pi, args=(theta_x, theta_y, sigma_max), full_output=1)[0]
 
     # Rotate and return the frame order matrix.
     return rotate_daeg(matrix, R_eigen)
@@ -152,7 +155,7 @@
     tmax = tmax_pseudo_ellipse(phi, x, y)
 
     # The theta-sigma integral.
-    return sin(smax) * (2 * (1 - cos(tmax)) * sin(phi)**2 + cos(phi)**2 * 
sin(tmax)**2)
+    return cos(phi)**2 * sin(tmax)**2  -  2.0 * sin(phi)**2 * cos(tmax)
 
 
 def part_int_daeg1_pseudo_ellipse_yy(phi, x, y, smax):
@@ -174,7 +177,7 @@
     tmax = tmax_pseudo_ellipse(phi, x, y)
 
     # The theta-sigma integral.
-    return sin(smax) * (2.0 * cos(phi)**2 * (1.0 - cos(tmax)) + sin(phi)**2 
* sin(tmax)**2)
+    return sin(phi)**2 * sin(tmax)**2  -  2.0 * cos(phi)**2 * cos(tmax)
 
 
 def part_int_daeg1_pseudo_ellipse_zz(phi, x, y, smax):
@@ -196,7 +199,7 @@
     tmax = tmax_pseudo_ellipse(phi, x, y)
 
     # The theta-sigma integral.
-    return smax * sin(tmax)**2
+    return sin(tmax)**2
 
 
 def part_int_daeg2_pseudo_ellipse_00(phi, x, y, smax):